Compute EEG sources with sLoreta

Discussions
,
1

I have questions about the computation of the inverse solution (imaging Kernel) of EEG data using sLoreta. My explanation is long, but I have tried to be as clear as possible. In the explanation below, I am referencing to sLoreta implementation found in bst_inverse_linear_2018.m. My questions are at the end of this long explanation.

Let’s assume that M are raw EEG data recorded at Ne electrodes (with respect to an arbitrary reference electrode). Let Cn be the noise covariance matrix of the sensors (Cn has size Ne x Ne); for example let Cn be the identity matrix if we assume ‘no noise modeling’ (actually brainstorm sets the diagonal elements at 10^-10 as bst assumes noise = 10 microV). We want to estimate the activity of Nd dipoles assuming that they have fixed orientation orthogonal to the cortex. Let A be the lead field matrix having size Ne x Nd (i.e. the lead field matrix A already includes dipole orientation).

If I understand rightly, brainstorm applies the average reference (let’s say V, having size Ne x Ne) to both the data and the lead field matrix. In the following, I will call Mv (=V x M) the re-referenced EEG data, Av (=V x A) the re-referenced lead field matrix, Cnv (= V x Cn x V’) the covariance matrix of the re-referenced noise, where the symbol ‘ indicate the transpose.
In other words, we have:
V x M = (V x A) x S + V x n that is Mv = Av x S + nv
where M is the raw data matrix (size Ne x T time instants), S the source activities (size Nd x T) that we want to estimate, and n the raw noise matrix (size Ne x T).

Furthermore, if I understand rightly, brainstorm assumes to apply the inverse whitener operator iW, defined as:
iW = Cnv^(-1/2) (i.e iW is the covariance matrix Cnv elevated by -1/2)
So we have
iW x Mv = iW x Av x S + iW x nv that is Mvw = Avw x S + nvw
where Mvw = iW x Mv, Avw= iW x Av, nvw= iW x nv.

For consistency with bst_inverse_linear_2018.m, I will set Avw=L; accordingly we have
Mvw = L x S + nvw

Now, let’s compute the covariance matrix (Cd) of the whitened and re-referenced data Mvw:
Cd = L x Cs x L’ + I
where Cs is the “a priori” covariance matrix of the sources and I is the identity matrix (Ne x Ne) representing the covariance matrix of the whitened noise nvw.

The Cs is assumed equals to lambda x I and the scalar lambda is computed to achieve the desired SNR ratio. So we have,
Cd = lambda x L x L’ + I

In the minimum norm estimation, the inverse operator K is computed as
Kmne = Cs x L’ x (L x Cs x L’ + I)^-1 = lambda x L’ x (lambda x L x L’ + I)^-1 = lambda x L’ x (Cd)^-1.
Indeed brainstorm applies the above equation to compute the inverse operator Kmne (exploiting singular value decomposition to compute the inverse of Cd). So up to now – if there are no mistakes in my reasoning – I have no specific doubt.

However, to compute sLoreta solution, the above kernel Kmne must be modified since sLoreta standardizes the estimated sources by their standard deviation (accounting both for the a priori source variation and measurement noise). Precisely, according to the work of Pascual-Marqui 2002, we should have:
Se = Kmne x Mvw,
where Se are the estimated sources using the minimum norm operator.

Hence, the covariance matrix Cse of the estimated sources is
Cse = Kmne x Cd x Kmne’.

Therefore, based on the expression of Kmne, we have
Cse = lambda x L’ x ((Cd)^-1) x Cd x (lambda x L’ x (Cd)^-1)’= lambda x L’ x ((Cd)^-1) x lambda x L

That is
Cse = Kmne x lambda x L.

Then, the imaging kernel for sLoreta source estimation should be
Kslor = bsxfun(@rdivide, Kmne, sqrt(diag(Cse)) (Eq. a)

  1. My very first question is that it seems that in brainstorm Cse is computed as Cse=Kmne x L rather than Cse=lambda x Kmne x L as I have obtained in my computation. Is there some error in my computation or something that I have missed?

Finally, according to the previous explanation, the sLoreta image should obtained as
S(sLoreta)=Kslor x Mvw = Kslor x iW x V x M

Hence, to apply the imaging kernel directly to the raw data M, I would provide as output the kernel Kslor obtained by (Eq. a) multiplied by iW x V, i.e.
Kslor = Kslor x iW x V

  1. My second question is that the final imaging kernel provided by brainstorm (to be applied to the raw M data), is obtained multiplying only by the inverse whitener iW rather than by iW x V. Could you explain where I am wrong?

Thanks a lot in advance

2

@rey @John_Mosher @Sylvain?

3

@rey @John_Mosher @Sylvain @tmedani @juangpc ?

4

Hello Francois, nothing new about my questions?

5

Unfortunately, I'm not the right person to answer them.
I was hoping for a reply from @rey, @John_Mosher or @Sylvain ...

6

% Here's my two cents ElisaM.
% Sorry this took so long. I'm just very busy, and I haven't been working
% on sLORETA for a while, even though I did write the sLORETA code for BST.
% I now use eLORETA (not sLORETA) for simple distributed linear imaging instead,
% and I'm working on COMPLETELY different methods for inverse imaging.
% This is an actual program written in Matlab tongue, so you can copy paste
% and reproduce everything as many times as you want.

L=randn(10,100); % Toy problem, "lead field" matrix is sampled from a Gaussian distribution, only 10 channels, 100 dipoles.
Ln=randn(10,13); % Add a noise dictionary to get covariance matrix with non-diagonal structure (i.e., correlatted noise).
Jn=randn(13,1000); % The noise current density time series associated with the noise dictionary.
Bn=LnJn; % Simulated noise data.
Cn=cov(Bn')/100; % Noise covariance matrix.
Cn=(Cn+Cn')/2; % Making sure it's perfectly symmetric.
W=inv(sqrtm(Cn)); % Whitening matrix.
LW=WL; % Whitened lead field.
Sigma_J=eye(100,100)(10/trace(LWLW')); % A priori source covariance.

% Numerically more stable SVD approach. I think this is how I wrote for
% Brainstorm many years ago.
Rc=chol(Sigma_J,'lower'); % Cholesky used only in case we had nonzero off-diagonal elements, which we usually never do, at least in BST.
LWRc=LWRc; % This is the matrix we need to SVD.
[U,S,V]=svd(LWRc,'econ'); % SVDing.
s=diag(S); % vector of singular values.
SNR=3; % Just a common safe default, could be time-dependent based on whitened data.
lambda2=SNR.^(-2); % Would be time dependent if SNR was time-dependent.
ss=s./(s.^2 + lambda2); % Same as above.
Kernel=RcV*diag(ss)U'; % This is the MNE "imaging kernel".
R=KernelLW; % Resolution matrix for MNE "imaging kernel"
sloretadiag=sqrt(sum(Kernel.*LW',2)); % computing the square root of diagonal of R above. Equivalently, sloretadiag=sqrt(diag(R));
KernelSLO=Kernel./repmat(sloretadiag(:),[1 10]); % sLORETA operator for the simpler case of fixed dipole orientations.

% Gram numerically noisy approach (i.e., inner product might accumulate error when there's a lot of dipoles, but it doesn't really matter in noisy reality)
T=Sigma_JLW'inv(LWSigma_JLW' + lamda2*eye(10,10)); % Exactly equivalent to MNE Kernel via SVD,
% which is the same as:
T_alt=(1/lambda2)Sigma_JLW'inv(LW(1/lambda2)Sigma_JLW' + eye(10,10)); % Exactly equivalent to MNE Kernel via SVD
corrcoef(Kernel,T)
norm(Kernel)
norm(T)
corrcoef(T,T_alt)
norm(T)
norm(T_alt)

T_oldschool=Sigma_JL'inv(LSigma_JL' + lamda2Cn); % Matrix is not identical to Kernel (as it should be) but it will produce the same
% exact current density vector when multiplied by data vector that you get if you multiply Kernel
% by whitened data vector.
% Note that L is used not LW.
jK=KernelLW(:,1); % essentialy the point-spread function for first dipole
j_oldschool=T_oldschool*L(:,1); % essentialy the point-spread function for first dipole
corrcoef(j,jK) % perfect correlation
norm(jK) % same norms
norm(j) % same norms

RT=TLW; % Matches R as it should. Resolution matrix obviously has to since T and Kernel are identical.
R_oldschool=T_oldschoolL; % Also matches R as it should. Note that you have to multiply by L not LW.
corrcoef(R,R_oldschool)
norm(R)
norm(R_oldschool)
% So all these different ways of computing the MNE resolution matrix are going
% to give you the same exact thing. Therefore, they will give you the same
% exact sLORETA 'inverse' operator. The important matrix here again, I
% want to emphasize is the resolution matrxix (as far as I know).

% Now let's look at equation (15) of Pascual-Marqui 2002
S_j_PM=T_oldschool*(LSigma_JL' + lambda2Cn)T_oldschool';
% Is it exactly equal to the resolution matrix? Let's check.
corrcoef(R,S_j_PM) % well you should get a perfect correlation of 1 (just a fast check as we hack along)
% But is there a scale difference?
norm(R)
norm(S_j_PM) % Yes it seems there's a scale factor that is off.
% To get the exact same scale for the resolution matrix.
% You have to do this:
S_j_alt=(1/Sigma_J(1,1))T_oldschool(LSigma_JL' + lambda2*Cn)*T_oldschool'; % Note that the factor is 1/Sigma_J(1,1) not lambda2 as in your email, these are different.
% Let's check correlation
corrcoef(R,S_j_alt)
norm(R)
norm(S_j_alt) % Bingo they match now, which is what I think you had noticed. You are very very veru observant.

% But what's important for sLORETA is the resolution matrix.
% Let's tale a look at the old school resolution matrix.
R_oldschool=T_oldschool*L; % Note that you multiply by L not LW in this case.
corrcoef(R,R_oldschool) % correlation of 1.
norm(R) % checking scales
norm(R_oldschool) % yes they match, and they are identical.

% So my two cents is that as far as I know what really matters for computing sLORETA
% is the MNE resolution matrix, which was done here in several different
% ways and gets the same exact result. You were right that there's a
% scaling issue with equation (15) but that doesn't matter as far as I
% know. You have to use the resolution matrix (see Backus-Gilbert work for
% more info on resoution matrices and their optimization).

% AS A FINAL NOTE, I REALLY WANT TO EMPHASIZE THAT I PERSONALLY STRONGLY
% RECOMMEND ELORETA OVER ANY OTHER LINEAR DISTRIBURED METHOD INCLUDING
% SLORETA. WITH ELORETA YOU GET ZERO LOCALIZATION BIAS AND THEY ARE
% CURRENT DENISTY SOLUTIONS, MEANING REAL FEASIBLE SOLUTIONS, NOT ACTIVITY
% MAPS WHICH IS WHAT SLORETA OR BEAMFORMERS GIVE YOU, YUCK.

% About question #2 I think you have to ask Francois but I think basically
% he probably multiplies the imaging Kernel by the whitener so that the
% data is used instead of the whitened data. If this is the case, it
% actually uses more memory and more floating point operations, but there
% might be a practical reason. And as for the EEG case, well I have a
% feeling that the data is already rereferenced to the average 'electrode' so no need to
% keep doing that. The data is already probably stored properly referenced to the mean.
% Hope it makes sense. Best regards. And keep looking for
% mistakes as there could be some waiting to be corrected.
% In any case, if there had been an error, it would have been a simple
% scaling and that wouldn't matter much for sLORETA as those activity maps
% are in arbitrary units, not units of current densisty.
% These are only my two cents. Anyone should feel more than free to look
% into this some more and possibly correct me if I'm wrong this time.
% It's very tricky business, so the more brains the merrier. Brainstorm
% away. –Rey Rene Ramirez

7

ElisaM, some stuff got corrupted in the code I copy pasted.
I sent it to Francois, maybe he can post it. Otherwise, email me directly, and I'll send you the m-file
of the script. rrramir@uw.edu

8

test_sloreta.m (6.6 KB)

OK, ElisaM, I figured how to upload m files to the BST forum.

9

@rey
Thanks for the thorough response.

What would it take to get eLORETA into Brainstorm?

The current implementation of sLORETA is in this function:
https://github.com/brainstorm-tools/brainstorm3/blob/master/toolbox/inverse/bst_inverse_linear_2018.m#L901
Would it be just an extra case to add in there? Could you share code for it?

Last question:
I've just discussed this topic with Richard and Sylvain, and they asked if you could share some reference mentioning (or write) in simple terms the difference between sLORETA and eLORETA?

10

Yes I can share my Matlab code, of course. It might take some time for me to get it up to you since I'm very busy at the moment. But I can prioritize it, if there's interest in trying it out.
There's also eLORETA code (I think) from MNE and Fieldtrip (Alex had mentioned it a while back).
I haven't tried them or compared them with my code. That could be an interesting project for someone.
The eLORETA algorithm involves "learning" the a priori source covariance given the lead field matrix,
and a given regularization parameter. Thus, there's an extra step that is iterative, a nonlinear optimization. Once that's learned, it is fixed (at least for the regularization parameter and leadfield matrix), and the algorithm works like a regular L2-minimum norm approach.

11

Hi Rey,
thank you so much for your response and the uploaded file. I am just back from a conference and I will go through it in the next days.

Of course, it would be very useful to have eLoreta implemented in Brainstorm, too. Or even if you could share your matlab code.

Thank again for your help

12

@rey
That would be great if you could help us implementing this in Brainstorm, indeed!
Thanks!

@Sylvain @leahy @John_Mosher
Anybody in your labs who could be interested in testing it vs FieldTrip or MNE?

13

@Francois,
Here's the eLORETA code you asked for (attached).
I see that some changes were made on bst_inverse_linear_2028.m.
Hopefully none of those changes involved how the SNR, the noise regularization
parameter, and the whitening are done, as those are inter-related, and would affect this code if different
from what they should be. My advice for the future, is if it's not broken don't try to fix it..
I am attaching an m file that has only the eLORETA specific code for both fixed and free dipole
orientations.
It is assumed that you have whitened the lead field matrix and provided
a fix regularization parameter, but I hint in comments as to how to also use time-dependent SNR.
I have made no attempt to reconnect it with the larger program, but it should be clear
where the code should go (i.e., after you have whitened the lead field matrix and have computed a regularization parameter). The last line gives you the Kernel_eLORETA, which is what you are after and works just like the weighted minumum-L2-norm Kernel.
Note that for eLORETA the source covariance is diagonal but it's not a multiple of the identity matrix.
So some of the "new" edits that I did not introduce treat it as a scalar, and that is just not good, and it's incompatible with any method like eLORETA that does a bit more and actually learns a source a priori covariance that results in absolutely zero localization bias and real current density solutions (not arbitrary unit "activity").
Best,
Rey
eloreta_code_chunk.m (5.0 KB)

14

@John_Mosher (and @juangpc ?)
Could you check if this works with your current implementation of bst_inverse_linear_2018.m?

Thanks

15

I know there is a long "to-do" list for further improvements to Brainstorm. Just thought you might be interested in this new and very exhaustive paper comparing inverse solutions with MEG, including eLORETA. Interestingly, the latter wins on almost every front. Hopefully a reason to promote it higher up the list :wink:

https://www.biorxiv.org/content/10.1101/2020.01.12.903302v1.abstract