Resolution matrix

Dear Francois,

I want to generate metrics such as point spread and crosstalk functions to assess inverse solution accuracy, in my case sLORETA.

For that, I need to calculate the resolution matrix (Liu 2002, Pasqual Marqui 2002)

R = W*A; W-inverse operator, A- gain matrix

Could you confirm that using BST,
R = Kernel.Imaging * Headmodel.Gain

Thank you,

Octavian

1: Pascual-Marqui RD. Standardized low-resolution brain electromagnetic
tomography (sLORETA): technical details. Methods Find Exp Clin Pharmacol. 2002;24
Suppl D:5-12. PubMed PMID: 12575463.

2: Liu AK, Dale AM, Belliveau JW. Monte Carlo simulation studies of EEG and MEG
localization accuracy. Hum Brain Mapp. 2002 May;16(1):47-62.

Dear Octavian,

Good question, like usual. Yes, the resolution matrix will be:
R = Kernel.Imaging * Headmodel.Gain;
This matrix maps the space of ‘true’ currents to the space of estimated currents.
The point spread functions will be its column vectors. The ith column (i.e., psf) will be a map of how true activity in the ith dipole component spreads from that source point to all over the source space during estimation.
The resolution kernels will be the row vectors. These will be images of how activity from all over the source space contributes to the ith component estimate.

Now, keep in mind that in BST, we usually recommend for people to whitened the data. And if the data has been processed with SSS, SSP, or any other dimensionality reduction denoising algorithm, the noise covariance will not be full rank.
That means that the whitening operator will not be square. It will have more columns than rows. Thus, the whitened data will have less signals than the original data. This will also be true about the whitened lead field matrix. And it will also imply that the inverse operator will technically have less column vectors than there are sensors. BUT, in BST, we have taken the convention of pre-multiplying the inverse operator from the right side with the whitening operator. This is why the ImagingKernel matrix will have as many columns as sensors. Thus, when this imaging kernel is multiplied with the unwhitened data, it automatically whitens the data, and estimates the sLORETA activity. This will also be true of when you multiply the imaging kernel with the unwhitened lead field matrix, thereby whitening it.

I wrote this function to compute the sLORETA estimate to be equivalent to the Pascual-Marqui paper. But I should make a few comments to clarify details. First of all, we are using whitening to model and suppress correlated noise. Also, we are computing the MNE part of the sLORETA the same way that the MNE Suite does it. This means that the SNR2 can be modeled as the ratio of the whitened data power, taken into account the fact that a dimensionality reduction may have occurred. This implies that the regularization parameter lambda2=1/SNR2. Also, it implies that the source covariance matrix has been adjusted as a function of the trace of LC_sL’, so that the trace is equal to the rank of the noise covariance, where L is the lead field matrix and C_s is the source covariance, i.e., a multiple of the identity matrix. I recommend checking out the MNE Suite Manual for more details. Also, check at the matlab code. It’s well commented.

Also, I should mention that all the sLOERETA computations related to the inverse block resolution matrices have been taken care of in a for loop in such a way that when you multiply your data or your lead field matrix with the ImagingKernel, it is directly producing the sLORETA activity. In the case of fixed orientation constraints, taking the absolute value of the activity will give you Pascual-Marqui’s absolute standardized activity. But you don’t have to take that absolute value if you want to for example look at non-rectified oscillations in the source domain. In the case, of free dipole orientations, then you have 3 components to the sLORETA activity per source point. If you square them and sum them, you have the sLORETA power, and if you further take the square root you have a modulus of sLORETA activity. These will be rectified, so don’t do frequency analysis at this step. Do it at the sensor level, or at the sLORETA component level, before taking the power.

In BST, we allow sLORETA solutions to have loose orientation constraints, in addition to fixed and free. Personally, I would recommend to use free or loose orientation constraints. In both cases, you will have three components per source point. For statistical contrasts, the ERFs can be subtracted so that the sLORETA reflects a difference, for example, or likewise you can subtract the sLORETA activity components from the two conditions. I believe that in BST the effective number of averages is taken care of by dividing the noise covariance appropriately, taking into account the number of trials in the two conditions. If not, or to be sure, I would take the difference at the source level of sLORETA activity per component, and then take the modulus of the difference sLORETA map, before moving the group analyses. Contrasts in the sensor or source domain should give the same exact results, literally exactly the same. Of course, there are other statistical hypotheses that can be made, actually many many versions, and it all depends on the questions being asked, e.g., statistical differences within subject averaging across epochs, versus group statistics relative to pre-stimulus baseline, or relative to a different conditions, etc, etc, etc. It gets complicated. Just mentioning a few cases for people to be aware of some of the possibilities.

Please don’t hesitate to ask more questions, and I will try to answer as best as possible.

Cheers,

Rey

Dear Rey,

Thank you for your very prompt and detailed aswer.
2 related questions:

  1. For one specific application (regularization by crossvalidation, l curve), I need to get the actual whitened (measured) data; can I assume that
    whiteneddatamat.F = kernelmat.Whitener * datamat.F;?

  2. Is there a way/code already written in order to visualize in BST or matlab the sLORETA source orientation as a sheet/grid of minidipoles with the conventional sphere rod?
    Thank you, Octavian.

Dear Ray,

Speaking about the resolution matrix (R = Kernel.Imaging * Headmodel.Gain) in case of sLORETA, and looking at the bst_wmne code, which Headmodel.Gain should I use: Lk (line 318), L (whitened, line 358) or depth normalized (LW, line 413).
Thank you very much,

Octavian

I assume is the unwhitened leadfield Lk (line318), but please confirm. Another basic question related to the resolution matrix R = Kernel.Imaging * Headmodel.Gain: since I use sLORETA unconstrained (numdipcomp = 3), when calculating the point spread or cross-talk metric for a location i, should I sum the squares the 3 components per location i for both Kernel.Imaging and Headmodel.Gain, and thus compress the resolution matrix, or should just go ahead with multiplication. In the former case, I would lose all orientation info, so it does not make sense.
(I am trying to calculate Liu’s crosstalk metric for a location i due to activity from location j: Eij^2= [(WA)ij]^2/[(WA)ii]^2 , where W-inverse operator, A- gain matrix).

Please advise,

Octavian

Hi Octavian,

Good questions, like always. Here I answer a bunch of them that have accumulated.

  1. Yes, “whiteneddatamat.F = kernelmat.Whitener * datamat.F;” is the whitened data.

  2. I don’t think BST can visualize dipole orientations. This can be programmed easily though, and could be included in BST in the future.

  3. About which Leadfield Matrix to use … definitely use the unwhitened and unweighted, remember what I said in the last email, the inverse operator has been pre-multiplied with the whitenining operator from the right side. So, there’s no need to whitened the leadfied matrix.

  4. Also, it’s no coincidence that for sLORETA, BST’s default parameters for depth weighting are NO depth weighting. Actually, it is very strongly discouraged to use any depth weighting for sLORETA. The necessary depth weighting is done with the resolution matrix indeed, and if you try to use any depth weighting other than a multiple of the identity matrix as a source covariance matrix, you will end up with localization bias. This is of course not true for wMNE and dSPM, which always have localization bias, and do need a depth compensation prior source covariance matrix.

  5. About the psf and cross-talk metrics for the case of free or loose orientation, you could sum the 3 squares (a power measure) or also take the square root of the sum of the squares. As long as you specify what you did it’s OK. If you want to replicate Liu’s work and metrics, you should do exactly what he did.

Cheers,

Rey