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
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)
- 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
- 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