Europe PMC

This website requires cookies, and the limited processing of your personal data in order to function. By using the site you are agreeing to this as outlined in our privacy notice and cookie policy.

Abstract 


Two high resolution methods solving inverse problems potentially ill-posed, named 4-MUSIC and 4-RapMUSIC, are proposed. They allow for localization of brain current sources with unconstrained orientations from surface electro-or magneto-encephalographic data using spherical or realistic head geometries. The 4-MUSIC and 4-RapMUSIC methods are based on i) the separability of the data transfer matrix as a function of location and orientation parameters and ii) the fourth order (FO) virtual array theory. In addition, 4-RapMUSIC uses the deflation concept extended to FO statistics accounting for the presence of potentially but not totally coherent sources. Computer results display the superiority of the 4-RapMUSIC approach in different situations (two closed sources, additive Gaussian noise with unknown spatial covariance, ...) especially over classical algorithms.

Free full text 


Logo of halLink to Publisher's site
Conf Proc IEEE Eng Med Biol Soc. Author manuscript; available in PMC 2007 Dec 5.
Published in final edited form as:
PMCID: PMC2099315
HALMS: HALMS183642
PMID: 17946634

Fourth order approaches for localization of brain current sources

Abstract

Two high resolution methods solving inverse problems potentially ill-posed, named 4-MUSIC and 4-RapMUSIC, are proposed. They allow for localization of brain current sources with unconstrained orientations from surface electro- or magneto-encephalographic data using spherical or realistic head geometries. The 4-MUSIC and 4-RapMUSIC methods are based on i) the separability of the data transfer matrix as a function of location and orientation parameters and ii) the Fourth Order (FO) virtual array theory. In addition, 4-RapMUSIC uses the deflation concept extended to FO statistics accounting for the presence of potentially but not totally coherent sources. Computer results display the superiority of the 4-RapMUSIC approach in different situations (two closed sources, additive Gaussian noise with unknown spatial covariance, …) especially over classical algorithms.

Keywords: Source localization, Brain activity, Statistical Signal Processing, EEG, Epilepsy, Neurology

I. INTRODUCTION

Electroencephalography (EEG) and magnetoencephalography (EEG) are two complementary techniques respectively measuring electrical potentials and magnetic fields produced by neuronal activity, at the surface of the head. Localization of neuronal activity sources requires to solve the inverse problem which is underdetermined in theory, as the number of sources is generally larger than the number of sensors, and so ill-posed. Conversely, when the number of sources is assumed to be lower than the number of scalp measurements, the problem is overdetermined and has a unique solution. In order to solve the EEG/MEG inverse problem both a model of neuronal sources and a model of the head are required. The current dipole is the most commonly used model for a source of electrical activity in the brain. Head models aim at representing geometrical and electrical properties of the different tissues composing the volume conductor. Various models were proposed going from concentric homogeneous spheres with isotropic conductivities to realistically shaped models with refined tissue conductivity values.

Numerous array processing methods were developed to estimate multidimensional parameters of sources, allowing among other things the localization of brain current dipoles from scalp measurements. Among subspace approaches, the Second Order (SO) MUSIC (MUltiple SIgnal Classification) method [10] was proposed for overdetermined mixtures of sources. This approach gave rise to several variants aimed at improving performances. On the one hand, among time MUSIC-like methods, one can mention the extension of the original MUSIC algorithm to Fourth Order (FO) statistics (so-called MUSIC4) able to deal with the case of under-determined mixtures of sources [9]. One can also mention sequential approaches such [7], which are based on both SO statistics and deflation concept to increase localization resolution. In particular, the RapMUSIC algorithm [7] takes advantage of the factored matrix formulation of the transfer relationship between deep sources and scalp data by separating nonlinear (location) from quasilinear (orientation) source parameters in order to reduce computing time [3]. On the other hand, Time-Frequency (TF) approaches were proposed (see [1] for instance) to improve the resolution of the localization in the case of very closed sources with spectral nonstationary properties. Finally, besides subspace methods, other localization methods applied to EEG and MEG data were reported (review in [6]).

Three remarks can be made from this brief overview. First, most of the aforementioned array processing methods are based on SO statistics, which implicitly imply that sources are Gaussian. In practice, this assumption is not justified physiologically and information available at higher orders can be utilized. Second, TF approaches showed to be not applicable when sources have quasi-identical TF supports. Third, time SO techniques cannot deal with underdetermined mixtures of sources or with a Gaussian noise of unknown spatial covariance.

Higher Order (HO) methods inherently account for these limitations. However, to date, there is no attempt to propose a FO method that takes advantage of i) the separability of the data transfer matrix as a function of location and orientation parameters and/or ii) the deflation concept. The intent of the present paper is to describe two new FO MUSIC-like methods addressing these issues. These methods are referred to as 4-MUSIC and 4-RapMUSIC, respectively. Both are based on i) the factored matrix formulation of the data transfer function and on ii) the FO virtual array theory whose relevance has already been displayed in radiocommunications contexts [2], and both account for the presence of possibly but not totally coherent sources.

II. ASSUMPTIONS AND NOTATIONS

A. The Problem Formulation

We assume that a (N × K) realization of a N-dimensional random process {x(k)} is observed. Moreover, each random vector x(k) is given by:

x(k)=A(Θ)s(k)+ν(k)
(1)

where {s(k)}is a P-dimensional source vector process, which observations are the time courses of P current dipoles; A(Θ) = [a(θ1), …, a(θp)] is the (N × P) static mixing matrix that depends on Θ = {θ1, …, θP}, the collection of the multi-parameters of the P sources; and {ν(k)}is the noise vector process that is assumed to be Gaussian and statistically independent of the source vector process.

In EEG (or MEG) applications, each column vector a(θ) of the static mixing matrix represents the electrical (or magnetic) field generated at all scalp sensors by a current dipole with a unit time course localized at a given position ρ for a given orientation [var phi]. Vector a(θ) can be written as the product of a (N × 3) gain matrix G(ρ) and the orientation vector [var phi]:

a(θ)=G(ρ)φ
(2)

where the multi-parameter vector θ = [ρT [var phi]T]T of the considered current dipole includes the nonlinear location parameter ρ and the quasilinear orientation parameter [var phi]. Vector a(θ) will be referred to as the source localizing vector in the sequel.

Although both methods we developed can be applied to both EEG and MEG data, and to both spherical and realistic head models, the present work is focused on the EEG context with a spherical head model. An analytic expression for the gain matrix G(ρ) can be found in [8].

B. FO statistics

For the sake of convenience, the present work is limited to stationary and ergodic data. In that case, the (N2 × N2) quadricovariance matrix [9] [2], Qx, of process {x(k)}can be easily estimated from the scalp data. Given the multilinearity property of cumulants [4], Qx has a special algebraic structure, with several matrix redundancies. If statistical independency between sources and noise is assumed, this property can be expressed as follows:

Qx=[A2]Qs[A2]T=AQsAT
(3)

where Qs, A2def=AA and [multiply sign in circle] are the (P2 × P2) quadricovariance matrix of {s(k)}, the Kronecker square of matrix A and the Kronecker product operator, respectively.

III. ALGORITHM

A. The FO null-spectrum

Compute the Eigen Value Decomposition (EVD) of the symmetrical matrix Qx as follows:

Qx=[EsEν][Ls000][EsEν]T
(4)

where Ls, Es, Eν and R denote the ( R×R) real-valued diagonal matrix of the non zero eigenvalues of Qx, the ( N2×R) matrix of the associated eigenvectors, the (N2×(N2R)) matrix of the eigenvectors associated with the zero eigenvalues of Qx and the rank of Qx respectively. Since matrix Qx is symmetrical, each column of Es is orthogonal to each column of Eν. Now Span {A}=Span{Es}, therefore each column of A is orthogonal to each column of Eν. Given θP the location/orientation parameters of the p-th source and a(θp)2def=a(θp)a(θp) the FO virtual source localizing vector, that is, the true localizing vector of the p-th source for the corresponding FO virtual array [2]. The column vectors a(θp)[multiply sign in circle]2 (1 ≤ pP) appear in matrix A, and so they are orthogonal to each column of Eν. Thus the normalized FO null-spectrum (null-polyspectrum) criterion can be defined as follows:

J1(θ)=([a(θ)2]TΠ0[a(θ)2])/a(θ)22
(5)

where Π0=EνTEν is the noise subspace projector of matrix Qx, and the P roots of J1 correspond asymptotically to the P vectors θP. Using properties of the Kronecker product, J1, after inserting (2) into (5), becomes:

J2(θ)=ΦTG(ρ)TΠ0G(ρ)ΦΦTG(ρ)TG(ρ)Φ
(6)

where Φdef=φ2 and G(ρ)def=G(ρ)2 are the FO source orientation vector and the FO gain matrix, respectively. Parameters θP = [ρpT[var phi]pT]T (1 ≤ pP) can then be obtained by i) looking for the P global minima ρp of the function in ρ defined by the minimum eigenvalue of matrix G(ρ)TΠ0G(ρ) in the metric G(ρ)TG(ρ) and by ii) identifying each vector Φ(1 ≤ pP) as the eigenvector associated with the minimum eigenvalue of matrix G(ρp)TΠ0G(ρp) in the metric G(ρp)TG(ρp) [5]. An algorithm is proposed in section III-B in order to compute vector [var phi] from Φ. Consequently the orientation parameters are deduced from the location parameters, which modifies the 6D-optimization of criterion (6) to a 3D-optimization. Therefore, criterion (6) can be concentrated with respect to ρ, leading to:

J3(ρ)=λ{[G(ρ)TG(ρ)]1G(ρ)TΠ0G(ρ)}
(7)

where λ{B} denotes the minimum eigenvalue of matrix B. The computational cost can considerably be reduced again if criterion J3 is replaced by the following criterion:

J4(ρ)=det{G(ρ)TΠ0G(ρ)}det{G(ρ)TG(ρ)}
(8)

where det{B} denotes the determinant of matrix B.

B. The 4-MUSIC and 4-RapMUSIC approaches

In this section, we present the 4-MUSIC and 4-RapMUSIC methods based on i) the separability of the data transfer matrix as a function of location and orientation parameters, and ii) the FO virtual array theory [2], using the low cost function as defined in (8). In addition, the 4-RapMUSIC algorithm exploits the deflation concept that we have extended to FO statistics.

More particularly, the 4-MUSIC algorithm consists in searching for global minimizers of criterion J4. Indeed, if the noise subspace projector was estimated perfectly, such as asymptotically, then the P source location vectors ρp would be directly found as the P global minimizers of (8) over a sufficiently densely sampled grid of the nonlinear parameter space. Then each source orientation vector [var phi]p (1 ≤ pP) can be computed from the previous locations using both following steps. First, let the FO source orientation vector Φp be the normalized quasilinear parameter vector that must multiply G(ρp) on right to produce vector a(θP)[multiply sign in circle]2. It has to be derived from the eigenvector corresponding to the minimum eigenvalue of matrix [G(ρp)TG(ρp)]1G(ρp)TΠ0G(ρp). Secondly, the source orientation vector [var phi]p can be computed from Φp, by i) remodeling it into one (N×N) matrix Bp (the n-th column of Bp is made up from the N consecutive elements of Φp as from the [N(n − 1) + 1]-th one), and ii) diagonalizing it. Indeed, the eigenvector associated with the strongest eigenvalue of Bp, is, up to a sign factor, equal to [var phi]p.

Nevertheless, for a finite number of samples, errors in our statistic estimate reduce (8) to a function with i) a single global minimum that corresponds for instance to the source of maximum Signal-to-Noise Ratio (SNR), and ii) P−1 local minima. Although the global minimum is easily identifiable, it is more difficult to find the P−1 remaining local minima since nonlinear search techniques may miss shallow or adjacent peaks and return to a previous peak. Algorithms have been proposed to solve peak-picking problem (review in [7]), but they rapidly become complex and subjective as the number of sources and the dimensionality of vectors ρp increase [7]. So, a computation strategy such as the following FO deflation concept has to be established to avoid this peak-picking problem, giving rise to the 4-RapMUSIC method.

This latter consists in localizing recursively the P sources. Indeed, the p-th step of the 4-RapMUSIC method allows for identifying the ξ(p)-th source location and orientation vectors. The use of the permutation function ξ of {1,2, …, P} is necessary since the P source localizing vectors a(θp) may be recovered only in the disorder. Indeed, equation (1) shows that the order in which the components of s(k) and the associated columns of A(θ) are set does not change the expression of x(k). The ξ(p)-th source location vector, ρξ(p), can be achieved by searching for the global minimum root of (8) replacing G(ρ) and π0 by Ap12G(ρ) and the noise subspace projector, Πp−1, of matrix [Ap12]Qx[Ap12]T, respectively, where:

j,1jp1,a(θξ(j))=G(ρξ(j))φξ(j)Ap1=[a(θξ(1))a(θξ(p1))]Ap1=INAp1[(Ap1)TAp1]1[Ap1]T
(9)

Note that the rank of the data quadricovariance matrix [Ap12]Qx[Ap12]T is strictly smaller than R. Indeed, the multiplication of the initial quadricovariance Qx on left and right by Ap12 and [Ap12]T, respectively, comes down to remove the contribution of the p − 1 first localized sources from Qx, and consequently to increased the dimension of the noise subspace and so the capacity of source localization. As far as the source orientation vector [var phi]ξ(p) is concerned, it can be computed from the source location vector ρξ(p) as explained in the previous paragraph.

IV. RESULTS

Performances of four MUSIC-like algorithms (2-MUSIC [3], 2-RapMUSIC [7], 4-MUSIC and 4-RapMUSIC) were compared in an EEG context, in three different situations. Simulations were performed using a 3-shell spherical head model (radii were 8 cm for brain, 8.5 cm for skull and 9.2 cm for scalp; brain and scalp conductivities were 33.10−4 S/cm and skull conductivity was 40 times lower) and a set of ten electrodes from 10–20 standard (namely Fz, Cz, Pz, Fp2, F3, F4, P7, P8, T7 and O2). Depending on the situation, P= 1 or 2 independent sources were considered. A physiologically-relevant model, consisting in a network of coupled neuronal populations [12], was used to compute realistic source time courses. The P sources had the same SNR equal to 15 dB, they were arranged in the same z-plane, and their orientations [var phi]p (1 ≤ pP) were randomly fixed such as ||[var phi]p|| = 1. The background noise was considered as temporally and spatially white, except for situation IV-C. Simulation results were averaged over M = 200 realizations.

Two criteria were used to quantify the quality of the source localization for each method. The first one is the Probability of Non Localization (PNL), which is defined by the ratio between the number of realizations for which all the sources are not localized and the total number of realizations M. The second one is the averaged Root Mean Square Error (RMSE), which is defined, for source p, by:

RMSE(θp)=1Mm=1M(min1jP{θpθ^j(m)})
(10)

where M′(M′≤ M) is the number of realizations for which the localization method has succeeded in finding exactly P solutions, and θ^j(m) is the j-th source parameter vector estimated during the m-th realization.

A. Effect of the dipole depth on source localization

The behaviour of the four MUSIC-like algorithms was studied in the presence of a unique source using only six scalp electrodes. Figure 1 displays the variations of the RMSE criterion of the four methods as a function of the source location on the z-axis. It appears that both FO methods localize more precisely the source than both SO methods, whatever source depth. Note that PNL criterion is not reported here because it was quasi-zero for all methods whatever source depth.

An external file that holds a picture, illustration, etc.
Object name is halms183642f1.jpg

Localization of 1 source with 10 electrodes.

B. Localization of poorly spatialy separated sources

The behaviour of the four algorithms was also studied in the case of two poorly spatialy separated sources. Figure 2 presents the RMSE and PNL criteria of the four methods as a function of samples for two sources using ten surface electrodes. Source location parameters were chosen as equal to (ρ1 = [0,0,0.8]T and ρ2 = [0,0,1.12]T (values in centimeters). It clearly appears on figure 2(b) that the 2-MUSIC and 4-MUSIC methods have difficulty in localizing both sources, whereas the 2-RapMUSIC and 4-RapMUSIC algorithms succeed in finding two solutions. Nevertheless, as displayed by figure 2(a), the 2-RapMUSIC method does not localize as precisely both sources as the 4-RapMUSIC algorithm.

An external file that holds a picture, illustration, etc.
Object name is halms183642f2.jpg

Localization of 2 very close sources with 10 electrodes.

C. The case of colored noise

Both FO algorithms were compared to both SO algorithms in the presence of a Gaussian noise with unknown spatial covariance using ten surface electrodes. Two sources were positioned in depth such that their location vectors were given by (ρ1= [0,0,2]T and ρ2 = [0,0,4.4]T respectively. Figure 3 displays the variations of RMSE and PNL criteria for the four methods as a function of the noise spatial covariance factor ρ. Note that the Gaussian noise model employed in this simulation is the sum of an internal noise νin(k) and an external noise νout(k), of covariance matrices Rνin and Rνout respectively such that:

An external file that holds a picture, illustration, etc.
Object name is halms183642f3.jpg

Source localization in the presence of a colored noise.

Rνin(r,q)def=σ2δ(rq)/2Rνout(r,q)def=σ2ρrq/2
(11)

where σ2, ρ, Rν(r,q)def=Rνin(r,q)+Rνout(r,q) are the total noise variance per sensor, the noise spatial covariance factor and the (r, q)-th component of the total noise covariance matrix, respectively.

Figure 3(a) shows that both SO algorithms are sensitive to a Gaussian noise with unknown spatial covariance and are affected as soon as the noise spatial covariance increases beyond 0.1. Indeed, 2-MUSIC and 2-RapMUSIC would theoretically require a perfect knowledge of the noise covariance [10], while 4-MUSIC and 4-RapMUSIC, since they use FO cumulants, are asymptotically insensitive to Gaussian noise, regardless of its space/time color. Indeed, the computer results show that 4-MUSIC and 4-RapMUSIC localize both sources with precision whatever the noise spatial covariance is. Nevertheless, figure 3(b) displays that, for a given number of ten thousand samples, only 4-RapMUSIC succeeds in localizing both sources at each time.

V. CONCLUSION

We proposed in this paper two novel algorithms for brain current source localization, namely the 4-MUSIC and 4-RapMUSIC methods. Computer results showed the superiority of 4-RapMUSIC over 4-MUSIC and classical algorithms such as 2-MUSIC [3] and 2-RapMUSIC [7] for both overdetermined mixtures of sources and a small number of ten electrodes. Forthcoming works will display i) its superiority for larger number of electrodes as used in standard or higher resolution montages, ii) its capacity to localize more sources than surface observations, and iii) its behavior when applied to patients in whom strong hypotheses about localization of epileptic zones are available.

Acknowledgments

This work was supported in part by the Regional Council of Brittany (http://www.region-bretagne.fr) and it is protected by a patent whose reference is no. 05.11668, DV326 (SAIC), November 17 2005.

References

1. BELOUCHRANI A, ABED-MERAIM K, AMIN MG. Time-frequency music. IEEE Signal Processing Letters. 1999 May;6(5):109–110. [Google Scholar]
2. CHEVALIER P, FERREOL A. On the virtual array concept for the fourth-order direction finding problem. IEEE Transactions On Signal Processing. 1999 Sep;47(9):2592–2595. [Google Scholar]
3. FERRARA E, PARKS T. Direction finding with an array of antennas having diverse polarizations. IEEE Transactions On Antennas Propagation. 1983 Mar;31:231–236. [Google Scholar]
4. McCULLAGH P. Tensor Methods in Statistics. Chapman and Hall: Monographs on Statistics and Applied Probability; 1987. [Google Scholar]
5. GANTMACHER FR. The Theory of Matrices. Vol. 1. New York: Chelsea; 1959. [Google Scholar]
6. MICHEL CM, MURRAY MM, LANTZ G, GONZALEZ S, SPINELLI L, GRAVE DE PERALTA R. EEG source imaging. Clinical Neurophysiology. 2004 Oct;115(10):2195–2222. [Abstract] [Google Scholar]
7. MOSHER JC, LEAHY RM. Source localization using Recursively Applied and Projected (RAP) music. IEEE Transactions On Signal Processing. 1999 Feb;47(2):332–340. [Abstract] [Google Scholar]
8. MOSHER JC, LEAHY RM, LEWIS PS. EEG and MEG: Forward solutions for inverse methods. IEEE Transactions On Biomedical Engineering. 1999 Mar;46(3):245–259. [Abstract] [Google Scholar]
9. PORAT B, FRIEDLANDER B. Direction finding algorithms based on high-order statistics. IEEE Transactions On Signal Processing. 1991 Sep;39(9):2016–2024. [Google Scholar]
10. SCHMIDT RO. Multiple emitter location and signal parameter estimation. IEEE Transactions On Antennas Propagation. 1986 Mar;34(3):276–280. [Google Scholar]
11. STOICA P, HANDEL P, NEHORAI A. Improved sequential music. IEEE Transactions On Aerospace and Electronic Systems. 1995 Oct;31(4):1230–1239. [Google Scholar]
12. WENDLING F, BELLANGER J, BARTOLOMEI F, CHAUVEL P. Relevance of nonlinear lumped-parameter models in the analysis of depth-eeg epileptic signals. Biological Cybernetics. 2000;83:367–378. [Abstract] [Google Scholar]

Citations & impact 


Impact metrics

Jump to Citations

Article citations