Discrete, 3D distributed, linear imaging methods of electric neuronal activity. Part 1: exact, zero error localization

Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Discr et e, 3D distributed line ar imaging methods of electric neur onal activity . Part 1: exact , zer o err or localization Roberto D. Pascual-Marqui The KEY Institute for Brain-Mind Research University Hospital of Psychiatry Lenggstr. 31, CH-8032 Zurich, Switzerland pascualm at key.uzh.ch www.keyinst.uzh.ch/loreta 1. Abstract This paper deals with the EEG/MEG neuroimaging pro blem: given measurements of scalp electric potential differences (EEG: elec troencephalogram) an d extracranial magnetic fields (MEG: magnetoencephalogram), find the 3D distribution of the generating electric neuronal activity. This problem has no unique solution. Only particular solutions with “good” localization properties are of interest , since neuroimaging is concerned with the localization of brain function. In this paper, a general family of linear imaging methods with exact, zero error localization to point-test so urces is presented. One particular member of this family is sLORETA (standardized low resolution brain electromagnetic tomography; Pascual-Marqui, Methods Find. Exp. Clin. Pharmacol. 2002, 24D:5-12; http://www.unizh.ch/keyinst/NewLORETA/sLORETA/sLORETA-Math01.pdf ). It is shown here that sLORETA has no localization bias in the presence of measurem ent and biological noise. Another member of this family, denoted as eLORETA (exact low resolution brain electromagnetic tomography ; Pascual-Marqui 2005: http://www.research-projects.unizh.ch/p6990.htm ), is a genuine inverse solution (not merely a linear imaging method) with exact, zero error localization in the presence of measurement and structured biological noise. The general family of imaging methods is further extended to include da ta-dependent (adaptive) quasi-linear imaging methods, also with the exact, zero error localization property. 2. The forward equation Details on the electrophysiology and physics of EEG/MEG generation can be found in Mitzdorf (1985), Llinas (1988), Martin (1991), Hämäläinen et al (1993), Haalman and Vaadia (1997) , Sukov and Barth (1998), Dale et al (2000) , Baillet et al. (2001). The basic underlying physics can be studied in Sarvas (1987). Consider the forward EEG equation: Eq. 1: c =+ KJ 1 Φ where the vector contains instantaneous scalp electric potential differences measured at electrodes with respect to a single common reference electrode (e.g., the reference can be linked earlobes, the toe, or one of the elec trodes included in ); the matrix is the lead field matrix corresponding to voxels; is the current 1 E N × ∈ \ Φ E N ) Φ ) 1 ( 3 EV NN × ∈ K \ V N ( 3 V N × ∈ J \ Page 1 of 16 Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 2 of 16 density; c i s a s c a l a r a c c o u n t i n g f o r t h e p h y s i c s n a t u r e o f e l e c t r i c p o t e n t i a l s w h i c h a r e determined up to an arbitrary constant; and 1 denotes a vector of ones, in this case . Typically , and . 1 E N × ∈ 1 \ N E NN  ... = K V 19 E N ≥ 12 1 22 2 12 ... ... ... V EE TT T N TT T N TT NN k k k In Eq. 1 , the structure of K is: Eq. 2: 11 21 V E V T N N ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ kk kk kk where the superscript “ T ” denotes transposition; and , for and for 31 ij × ∈ k \ 1... E i = 1... V j N = , corresponds to the scalp potentials at the i -th electrode due to three ortho gonal unit strength dipoles at voxel j , each one oriented along the coordinate axes x , y , and z . For instance, in infinite homogeneous medium with conductivity : σ Eq. 3: () 3 1 Ei Vj Ei Vj πσ − − rr rr 4 ij = k 31 , Vj × ∈ rr \ where are position vectors for the i -th scalp electrode and for the j -th voxel, respectively. As another example, for the case of a homogeneous conducting sphere in air, the lead field is: Ei Eq. 4: () () () 3 1 2 4 ij Ei Vj E i Ei Vj Ei Vj Ei T Ei Ei Vj Ei Ei Vj Ei Ei V j Ei V j πσ ⎡ ⎤ −− + − ⎢ ⎥ =+ ⎢ ⎥ ⎡ ⎤ −− + − − ⎣ ⎦ ⎣ ⎦ rr r rr rr r rr r rr r r r r rr k In the previous equations, the following notation was used: Eq. 5: () ( 2 TT tr tr == XX X X X ) where tr denotes the trace, and X is any matrix or vector. If X is a vector, then this is the squared Euclidean no rm; if X is a matrix, then this is the squared Frobenius norm. 2 L Note that K can also be conveniently written as: Eq. 6: () 123 , , , ... , V N = KK K K K where , for 3 E N j × ∈ K \ 1... V j N = , is defined as: Eq. 7: 1 2 ... ... E T j T j j T Nj ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎝⎠ k k K k Ideally, the lead field should correspond to the real head (with re alistic geometry and conductivities). For the EEG prob lem, the voxe ls should correspond to cortical grey matter. For other situations (e.g. EKG), appropria te vo lume conductor models and solution spaces should be used. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 3 of 16 In Eq. 1 , J is structured as: Eq. 8: 1 2 ... ... ⎛⎞ ⎜⎟ ⎜⎟ ⎜⎟ = ⎜⎟ ⎜⎟ ⎜⎟ V N ⎝⎠ j j J j where denotes the current density at the i -th voxel. 31 i × ∈ j \ 3. The reference electrode problem As a first step, before even stating the inverse problem, the reference electrode problem will be solved , by estimating “ c ” in Eq. 1 . Given and , the reference electrode problem is: Φ KJ Eq. 9: 2 min c c −− KJ 1 Φ The solution is: Eq. 10: () T T c =− 1 KJ 11 Φ Plugging Eq. 10 into Eq. 1 gives: Eq. 11: = HH K Φ J where: Eq. 12: T T =− 11 HI 11 is the average reference operator, also known as the centering matrix, and is the identity matrix. EE NN × ∈ I \ This result establishes the fact that any in verse solution (of any form, not necessarily linear) will not depend on the reference electrode. Henceforth, it will be assumed that the EEG measurements and the lead field are average reference transformed, i.e.: Eq. 13: ← ⎧⎫ ⎨⎬ ← ⎩⎭ H KH K ΦΦ and Eq. 1 will be rewritten as: Eq. 14: = KJ Φ Note that H plays the role of the identity matr ix for EEG data. It actually is the identity matrix, except for a null eigenvalue corresponding to an ei genvector of ones (see Eq. 12 ), accounting for the reference electrode constant. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 4 of 16 4. A family of discrete, 3D distribu ted linear imaging methods with exact, zero error localization The family of linear imaging methods co nsidered here is parameterized by a symmetric matrix , such that: EE NN × ∈ C \ Eq. 15: () 12 ˆ TT ii i i Φ − ⎡⎤ = ⎢⎥ ⎣⎦ jK C K K C where is any estimator for the electric neuronal activity at the i -th voxel, not necessarily current density (e.g. it can be standardized current de nsity, as in Pascual-Marqui 2002). 31 ˆ i × ∈ j \ Note that in the case of MEG, C must be non-singular . In the case of EEG, C must be of rank ( , with its null eigenvector equal to a vector of ones (accounting for the reference constant). ) ) 1 E N − Note that in Eq. 15 , the symmetric matrix ( is of dimension , and the notation T ii KC K 33 × () 12 T ii − ⎡ ⎢ ⎣⎦ KC K ⎤ ⎥ indicates the symmetric square root inverse. In the particular case of MEG in a spherical head model, the matrix is of rank two, and its symmetric square root pseudo-inverse must be used. () ii K T KC Localization inference in neuroimaging is typically based on the search for large values of the power (squared amplitude) of the estimator for e lectric neuronal activity, i.e. 2 ˆ i j . In order to test the localization properties of a linear imaging method, consider the case when the actual source is an arbitrary point-test source at th e j -th voxel. This mea ns that: Eq. 16: j Φ = KA where is defined in j K Eq. 7 above, and is an arbitrary non-zero vector (containing the dipole moments). 31 × ∈ A \ Plugging Eq. 16 into Eq. 15 and taking the squared amplitude gives: Eq. 17: () 2 ˆ TT T T ij i i i i + = j j AKC K K C K K C K A where the superscript “+” denotes the Moore-Penrose pseudoinverse (which is equal to the common inverse if the matrix is non-singular). Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Following the same type of derivations as in Greenblatt et al ( 2005), the derivative of 2 ˆ i j in Eq. 17 with respect to is: i K Eq. 18: () () () 2 ˆ 2 2 i TT T jj i i i i TT T T T ii i i j j ii i + ++ ⎧⎫ ∂ ⎪⎪ = ⎪⎪ ∂ ⎨⎬ ⎪⎪ − ⎪⎪ ⎩⎭ j CK AA K C K K CK K C K KC K KC K A AKC K KC K It can be easily shown that th is derivative is zero when is equal to , demonstrating that this family of methods prod uces exactly localized maxima to point-test sources anywhere in the brain, i.e. this family of linear imaging methods attains exact, zero error localization. i K j K Note that the choice: Eq. 19: () T α + =+ CK K H gives the sLORETA method (Pascual-Marqui 2002), where is the regularization parameter. 0 α ≥ Note that these results can be applied in a straightforward manner to the case where the current density orientation is known (i.e. known cortical geometry), but with unknown current density amplitude. 5. Unbiase d localization for sLORETA As in the previous section, consider the ca se when the actual source is any arbitrary point-test source at the j -th voxel, but now the measurements are contaminated with measurement and biological noise. This means that: Eq. 20: j εε Φ Φ =+ + J KA K where represents the measurement noise and the biological noise. It will be assumed that both noise sources are zero mean an d independent, with covariance matrices: ε Φ ε J Eq. 21: () cov εσ ΦΦ = H Eq. 22: () cov εσ = JJ I This gives the following expected covariance matrix for the measurements: Eq. 23: () cov TT T jj σσ ΦΦ ΦΣ == + + J KA A K H K K The corresponding expected square amplitude then is: Eq. 24: () () () () () 2 12 12 ˆ TT T T T T ii i i j j i i i TT T T T jj i i i i Et r tr σσ σσ Φ Φ − − + ⎧⎫ ⎛⎞ ⎡⎤ =+ + ⎜⎟ ⎪⎪ ⎢⎥ ⎪⎪ ⎣⎦ ⎝⎠ ⎨⎬ ⎡⎤ ⎪⎪ =+ + ⎢⎥ ⎪⎪ ⎣⎦ ⎩⎭ J J j K CK K C K A A K H KK CK K CK KA A K H K K C K K C K K C Page 5 of 16 Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 6 of 16 The derivative of 2 ˆ i E ⎛⎞ ⎜⎟ ⎝⎠ j in Eq. 24 with respect to is: i K Eq. 25: () () () () () 2 ˆ 2 2 i TT T T jj i i i i TT T T T T ii i i j j ii i E σσ σσ Φ Φ + + + ⎧⎫ ⎛⎞ ∂ ⎜⎟ ⎪⎪ ⎝⎠ ⎪⎪ =+ + ⎨⎬ ∂ ⎪⎪ ⎪⎪ −+ + ⎩⎭ J J j CK A AK H K K C K KC K K C K KC K KCK A A K H K K C K KC K It can be easily shown that the derivative in Eq. 25 is zero for the sLORETA case, when the parameter matrix is: Eq. 26: T σ σ Φ + ⎛⎞ =+ ⎜⎟ ⎜⎟ ⎝⎠ J CK K H and when is equal to , thus demonstrating that sLOR ETA produces exactly localized maxima to point-test sources anywhere in the brain, even in the presenc e of noise, i.e. sLORETA is unbiased. i K j K T h i s n e w r e s u l t i s t o b e c o n t r a s t e d w i t h t h o s e p u b l i s h e d b y S e k i h a r a e t a l ( 2 0 0 5 ) a n d Greenblatt et al (2005). They showed that under pure measurement noise, sLORETA is biased, and only attains exact localization under ideal conditions of no noise. They did not consider the more realistic case where the brain in general is always active, as modeled here by the biological noise. Under these arguably much more realistic conditions, sLORETA is unbiased. 6. eLORETA: exact low resolution br ain electr omagnetic tomography The eLORETA method was developed and offici ally recorded as a working project at the University of Zurich in Ma rch 2005. A description (includi ng the official registration date) can be obtained from the University of Zurich server at: http://www.research-projects.unizh.ch/p6990.htm An additional reference to eLORETA is: Roberto D. Pascual-Marqui, Alberto D. Pascua l-Montano, Dietrich Lehmann, Kieko Kochi, Michaela Esslen, Lutz Jancke, Peter Anderer, Be rnd Saletu, Hideaki Tanaka, Koichi Hirata, E. Roy John, Leslie Prichep. Exact low res olution brain electromagnetic tomography (eLORETA). Neuroimage 2006, Vol 31, Suppl. 1, page:S86 Consider the general weight ed minimum norm solution (see, e.g. Pascual-Marqui 1999): Eq. 27: ˆ Φ = J T with: Eq. 28: () 11 TT α + −− =+ TW KK W K H where denotes the symmetric weight ma trix, and denotes the regularization parameter. () ( 33 VV NN × ∈ W \ ) 0 α ≥ Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” T h e p a r t i c u l a r c a s e o f i n t e r e s t h e r e w i l l only consider a structured block-diagonal weight matrix W , where all matrix elements are zero except for the diagonal sub-blocks denoted as , the i -th voxel, with . 33 i × ∈ W \ 1... V iN = Note that for , this is a genuine solution, in the sense that 0 α = ˆ J is a direct estimator for the current density, and it re produces exactly the measurements. In other words, for : 0 α = Eq. 29: ˆ ΦΦ ⎧⎫ == ⎪⎪ ⎨⎬ = ⎪⎪ ⎩⎭ KJ KT HK T The current density estimator at the i -th voxel then is: Eq. 30: () 11 ˆ TT ii i α Φ + −− =+ jW K K W K H Based on the results of the previous section (entitled “ A family of discrete, 3D distributed linear imaging methods wi th exact, zero error localization ”), by comparing Eq. 30 with Eq. 15 , exact, zero error localization is attained with weights satisfying: Eq. 31: () 12 1 TT ii i α + − ⎡⎤ =+ ⎢⎥ ⎣⎦ WK K W K H K This result is easily derived by noting that Eq. 30 matches Eq. 15 when: Eq. 32: () 1 T α + − ←+ CK W K H and: Eq. 33: () 12 1 T ii i − − ← WK C K The weights satisfying the system of equations given by Eq. 31 define the eLORETA method, which is a genuine solution to the in verse problem (not merely a linear imaging method), and attains exact, zero error locali zation. Additionally, eLORETA is standardized by definition, meaning that its theore tical expected variance is unity. Furthermore, following the derivations as in the previous section entitled “ Unbiased localization for sLORETA ”, it can easily be shown that eL ORETA is unbiased in the presence of measurement and structured biological noise of the form: Eq. 34: () 1 cov εσ − = JJ W Unfortunately, such a structure on background brain activity (the so-called biological noise) is determined by the physics properties of the head model and the laws of electrodynamics, and migh t have little relation to electrophysiological reality. This might be seen as a disadvantage of eLORETA as compared to sLORETA. Page 7 of 16 Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 8 of 16 7. An alternative theoretical a ppr oach to eLORETA, including numerical methods 7.1. The classical weighted minimum norm tomography Consider the regularized, weig hted minimum norm problem: Eq. 35: 2 min T α ⎡⎤ −+ ⎣⎦ J KJ J WJ Φ ) where denotes a given symmetr ic weight matrix, and denotes the regularization parameter. () ( 33 VV NN × ∈ W \ 0 α ≥ The solution is linear: Eq. 36: ˆ = J T Φ with: Eq. 37: () 11 TT α + −− =+ TW KK W K H where the superscript “+” denotes the Moore-Penrose pseudoinverse (which is equal to the common inverse if the matrix is non-singular). The choice gives the classical minimum norm solution. This was t he first 2D distributed linear solution introd uced in ME G by Hämäläinen and Ilmoniemi (1984). Some of the images in that publication show that when the solution space is parallel to th e measurement space, point-test sources are corre ctly localized, albeit with low resolution. = WI However, when the solution space is extend ed to 3D, the minimum norm solution is utterly incapable of correct localization of depth. This was clarified in Pascual-Marqui (1999), where it was shown that the minimum norm solution is harmonic, and harmonic functions attain their extreme va lues on the boundary of thei r domain of definition. This means that deep sources are always incorre ctly localized to the outermost cortex. Another popular choice is depth weighting for the 3D solution space, i.e. larger weights are assigned to deeper sources, with the hope of correcting depth localization erro r. These solutions achieve lower localization erro r than the classical minimum norm, but their errors are still significant, no matter what inverse power for dept h weighting is used. The weighted minimum norm method that uses combined depth weighting and Laplacian smoothing, known as LORETA (l ow resolution brain electromagnetic tomography; Pascual-Marqui et al 1994), achiev ed the lowest localization error up to the present, among linear solutions. Yet, the m eth od has non-zero error, but quite lower than the two previous methods. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 9 of 16 7.2. eLOR ETA: optimal weights that produ ce exact localization The regularized problem in Eq. 35 was presented from a “functional analysis” point of view. Alternatively, a Bayesian point of vi ew renders the same formulation, where the quadratic functional in Eq. 35 is part of the posterior density, with: Eq. 38: noise α = H Φ Σ being the covariance matrix for the noise in the measurements, and: Eq. 39: 1 − = J W Σ being the “ a priori ” covariance matrix for the current density J . Based on the linear relation in Eq. 14 , extending it to include possible additive noise in the measurements, making use of Eq. 38 and Eq. 39 , and assuming independence of neuronal activity and measurement noise, the co variance matrix for the electric potential is: Eq. 40: 1 Tn o i s e T α − =+ = + J KK K W K H ΦΦ ΣΣ Σ Based on the linear relation in Eq. 36 , and making use of Eq. 40 , the covariance matrix for the estimated current density is: Eq. 41: () () () () () ˆ 1 11 1 1 11 1 T TT TT T T TT α αα α α − ++ −− − − + −− − ⎧⎫ = ⎪⎪ ⎪⎪ =+ ⎪⎪ ⎨⎬ =+ + + ⎪⎪ ⎪⎪ =+ ⎪⎪ ⎩⎭ J TT TK W K HT W K KW K H KW K H KW K H KW W K KW K H KW Φ ΣΣ 1 − When W is restricted to be a block-diagonal matrix, with the j -th block denoted as , for 33 j × ∈ W \ 1... V j N = , then the solution to the problem: Eq. 42: () 2 2 11 ˆ min min TT α + −− −= − + J WW I I W K KW K H KW Σ 1 − produces an inverse solution ( Eq. 36 and Eq. 37 ) with zero localization error. Zero localization error is defined in this st udy as follows: For a gi ven point-test source anywhere in the solution space, with a rbitrary orientation, compute th e extracranial EEG/MEG measurements, give them to the line ar inverse solution, threshold the inverse solution to the absolute maximum of the amplit ude of the current density vector field, and compute as localization error the distance be tween the actual point-test source and the position of the absolute maximum. This property has not been achieved by any previously published discrete 3D distributed linear solution. Note that the covariance matrix for the estimated current density ( Eq. 41 ) is not the resolution matrix of Backus and Gilbert. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 10 of 16 The solution to the problem in Eq. 42 satisfies the following set of matrix equations: Eq. 43: () 21 , 1... TT jj j V f or j N α + − =+ = WK K W K H K where the matrix is defined in j K Eq. 7 . The following simple iterative algorithm (i n pseudo-code) converges to the block- diagonal weights W that solve the problem in Eq. 42 and equivalently satisfies Eq. 43 : 1. Given the average reference lead field K and a regularization parameter , initialize 0 α ≥ the block-diagonal weight matrix W as the identity matrix. 2. Set: Eq. 44: () 1 T α + − =+ MK W K H 3. For 1... V j N = do: Eq. 45: SymmSqrt T jj j ⎡⎤ = ⎣⎦ WK M K Comment: denotes the symmetric square root of the matrix . SymmSqrt T jj ⎡⎤ ⎣⎦ KM K T jj ⎡⎤ ⎣⎦ KM K 4. Go to step 2 until convergence (negligible changes in W ). Finally, the block-diagonal matrix W produced by this algorithm should be plugged into the pseudoinverse matrix T (in Eq. 37 ) . T h i s i s d e n o t e d a s t h e e LORETA inverse solution. 7.3. eLORETA for EEG with known current densi ty vector orientation, unknown amplitude The average reference forward EEG equation ( Eq. 14 ) is now written as: Eq. 46: Φ == KJ KNL with: Eq. 47: = J NL where contains the current density amplitudes at each voxel, and contains the outward normal vectors to the cort ical surface at each voxel. Note that the columns of N , denoted as for 1 V N × ∈ L \ () 3 VV NN × ∈ N \ () 31 V N j × ∈ N \ 1... V j N = are: Eq. 48: 1 2 3 12 3 ; ; ; ... ; ... ... ... ... V V N ⎛⎞ ⎛⎞ ⎛⎞ ⎛ ⎞ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ == = = ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜ ⎟ ⎜⎟ ⎜⎟ ⎝⎠ ⎝ ⎠ ⎝⎠ ⎝⎠ N 0 0 n0 0 0 0n 0 n 00 NN N N 0 00 0 n 00 0 where is a vector of zeros, and is the normal vector at the j -th voxel, i.e.: 31 × ∈ 0 \ 31 j × ∈ n \ Eq. 49: 1 T jj = nn In this section, N is assumed to be known. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 11 of 16 The regularized, weighted minimum norm problem is: Eq. 50: 2 min T α Φ ⎡⎤ −+ ⎣⎦ L KNL L WL + 1 − where in this case denotes a given symmetric weight matrix, and denotes the regularization parameter. VV NN × ∈ W \ 0 α > The solution is linear: Eq. 51: ˆ Φ = LT with: Eq. 52: () () () () 11 TT α + −− =+ T W KN KN W KN H Following similar lines of reasoning as in th e previous section, the covariance matrix for the electric potential is: Eq. 53: () () () () 1 TT noise α ΦΦ ΣΣ Σ − =+ = L KN KN KN W KN H where: Eq. 54: 1 Σ − = L W is the “ a priori ” covariance matrix for the current density amplitudes L . In addition, the covariance matrix for the estimated current density is: Eq. 55: () () () () () 11 ˆ TT α Σ + −− =+ L W K N K NW K N H K NW When W is restricted to be a diagonal matrix, with the j -th element denoted as , for j W 1... V j N = , then the solution to the problem: Eq. 56: () () () () () 2 2 11 ˆ min min TT α Σ + −− −= − + L WW I I W K N K NW K N H K NW 1 − produces an inverse solution ( Eq. 51 and Eq. 52 ) with zero localization error. The solution to the problem in Eq. 56 satisfies the following set of equations: Eq. 57: () () () 1 , 1... T T j V jj f or j N α + − =+ WK N K W K H K N = where the vector corresponds to the j -th column of . () 1 E N j × ∈ KN \ () KN The following simple iterativ e algorithm (in pseudo-code) converges to the diagonal weights W that solve the problem in Eq. 56 and equivalently satisfies Eq. 57 : 1. Given the average reference lead field K , the cortical normal vectors N , and a regularization parameter , initialize the diagonal weight matrix W as the identity 0 α ≥ matrix. 2. Set: Eq. 58: () () () 1 T α + − =+ MK N W K N H 3. For 1... V j N = do: Eq. 59: () () T j jj = WK N M K N 4. Go to step 2 until convergence (negligible changes in W ). Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 12 of 16 Finally, the diagonal matrix W produc ed by this algorith m should be plugged into the pseudoinverse matrix T (in Eq. 52 ). This is denoted as the e LORETA inverse solution. 7.4. eLORETA for MEG with fully unknown current density vector field This case follows the same derivations as given above for the case “EEG with fully unknown current density vector field”. The forward MEG equation has similar form to Eq. 14 . For the MEG case, would represent the magnetometer or gradiometer measurements, would represent the magnetic lead field, and is exactly the same current density vector field (common to both EEG and MEG). Φ K J In the MEG case, there is no reference electrode constant to be a ccounted for. The consequence is that the EEG regularization term appearing in some of the equations above ( ( α H ) i 3 Eq. 37 to Eq. 44 ) must be changed to the MEG regularization term ( , where I is the identity matrix. ) α I In the case of spherical head models, care must be taken in the MEG case because only the tangential part of the current density vector field is non-sile nt. The same occurs in realistic head models, in areas that are quasi-sphe rical. This implies that all calculations at the voxel level have only rank=2 for MEG. Therefore, inverse and symmetric square-root matrix computations should be made via the singular value decomposition (SV D), ignoring the smallest eigenvalue if it is numerically negligible relative to the largest eigenvalue. In particular, consider the algorithm involving Eq. 44 and Eq. 45 . Note that Eq. 44 makes use of the inverse of the weight matrix, which consists of the inverses of all block-diagonal submatrices. In the quasi- spherical MEG case, these submatrices have rank=2. Referring to 33 × Eq. 45 , consider the SVD of the matrix of interest: Eq. 60: 3 1 TT jj i i i λ ΓΓ = ⎡⎤ = ⎣⎦ ∑ KM K where are the orthonormal eigenvectors, and are the eigenvalues. Then 31 i Γ × ∈ \ 12 λλ λ ≥≥ Eq. 45 should be replaced by: Eq. 61: () 2 31 1 3 1 , , T ii i Symm Sqrt i T jj j T ii i i if oth erwi se λλ λ λ ΓΓ ΓΓ = = ⎧ < ⎪ ⎪ ⎡⎤ == ⎨ ⎣⎦ ⎪ ⎪ ⎩ ∑ ∑ WK M K ε ε where depends on the numerical precision of the calculations (typically ). ε 5 10 − ≤ Moreover, the inverse of (see j W Eq. 61), which is needed in Eq. 44 , and later on after convergence in Eq. 37 for the final inverse solution, should be calculated as the Moore- Penrose pseudoinverse (ignoring the smallest ei genvalue if it is numerically negligible relative to the largest eigenvalue). Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 13 of 16 Given these provisions and modifications, th e discrete 3D distributed linear solution known as e LORETA is given by Eq. 36 and Eq. 37 with the weights defined by the solution the problem in Eq. 42 , obtained with the algorithm specified by Eq. 44 and Eq. 45 . 7.5. eLORETA for MEG with known curr en t density vector orientation, unknown amplitude This case follows the same derivations as given above for the case “EEG with known current density vector orient ation, unknown amplitude”. The forward MEG equation in this case has similar form to Eq. 46. For the MEG case, would represent the magnetometer or gradiometer measurements, K would represent the magnetic lead field, N would contain the outward normal vectors to the cortical surface at each voxel (assumed known), and L contains exactly the same current density amplitudes (common to both EEG and MEG). Φ As explained previously, the EEG regularization term ( appearing in some of the equations above ( ) α H Eq. 52 to Eq. 58 ) must be changed to the MEG regularization term , where I is the identity matrix. () α I The existence of silent MEG sources might occur in practice, especially for quasi- radial sources in quasi-spherical head geometry. Care should be taken to exclude these possible silent sources from the solution space, even if this implies that there are missing cortical patches for the MEG solution space. Given these provisions and modifications, th e discrete 3D distributed linear solution known as e LORETA is given by Eq. 51 and Eq. 52 with the weights defined by the solution to the problem in Eq. 56 , obtained with the algorithm specified by Eq. 58 and Eq. 59 . The general family of linear imaging meth ods is further extended to include data- dependent (adaptive) quasi-linear imaging me thods, also with the exact, zero error localization property. 8. A family of discrete, 3D distributed quasi-linear imaging metho ds with exact, zero error localiza tion: data-dependent ( adaptive) methods Formally, this family of methods is identical to the one defined by Eq. 15 , except that now, instead of defining the parameter matrix C as a g i ve n , fi x ed m at r i x , i t c an , f o r e x a mp le , be taken as the inverse covariance matrix for the measurements, i.e: Eq. 62: () () 1 1 K N T kk k K N ΦΦ ΦΦ + = ⎡⎤ =− − ⎢⎥ ⎣⎦ ∑ C where the subscript “ k ” may index time or any other form of repeated measurements for the data. Note that in the case of MEG, C must be non-singular. In the case of EEG, C must be of rank , with its null eigenvector equal to a vector of ones (accounting for th e ( 1 E N − ) Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 14 of 16 reference constant). If the data so happens to be insufficient, i.e. , or the data happens to be almost deterministic, resulting in a low rank matrix C , then the method will not have exact, zero error localization. K NN ≤ E Eq. 62 corresponds to a single example illustra ting the “adaptive” character of this family of methods. Any data dependent matrix C can be used, such as, for example, the squared inverse covariance matrix for the measurements. Rigorously speaking, this method is not linear because the transformation depends on the data on which imaging is being carried out. 9. Conclusions In Pascual-Marqui (1995, 1999, and 2002), the following arguments were used for selecting the best discrete, 3D distributed, linear tomography: 1. The aim of functional imagin g is localization. Therefore, the best tomography is the one with minimum localization error. 2. In a linear tomography, the localization prop erties can be determined by using point-test sources, based on the principles of linearity and superposition. 3 . If a l i n ea r t om o g ra p h y is in c a pa bl e o f z er o error localization to point-test so urces that are active one at a time, then the tomography will certainly be incapable of zero er ror localization to two or more simultaneously active sources. Here we present a general family of linear imaging methods with exact, zero error localization to point-test sources. We show that one particular member of this family, sLORETA (Pascual-Marqui 2002) has no localization bias in the presence of measurement noise and biological noise. We introduce a new particular member of this family, denoted eLORETA. This is a genuine inverse solution and not merely a li near imaging method. We show that it has exact, zero error localization in the presence of measurement and structured biological noise. We derive and construct the method using two different approaches, and give practical algorithms for its estimation. We present a general family of quasi-linear imaging methods that are data-dependent (adaptive). We also show that they are endo wed with the exact, zero error localization property. These results are expected to be of value to the EEG/MEG neuroimaging community. Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” 10. References Baillet S, Mosher JC, Leahy RM: Elec tromagnetic Brain Ma pping. IEEE Signal Processing Magazine 18:14-30, 2001. Dale AM, Liu AK, Fischl BR, Buckner RL , Belliveau JW, Lewine JD, Halgren E: Dynamic statistical parametric mapping: co mbining fMRI and MEG for high-resolution imaging of cortical activi ty. Neuron 26: 55-67, 2000. Greenblatt, R.E. Ossadtchi, A. Pflieger, M.E. Local Linear Estimators for th e Bioelectromagnetic Inverse Problem. IEEE transa ctions on signa l processing 2005, 53: 3403- 3412. Haalman I, Vaadia E: Dynamics of neuronal interactions: relation to behavior, firing rates, and distance between neurons. Human Brain Mapping 5: 249-253, 1997. Hämäläinen, M.S., and Ilmoniemi, R.J. Interp reting measured magnetic fields of the brain: estimates of current di stributions. Tech. Rep. TKK-F- A559, Helsinki University of Technology, Espoo, 1984. Hämäläinen MS, Hari R, Ilmoniem i RJ, Knuutila J, Lounasmaa OV: Magnetoencephalography - theory, instrumentatio n, and applications to noninvasive studies of the working human brain. Rev. Mod. Phys. 65: 413–497, 1993. Llinas RR: The intrinsic electrophysiolog ical properties of mammalian neurons: insights into central nervous system function. Science 242: 1654-1664, 1988. Martin JH: The collective electrical behavior of cortical neurons: The electroencephalogram and the mechanisms of epilepsy. In Kandel ER, Schwartz JH, Jessell TM (Eds.) Principles of Neural Science, Pr entice Hall International, London, pp 777-791, 1991. Mitzdorf U. Current source-density method and application in cat cerebral cortex: investigation of evoked potentials and EEG phenomena. Physiol Rev 1985;65:37-100. Pascual-Marqui RD, Mishel CM, Lehmann D: Low resolution electromagnetic tomography: a new method for localizing electr ical activity in the brain. International Journal of Psychophysiology. 1994, 18: 49-65. Pascual-Marqui RD. Review of methods for solving the EEG inverse problem. International Journal of Bioelectromagnetism 1999; 1:75–86. ( http://www.unizh.ch/keyinst/NewLORETA/T echnicalDetail s/TechnicalDetails.pdf ) and ( http ://www.ijbem.org/ ). Pascual-Marqui RD: Reply to comments by Hämäläinen, Ilmoniemi and Nunez. In Source Localization: Continuing Discussion of the Inverse Pr oblem (W. Skrandies, Ed.), pp. 16-28, ISBET Newsletter No.6 (ISSN 0947-5133), 1995. ( http://www.unizh.ch/keyinst/NewLORETA /BriefHistory/LORETA-NewsLett2b.pdf ). Page 15 of 16 Cite as: “R.D. Pascual-Marqui: D iscrete, 3D distributed, line ar imaging me thods of elec tric neuronal activity. Part 1: exact, z ero error localization. arXiv:0710.3341 [math-ph], 2007-October-17, ” Page 16 of 16 Pascual-Marqui, R.D., 2002. Standardized low-resolution brain electromagnetic tomography (sLORETA): technical details. Methods Find. Exp. Clin. Pharmacol. 24 (Suppl D.), 5 –12. ( http://www.unizh.ch/keyinst/NewLORETA/sLORETA/sLORETA-Math01.pdf ). Roberto D. Pascual-Marqui, Alberto D. Pa scual-Montano, Dietrich Lehmann, Kieko Kochi, Michaela Esslen, Lutz Jancke, Peter An derer, Bernd Saletu, Hideaki Tanaka, Koichi Hirata, E. Roy John, Leslie Prichep. Exact lo w resolution brain elec tromagnetic tomography (eLORETA). Neuroimage 2006, Vol 31, Suppl. 1, page:S86 J. Sarvas, “Basic mathematical and elec tromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11-22, 1987. K. Sekihara, M. Sahani, and S. S. Nagarajan, “Localization bias and spatial reso lution of adaptive and nonadaptive sp atial filters for MEG source reconstruction,” Neuroim ag. 2005, 25:1056-1067. Sukov W, Barth DS: Three-dimensional analysis of spontaneous and thalamically evoked gamma oscillations in auditory co rtex. J. Neurophysiol. 79: 2875–2884, 1998.