Keywords: Source localization, distributed sources, EEG, MEG, MUSIC, higher order statistics

simulated EEG signals obtained with physiologically-relevant models of both the sources and the volume conductor show a highly increased performance of our 2q-ExSo-MUSIC method as compared to the classical 2q-MUSIC algorithms.

Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach Gwénaël Birot a,b , Laurent Albera a,b , Fabrice Wendling a,b , Isabelle Merlet a,b, * = Inserm -UMR 642, Campus de Beaulieu -Bat. 22, Rennes, F-35042, France > Université de Rennes 1, LTSI, Campus de Beaulieu -Bat. 22, Rennes, F-35042, France

)>IJH=?J

We propose a new MUSIC-like method, called 2q-ExSo-MUSIC (q ≥ 1).

This method is an extension of the 2q-MUSIC (q ≥ 1) approach for solving the EEG/MEG inverse problem, when spatially-extended neocortical sources (ExSo) are considered. It introduces a novel ExSo-MUSIC principle. The novelty is two-fold: i) the parameterization of the spatial source distribution that leads to an appropriate metric in the context of distributed brain sources and ii) the introduction of an original, ecient and low-cost way of optimizing this metric. In 2q-ExSo-MUSIC, the possible use of higher order statistics (q ≥ 2) oers a better robustness with respect to Gaussian noise of unknown spatial coherence and modeling errors. As a result we reduced the penalizing eects of both the background cerebral activity that can be seen as a Gaussian and spatially correlated noise, and the modeling errors induced by the non-exact resolution of the forward problem. Computer results on 1. Introduction Since the 1980's, source localization methods based on advanced signal processing techniques have been proposed to identify the intracerebral generators underlying surface ElectroEncepghaloGraphic (EEG) or MagnetoEn-cephaloGraphic (MEG) signals. Particularly, in the context of drug resistant partial epilepsy, source localization methods have been widely applied to interictal spikes. Indeed, these paroxysms witness the patient's epilepsy and, compared to ictal discharges, appear in the large majority of cases on surface recordings with a higher Signal-to-Noise Ratio (SNR). For reviews on this topic, the reader may refer to [START_REF] Ebersole | Noninvasive localization of epileptogenic foci by EEG source modeling[END_REF][START_REF] Funke | Magnetoencephalography and magnetic source imaging in epilepsy[END_REF][START_REF] Merlet | Dipole modeling of interictal and ictal EEG and MEG paroxysms[END_REF][START_REF] Michel | EEG source imaging[END_REF][START_REF] Plummer | EEG source localization in focal epilepsy: where are we now?[END_REF][START_REF] Verhellen | EEG source localization of the epileptogenic focus in patients with refractory temporal lobe epilepsy, dipole modelling revisited[END_REF].

In the large family of techniques aimed at solving the inverse problem, methods based on Minimum Norm Estimates (MNE) reconstruct the source activity at each location of the brain volume in a tomographic way. They are particularly suited for the localization of distributed generators of epileptic interictal spikes and their clinical yield was recently adressed [START_REF] Plummer | Clinical utility of distributed source modelling of interictal scalp EEG in focal epilepsy[END_REF]. While these methods proved to be very ecient under certain conditions, they suer from certain drawbacks (see review in [START_REF] Grech | Review on solving the inverse problem in EEG source analysis[END_REF]. In particular, their resolution is blurred when a L 2 -norm of the dipole intensity gradient is used in order to regularize the inverse problem. Other approaches such as Second Order (SO) MUSIC-like algorithms [START_REF] De Munck | Three-way matrix analysis, the MUSIC algorithm and the coupled dipole model[END_REF][START_REF] Mosher | Source localization using Recursively Applied and Projected (RAP) MUSIC[END_REF][START_REF] Mosher | Multiple dipole modeling and localization from spatio-temporal MEG data[END_REF][START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF] search for the best dipoles by exploring the parameter space (location and orientation). These SO methods are very powerful in multiple sources environments under the three following hypotheses: i) the number of sources is less than the number of sensors, ii) the spatial coherence of the noise is known and iii) there is no errors in the forward model. Indeed, SO approaches innite asymptotic precision and resolution whatever the Signal to Noise Ratio (SNR) [START_REF] Kaveh | The statistical performance of the MU-SIC and the minimum norm algorithms in resolving plane waves in noise[END_REF][START_REF] Porat | Analysis of the asymptotic relative eciency of the MUSIC algorithm[END_REF]Stoica andNehorai, 1989, 1990). However, in the particular context of localizing epileptic cerebral sources, none of the three aforementioned assumptions holds. Indeed, i) the sources of EEG epileptic activity are usually distributed over a relatively wide cortical area that involves much more dipoles than the number of scalp electrodes, ii) the EEG background activity behaves at as a Gaussian noise of unknown spatial coherence, and iii) the forward model is corrupted by errors induced by a possible misplacement of electrodes, by incorrect conductivity values for the headmodel, and by the estimation of brain, skull and skin surfaces. Contrarily to SO MUSIC-like methods, HO approaches (Albera et al., 2008;Birot et al., 2010;[START_REF] Chevalier | High-resolution direction nding from higher order statistics: The 2q-MUSIC algorithm[END_REF][START_REF] Chevalier | Higher order direction nding from arrays with diversely polarized antennas: The PD-2q-MUSIC algorithms[END_REF] are asymptotically insensitive to any Gaussian noise, since HO statistics of Gaussian random variables are zero [START_REF] Mccullagh | Tensor Methods in Statistics[END_REF]. They are also shown to be much less sensitive to modelling errors since, from the virtual array theory, the use of HO statistics virtually increases the number of sen-sors. Moreover, HO methods are able to process more independent sources than the number of sensors but this still remains insucient to localize all the possibly-correlated dipoles involved in epilepic distributed sources.

Under certain conditions (evoked potentials or epilepsy for instance), the spatial correlation within activated cortical areas can be exploited in order to summarize the high number of dipoles into a unique entity called Extended Source (ExSo) and consequently to return to an overdetermined mixture.

This property is particularly appropriate in the context of epileptic events. Indeed, several studies showed that these events can be detected on surface EEG recordings because they are likely to involve large cortical areas, in a highly synchronized manner [START_REF] Cosandier-Rimélé | The neuronal sources of EEG: Modeling of simultaneous scalp and intracerebral recordings in epilepsy[END_REF][START_REF] Ebersole | Magnetoencephalography/magnetic source imaging in the assessment of patients with epilepsy[END_REF][START_REF] Gavaret | Electric source imaging in frontal lobe epilepsy[END_REF][START_REF] Merlet | Reliability of dipole models of epileptic spikes[END_REF][START_REF] Tao | Intracranial EEG substrates of scalp EEG interictal spikes[END_REF].

In this paper, we propose a new source localization method specically adapted to such situations where the synchronization of dipole sources within activated areas is high. This novel MUSIC-like method, called 2q-ExSo-MUSIC (q ≥ 1), uses a new principle for the localisation of ExSo's from EEG/MEG data (ExSo-MUSIC) and make use of SO or HO statistics. The ExSo-MUSIC principle i) take advantage of the synchronization of dipoles within each ExSo's and ii) uses a MUSIC-like criterion in order to search for several simple-shaped surfaces which reunion ts the actual ExSo's. The generic method is proposed at any even statistical order and therefore permits the use of higher orders when robustness with respect to the background activity and to modelling errors is needed. In order to quantitatively validate our method, we evaluated its performance at second and fourth order on simulated EEG epileptic spikes obtained with physiologically relevant models of both sources and volume conduction [START_REF] Cosandier-Rimélé | A physiologically plausible spatio-temporal model for depth-EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF], 2008;[START_REF] Wendling | Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals[END_REF] and compared it to the MUSIC algorithm as used in distributed brain source localization [START_REF] Küçükaltun-Yildirim | Taskbased comparison of inverse methods in magnetoencephalography[END_REF], and its extension to Fourth Order (FO).

Materials and methods

Forward model and problem formulation

It is commonly admitted that the signals recorded at the head surface mostly reect the activity of pyramidal cells within the cortex. These activities can be modelled by current dipoles oriented orthogonally to the cortical surface. Consequently, in our study, the source space Ω is restricted to the cortical surface and the source distribution can be seen as a eld of current dipoles over Ω. Each of these dipoles is dened by its location ρ ∈ Ω and its electrical activity s(ρ, t). Under these considerations, the N -channels-EEG/MEG model used throughout this paper has the classical linear form

x(t) = ∑ ρ∈Ω a(ρ)s(ρ, t) + ν i (t) (1)
where {x(t)} is the EEG/MEG data, {ν i (t)} is a Gaussian instrument noise independent of dipole activity {s(ρ, t)}. The noise due to artefacts (muscular, cardiac or ocular) is not represented here since specic techniques have been developed to remove it (Albera et al., 2010;[START_REF] Congedo | On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics[END_REF][START_REF] James | Independent component analysis for biomedical signals[END_REF][START_REF] Vigario | BSS and ICA in neuroinformatics: From current practices to open challenges[END_REF], and denoising is outside the scope of this study. The N -dimensional known leadeld vector a(ρ) transcripts the contribution of a unitary dipole located at ρ and oriented orthogonally to the cortical surface to the set of sensors.

Following the superposition principle, interictal epileptic spikes appearing on EEG data can be considered as the sum of i) epileptic activity rising from epileptic sources and ii) normal background activity distributed all over the source space (see gure 1(a)). Denoting by Θ the origin of epileptic activity, the EEG/MEG vector (1) can be split as

x(t) = ∑ ρ∈Θ a(ρ)s e (ρ, t)
signal due to epileptic activity

+ ∑ ρ∈Ω a(ρ)s b (ρ, t)
signal due to background activity

+ ν i (t) (2)
where s e (ρ, t) is the epileptic activity at location ρ, i.e. the signal of interest and s b (ρ, t) is the background activity at location ρ, i.e. the signal of noninterest. The problem dealt with by the proposed 2q-ExSo-MUSIC method is the identication of the origin Θ of epileptic activity.

Extended source

The conguration of cerebral sources of EEG/MEG signals can be complex depending on space/time parameters under consideration. Nevertheless, is it well admitted that for EEG events in general, and epileptic spikes in particular, to be detectable on scalp EEG, a high synchronization level of intracerebral activity within a large area of cortex is necessary. The 2q-ExSo-MUSIC is built upon this property irrespectively of the spatial conguration of sources. Consequently, an extended source or ExSo (used synonymously to distributed source) is dened in the following as a set of contiguous dipoles (i.e patches) exhibiting a highly synchronized activity. When the activity between two (or more) spatially distinct patches is also highly synchronized, then the presented 2q-ExSo-MUSIC method will consider it as a single ExSo. The entire cortex (including gray regions) generates also background activity. (b) An example of ExSo's conguration. On this schematic representation, the cortical area Θ exhibiting epileptic activity can be partitioned in two extended sources θ 1 (in green) and θ 2 (in red). In θ 1 , two spatially distinct patches are considered. The activity within and between the two patches are synchronized. Therefore θ 1 will be considered as a single ExSo. θ 2 is made of a single patch. The activity within this patch is synchronized but is independent from the activity in θ 1 . Therefore θ 2 is considered as another ExSo.

%

When the activity between two (or more) spatially distinct patches is desynchronized, then the method will regard them as two (or more) ExSos (see gure 1(b)). As the activity is highly synchronized within one ExSo, 2q-ExSo-MUSIC method makes the assumption that these activity is identical.

Thus, the origin of epileptic activity Θ is partitioned into P spatially distinct subsets θ p such that all activity rising from θ p is/are identical. The P sets θ p represent the P extended sources and we have the property: for all 1 ≤ p ≤ P , for all ρ in θ p , {s A (ρ, t)} = {s p (t)}, where {s p (t)} is the activity associated with the pth ExSo. Accordingly, in the EEG/MEG model ( 2), the epileptic activity can be factorized as

x(t) = P ∑ p=1 s p (t) ∑ ρ∈θp a(ρ) + ∑ ρ∈Ω a(ρ)s > (ρ, t) + ν i (t) (3)
The second term of this equation can be seen as a Gaussian noise of unknown spatial coherence. The resulting global noise t) is consequently a Gaussian noise of unknown spatial coherence and independent from the epileptic sources. In order to compact (3), we introduce the ExSo leadeld vector h(θ) = ∑ ρ∈θ a(ρ) that transcripts the contribution of an extended source θ with unitary activity to the scalp electrodes. An ExSo leadeld matrix H(Θ) = [h(θ 1 ), . . . , h(θ P )] and an ExSo activity vector s(t) = [s 1 (t), . . . , s P (t)] T are also introduced. Finally, equation (3) can be written

ν(t) = ∑ ρ∈Ω a(ρ)s > (ρ, t)dρ + ν i (
x(t) = P ∑ p=1 h(θ p )s p (t) + ν(t) = H(Θ)s(t) + ν(t) (4) 2.3. The 2q-ExSo-MUSIC method
The 2q-ExSo-MUSIC method uses the 2q-th order (q ≥ 1) cumulants of the EEG/MEG data {x(t)}. These cumulants are ordered in a (N q × N q ) matrix C 2q, x also called the 2q-th order statistical matrix of the data (see Appendix A for technical details about the estimation of statistical matrices).

In agreement with (4) and using the multilinearity property enjoyed by the cumulants, the 2q-th order statistical matrix of the data can be written as a function of Θ as [START_REF] Chevalier | High-resolution direction nding from higher order statistics: The 2q-MUSIC algorithm[END_REF])

C 2q, x = H(Θ) ⊗q C 2q, s [H(Θ) ⊗q ] T +C 2, ν δ[q-1] (5)
where H(Θ) ⊗q is the Kronecker product of the matrix H(Θ) by itself q -1 times, C 2q, s is the (P q ×P q ) 2q-th order statistical matrix of {s(t)} and C 2, ν is the covariance matrix (second order statistical matrix) of the noise {ν(t)}.

As HO statistics of Gaussian processes are zero, this covariance matrix is weighted by the Kronecker symbol δ[q -1]. We dene the signal subspace

as Span{H(Θ) ⊗q C 2q, s [H(Θ) ⊗q ] T }, which can be directly estimated from the data if q > 1 since Span{H(Θ) ⊗q C 2q, s [H(Θ) ⊗q ] T } = Span{C 2q, x },

and from the data and the noise covariance if

q = 1 since Span{H(Θ)C 2, s [H(Θ)] T } = Span{C 2, x -C 2, ν }.
To be able to identify the origin of epileptic activity, the 2q-ExSo-MUSIC method requires the following assumptions: i) the P processes {s p (t)} are not Gaussian, ii) the dimension of the signal subspace has not reached its maximum (in the worst case, this surely holds if the number P of extended sources is strictly inferior to the number N of electrodes) and iii) the electrode setting has no ambiguities (meaning that each sensor brings its own information to the data). Under these hypotheses the vectors h(θ) ⊗q , for which θ = θ p , belong to the signal subspace. Moreover, for q = 1, all vectors h(θ) such that θ is a union of two or more sets θ p also belong to the signal subspace. For q > 1, the unions that enjoy this property depends on the correlation between the ExSo activity {s p (t)}. But whatever q and the correlation state of the sources, we are able to identify the origin of epileptic activity Θ = ∪ 1≤p≤P θ p by amalgamating all the candidate ExSo's θ such h(θ) ⊗q belongs to the signal subspace.

To check the belongness of h(θ) ⊗q to the signal subspace, a signal sub-

space projector E 2q, s E 2q, s T is built by EigenValue Decomposition (EVD) of C 2q, x -C 2, ν δ[q -1].
More precisely, for q > 1, E 2q, s is a matrix which column vectors are the eigenvector associated with the non-zero eigenvalues of C 2q, x . For q = 1, the covariance matrix of the noise has to be known since E 2, s is then the matrix associated with the non-zero eigenvalues of C 2, x -C 2, ν . Finally, in line with the MUSIC-like methods [START_REF] Chevalier | High-resolution direction nding from higher order statistics: The 2q-MUSIC algorithm[END_REF][START_REF] Schmidt | Multiple emitter location and signal parameter estimation[END_REF], the roots of the criterion

Υ(E 2q, s , h(θ) ⊗q ) = 1 -([h(θ) ⊗q ] T E 2q, s E 2q, s T h(θ) ⊗q )/ h(θ) ⊗q 2
2 , also called metric, allow us to give an estimation, Θ, of Θ:

Θ = ∪ {θ, Υ(h(θ) ⊗q , E 2q, s ) = 0} (6)
However, in practice, the identication of Θ is dicult. Indeed, extended sources can theoretically be a free subset of potentially non-contiguous areas of the source space Ω. Therefore, all the candidate sets θ representing all the possible topologies of extended sources cannot be exhaustively scanned and, as a correlate the comprehensive search of the roots of the metric Υ(h(θ) ⊗q , E 2q, s ) cannot be achieved with a reasonable computational cost. In order to overcome this diculty, we propose a low-cost optimization of the metric. This requires that the complex geometry of θ is parameterized with a relatively small number of parameters. An extended source of any shape can be approximated by the reunion of elementary simple-shaped subsets. Subsets with circular shapes are good candidates since they can be entirely described by only two parameters, say the center and the radius. We propose a non-trivial identication of these subsets based on the thresholding of the metric Υ(E 2q, s , h( θ) ⊗q ) where θ is a circular-shaped candidate subset.

More precisely, the candidate subsets θ for which the metric is lower than a given threshold λ are united in order to estimate simultaneously all the actual extended sources. A schematic representation of this concept is given in gure 2(a). In this situation, the resulting 2q-ExSo-MUSIC becomes

Θ 2q-ExSo-MUSIC (λ) = ∪ { θ, Υ(h( θ) ⊗q , E 2q, s ) ≤ λ} (7)
The validity of such a result is provided in a probabilistic sense in Appendix B. The estimated set of extended sources Θ 2q-ExSo-MUSIC (λ) depends on the threshold λ. This parameter can be adjusted for a good compromise between the specicity and sensitivity of the estimator. As shown in gure 2(b), for low values of λ the estimator has a low sensitivity but a high specicity. This means that the estimated set of extended sources is included in the actual set Θ (true positive) and does not encompass it (false positive). When λ increases, the sensitivity increases while the specicity decreases.

Simulation protocol

The 2q-ExSo-MUSIC method (q = 1, 2) was evaluated on 31 channels EEG data generated from a model developed in our team [START_REF] Cosandier-Rimélé | A physiologically plausible spatio-temporal model for depth-EEG signals recorded with intracerebral electrodes in human partial epilepsy[END_REF], 2008;[START_REF] Wendling | Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals[END_REF], which main features are briey summarized below. In order to construct a source space for our simulations, we built a mesh of the cortical surface from a 3D MRI T1 image of a patient (BrainVisa, SHFJ, Orsay, France). This mesh was composed of source Θ is red. The true estimated parts Θ4-ExSo-MUSIC (λ)∩Θ are beige, and the false estimated parts Θ4-ExSo-MUSIC (λ) -Θ4-ExSo-MUSIC (λ)∩Θ are blue. In practice, as the source space is a triangular mesh of the cortical surface, we cannot construct perfectly circular disks. Instead we consider pseudo-disks, made of one or several triangles and parameterized by the disk area and the position of its center, as dened by the position of the germ triangle. For a mesh with suciently high resolution, any ExSo can be approximated by the union of these pseudo-disks.

the barycenter ρ m of each triangle and oriented orthogonally to the triangle surface, leading to a eld of current dipoles. Our simulation source space was dened by {ρ m , 1 ≤ m ≤ M } and corresponded to a cortical mesh Ω with a resolution of about 1 mm 2 . Time activities {s(ρ m , t), 1 ≤ m ≤ M } of all dipoles of the source space were generated with a realistic neuronal population model [START_REF] Wendling | Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals[END_REF], in which excitation and inhibition parameters of each population can be adjusted to obtain either interictal epileptic spikes or background activity.

With an appropriate setting of coupling parameters bewteen populations, we built a set of highly-correlated interictal spike activities. As illustrated in gure 3(a), these activities were assigned to a set of triangle (actual ExSo) manually outlined with a mesh vizualisation software (Paraview, Kitware Inc., NY, US). Uncorrelated background activities were attributed to the other triangles (gure 3(a)). According to the scenario considered, we varied the size, location, number or conguration of the source(s). From this setup, we built a spatio-temporal source matrix S where the mth line stands for the activity of the mth dipole for K = 9984 temporal samples (39 s with sampling frequency of 256 Hz).

In order to generate the simulated EEG data, we used a realistic head model, that consists in three nested homogenous volume realistically shaping the brain, the skull and the scalp (gure 3(b)). The conductivity of each medium was xed to 0.33 S/m, 0.0082 S/m and 0.33 S/m, respectively [START_REF] Gonçalves | In vivo measurement of the brain and skull resistivities using an EIT-based method and realistic models for the head[END_REF]. The inter-medium surfaces were extracted from the segmentation of a T1-weighted 3D-MRI and meshed by 2440 triangles each Dipoles are placed all over the mesh at the barycenter of each triangle and oriented orthogonal to the triangle surface. Within the source, each dipole is assigned with a spike-like activity. For the whole source, a series of 1000 highly synchronized activities is generated using the neuronal population model (50 are displayed here). In the remaining part of the mesh dipoles are assigned with uncorrelated background activity. (b) Using this source conguration, a realistic head model and a given set of 31 electrodes, EEG data were computed solving numerically the forward problem. In this particular case, epileptic spikes generated from the left superior frontal ExSo are clearly identied on F3, FC1 and Fz.

"

(ASA, ANT, Enschede, Netherlands). From this head model, the forward problem was then numerically calculated for each triangle (ASA, ANT, Enschede, Netherlands), to obtain the leadeld matrix

A = [∆ 1 a(ρ 1 ), . . . , ∆ m a(ρ M )]
where a(ρ m ) is the leadeld vector associated with the m-th triangle and ∆ m the triangle surface. The A matrix gives the contribution of each dipole of the mesh at the level of 31 scalp electrode position (19-20 standard 10-20 electrodes plus additional electrodes at FC1, FC2, FC5, FC6, CP1, CP2, CP5, FT9, FT10, P9 and P10). The spatio-temporel matrix of simulated EEG/MEG data is given by X = AS+V , where V is a zero-mean Gaussian spatio-temporal matrix of the instrumental noise ν i (t).

As previously reported [START_REF] Grova | Evaluation of EEG localization methods using realistic simulations of interictal spikes[END_REF], we used ROC curves [START_REF] Metz | Roc methodology in radiologic imaging[END_REF] to measure the performance of the 2q-ExSo-MUSIC method (q = 1, 2).

This criterion represents

TPF ℓ (λ) = surf{ Θℓ (λ) ∩ Θ} surf{Θ} (8) FPF ℓ (λ) = surf{ Θℓ (λ)} -surf{ Θℓ (λ) ∩ Θ} surf{Ω} -surf{Θ} (9)
where surf{.} is the area of a given cortical region and Θℓ (λ) is the estimated set of extended sources at realization ℓ. In order to display the results in a compact way, only the normalized Area under the ROC Curve (AuC), for FPF from 0 to 0.1, is represented. AuC values range between 0 (the worst result) and 1 (the best result).

In order to analyze the performance of the methods as a function of the spike amplitude on the EEG data, the Mean Spike-to-Background Ratio (MSBR) is introduced as

MSBR = 10 log E[ x s (t) 2 ]/E[ x b (t) 2 ] ( 10
)
where {x s (t)} stands for an EEG epoch containing a spike (negative peak from baseline to baseline) and {x b (t)} is an EEG epoch containing only background. In the EEG generation, the MSBR was tuned by setting the amplitudes of the spikes at the level of the sources S.

3. Results

Inuence of MSBR

In order to study the performance of the methods as a function of the MSBR in simulated EEG data, an ExSo made of a single patch of 1000 mm 2 was placed in the left occipital gyrus. EEG epochs for each considered MSBR are displayed in gure 4(a). The performance of the methods is reported in gure 4(b). For the highest value of MSBR (18.3 dB) 2-ExSo-MUSIC and 4-ExSo-MUSIC showed the best performance. As the MSBR decreased, the AuC of 4-ExSo-MUSIC stayed stable . On the contrary, the performance of the 2-ExSo-MUSIC method decreased quickly from a MSBR equal to 10.7 dB and reached zero for 3.3 dB. The AuC of 2-MUSIC and 4-MUSIC followed 2-

ExSo-MUSIC and 4-ExSo-MUSIC ones, respectively, but with an important shift to the bottom. This result showed i) the better robustness of HO methods with respect to the MSBR and ii) the better accuracy of the 2q-ExSo-MUSIC algorithm with respect to the 2q-MUSIC approach, for a given value of q.

Inuence of the ExSo area

In For all the ExSo areas, the performance of 4-ExSo-MUSIC was always higher than that of other methods. However, ExSo's with a size less than 400 mm 2 (MSBR < 2 dB) remained poorly estimated by 4-ExSo-MUSIC (AuC< 0.3)

and by 4-MUSIC (AuC < 0.1). SO methods did not retrieve any of the ExSo's when their area was inferior to 800 mm 2 (MSBR < 4 dB) while 2-MUSIC, 4-MUSIC and 2-ExSo-MUSIC displayed a stable eciency for ExSo areas ranging from 800 to 1200 mm 2 . The ability of 4-ExSo-MUSIC to estimate the actual ExSo was considerably increased up to AuC values around 0.95.

4-MUSIC was more ecient than SO algorithms in all instances except for ExSo's over 1750 mm 2 for which 2-ExSo-MUSIC slightly exceeded 4-MUSIC.

FInally, 2-MUSIC always showed the poorest performance.

Inuence of the ExSo location

In order to study the inuence of the ExSo location on the eciency of the four algorithms, ExSo's made of single patch were placed successively Illustrative examples of the localization results are given in gure 7 for a xed FPF value of 5.10 -3 . This shows that the localization accuracy of the 4-ExSo-MUSIC method was highly superior to that of the other approaches.

Interestingly, with 4-ExSo-MUSIC, mislocated dipoles were contiguous with the external edge of the actual ExSo, while with classical MUSIC-like methods, mislocated dipoles were more scattered, some of them being not contiguous with the actual source.

Inuence of the number of ExSos and of their conguration

For this study, we have considered three dierent ExSo congurations: i) two spatially distinct patches, located in left inferior and superior temporal gyri, with non-synchronized dipole activities between the two patches (i.e.

two single-patch ExSos, see gure 8(a)), ii) the two same patches but with synchronized dipole activities between patches (i.e. one double-patch ExSo, see gure 8(b)), and iii) the two same patches with synchronized dipole activities between patches and an additional patch, located in the inferior frontal gyrus, which dipole activities are independent from those of the two other patches (i.e. one double-patch ExSo and one single-patch ExSo, see gure

8(c))

.

Table 1 and gure 8 provide the AuC values and illustrate the localization accuracy of the four approaches, respectively. For all three source congurations considered, the 4-ExSo-MUSIC approach was able to retrieve accurately the actual patches, whether interpatch activities have been set synchronous or not. As for the previous study, mislocated dipoles were contiguous with the external edge of the actual ExSos. A smaller area of the actual ExSos was retrieved with the 4-MUSIC approach. This method also mislocated dipoles

Inuence of modeling errors

In practice, some errors may occur in the forward model. The HO MUSIC-like methods have been shown to be more robust with respect to modeling errors than SO-MUSIC-like methods in particular in a multi-sources context (Birot et al., 2010;[START_REF] Chevalier | High-resolution direction nding from higher order statistics: The 2q-MUSIC algorithm[END_REF][START_REF] Chevalier | Higher order direction nding from arrays with diversely polarized antennas: The PD-2q-MUSIC algorithms[END_REF]. We studied the inuence modeling error by displacing scalp electrodes from their original 10-20 position (9(a)). For each electrode position error (0.00 ± 0.00, 2.02 ± 0.65, 6.92 ± 2.60 and 11.58 ± 5. The performance of the FO methods stayed almost stable until an electrode position errors reached 7 mm. On the contrary, AuC values of SO methods signicantly decreased as soon as errors (even of small amplitude) in the electrode position occurred. For a given statistical order, the 2q-ExSo-MUSIC algorithm remained more accurate than the 2q-MUSIC method for all values of electrode position error.

Discussion and Conclusion

The proposed 2q-ExSo-MUSIC method allows for the localization of spatially extended sources from non-invasive electromagnetic recordings, by combining both the use of a novel principle (ExSo-MUSIC) and that of HO The originality of the ExSo principle introduced herein is twofold. Firstly, it relies on the factorization of the multiple dipole activity within the extended source in order to summarize the large number of dipoles into a unique entity and therefore return to an over-determined mixture. This is made possible by the synchronization property of neuronal populations within the source. Indeed, it is very unlikely that uncorrelated or weakly correlated activities within a cortical area could generate signals like epileptic spikes on surface recordings. Indeed, the area and the geometry of the cortical source, along with the synchronisation level of activities within the source are mandatory factors inuencing the way cortical activity is reected at the level of scalp EEG recordings [START_REF] Ebersole | Magnetoencephalography/magnetic source imaging in the assessment of patients with epilepsy[END_REF][START_REF] Gloor | Neuronal generators and the problem of localization in electroencephalography: application of volume conductor theory to electroencephalography[END_REF][START_REF] Gloor | Preoperative electroencephalographic investigation in temporal lobe epilepsy: extracranial and intracranial recordings[END_REF].

Accordingly, simultaneous scalp and intracranial (ECoG) recordings have clearly illustrated that scalp-recognizable spikes are associated with intracranial events that arise in a highly-correlated manner, i.e. synchronously, with very similar shapes and amplitude, on several ECoG channels. Contrarily, no clear paroxysm emerges from background EEG activity when the synchrony and/or the cortical area are insucient (see [START_REF] Tao | Intracranial EEG substrates of scalp EEG interictal spikes[END_REF] for illustrative examples).

In the specic context of non-pathological (such as somatosensory evoked eld/potentials) signals as well, non parametric localization methods such as the Time-Coherent Expansion [START_REF] David | Time-coherent expansion of MEG/EEG cortical sources[END_REF] have exploited uniform current distributions within extended cortical sources. Parametric approaches have also been able to localize such synchronous sources with a small number of equivalent current dipoles [START_REF] Liu | Ecient localization of synchronous EEG source activities using a modied RAP-MUSIC algorithm[END_REF][START_REF] Mosher | Recursive MUSIC: A framework for EEG and MEG source localization[END_REF]. However, these methods took advantage of the high correlation between distant point-like sources and have mostly addressed the question of the interaction across brain areas rather than the geometry of the extended sources. To account for cortical distributed sources, multipolar models (which can be considered as a generalization of the dipolar model) have also been proposed [START_REF] Jerbi | On MEG forward modelling using multipolar expansions[END_REF][START_REF] Nolte | On the calculation of magnetic elds based on multipole modeling of focal biological current sources[END_REF]. One major advantage of the multipolar model is that an extended source of arbitrary current distribution at a given time can be modelled with a relatively few number of parameters. Following this approach, rst-order multipolar models proved to be more ecient than classical dipole models to localize the centroid of extended cortical sources [START_REF] Jerbi | Localization of realistic cortical activity in MEG using current multipoles[END_REF]. Unfortunately, the source geometry is hardly retrieved from multipolar parameters unless assumptions on the source shape are made and high signal-to-noise ratio in data are considered [START_REF] Nolte | Current multipole expansion to estimate lateral extent of neuronal activity: a theoretical analysis[END_REF].

Therefore in order to recover the exact shape of the source a novel parameterization of the source distribution is described. This constitutes the second novel aspect of the ExSo principle, along with the introduction of an ecient and low-cost algorithm for minimizing the metric. In the context of radiocommunications, [START_REF] Valaee | Parametric localization of distributed sources[END_REF] have introduced the Distributed Source Parameter Estimator (DSPE) to estimate the parameters controlling the spatial distribution of the source. However, as is, this method cannot be used in the context of brain source localization of EEG/MEG activities.

Indeed, the denition of the complex source topology would require many parameters and thus would involve a costly optimization of the DSPE criterion.

To overcome this diculty and to reduce the high computation complexity in the resolution of the inverse problem, we chose to make use of the fact that any extended source can be expressed as the union of elementary geometric surfaces like disks, that are described with two parameters only. In practice, only a nite number of circular patches is sucient to precisely recover the considered extended source. In other words, the search for an extended source of unknown geometry is transformed into the search for a nite number of elementary surfaces of only two parameters.

The novelty of ExSo-MUSIC also resides in the use of HO statistics. In the context of radiocommunications, the properties of HO statistics have been successfully exploited in various methods specically dedicated to the estimation of the Direction-Of-Arrival (DOA) for radio waves (see [START_REF] Chevalier | High-resolution DOA estimation with higher order statistics[END_REF], for a review). In particular, these methods were shown to be less sensitive than SO ones, to perturbations induced by a Gaussian noise with an unknown spatial coherence, to modelling errors, or to a small number of sensors. As previously reported in several studies, background EEG activity can reasonably be considered as Gaussian [START_REF] Elul | Gaussian behavior of the electroencephalogram: changes during performance of mental task[END_REF][START_REF] Knuutila | Characterization of brain noise using a high sensitivity 7-channel magnetometer[END_REF][START_REF] Mcewen | Modeling the stationarity and gaussianity of spontaneous electroencephalographic activity[END_REF]. As higher order cumulants of a Gaussian process are null, HO approaches are suited as they are asymptotically insensitive to background EEG activity. Moreover, in real data, background activity is expected to be spatially correlated, this correlation being unknown. Such property penalizes most of SO approaches (such as MUSIC) whereas it is likely to have no impact on HO approaches as the contribution of background activity is suppressed anyways. We recently made use of HO statistics in source localization of EEG/MEG signals extending the concept of deation source localization (RAP-MUSIC, Mosher and Leahy, 1999) to fourth order statistics and showed that their resolution was clearly improved in the case of very close sources, small number of electrodes, or strong background activity (Albera et al., 2008). Similarly, in the present study, the use of HO statistics highly increases the performance of our method, especially when low MSBR's are considered. This is an interesting property in the particular context of epileptic paroxysms, for which the involvement of deep mesial brain sources is common and gives rise to weak amplitude scalp signals compared to the background EEG. Accordingly, in many instances of our simulation study, sources in the cingulate, parahippocampal, or orbitofrontal gyri, as well as in the temporal pole and temporal basal region were retrieved only with the algorithms upgraded to HO statistics.

It is noteworthy that for our simulation study, we assumed that the neuronal populations with highly-correlated activities were contiguous within the extended source. Physiologically other congurations may appear where non-contiguous, yet highly-correlated, populations are involved (distributed sources). This kind of conguration is not unfavourable in our case, since according to our denition, an ExSo is not dened by any spatial constraint but only by the correlation properties of the source activities. Therefore, 2q-ExSo-MUSIC is also able to localize an extended source made up of several non-contiguous small groups of dipoles with correlated activities.

In this study, the performance of 2q-ExSo-MUSIC was evaluated on simulated data. Indeed, the validation of source localization methods implies that the ground truth is known. Unfortunately, this is dicult to achieve when working with real data, and more particularly in the context of partial epilepsies, as the exact localisation of the region(s) responsible for interictal spikes or ictal discharges cannot be dened with certainty. At best, it can be approached only when all functional (including invasive recordings) and structural explorations converge towards the same hypothesis. Alternatively, simulations for which both the geometry and the localisation of sources are dened a priori can assist the validation procedure [START_REF] Grova | Evaluation of EEG localization methods using realistic simulations of interictal spikes[END_REF]. In this study, we focused on the realism of our simulations with regard to the source location and geometry. Firstly, our model of EEG generation combines a model of distributed current dipoles and a model of coupled neuronal populations in order to describe as accurately as possible the spatio-temporal properties of sources of brain electrical activity. Secondly, cortical sources of dierent size were dened in various lobes on a realistic mesh of the cortical surface. These sources were used to generate synthetic data aimed at resembling as much as possible real interictal EEGs recorded in adult partial epilepsies. Thirdly, we considered highly correlated dipole activities within the source to reect the hypersynchronisation property of the epileptic tissue.

In these situations where the characteristics of the cerebral source is known, our method implemented for q = 2 clearly showed higher performance, in terms of TPF/FPF ratio, compared to the 2q-MUSIC (q = 1, 2) algorithm. This was observed for almost all considered ExSo sizes and locations, and was particularly marked for signals with low MSBR (i.e. interictal spikes resulting from the activity of ExSo's of small size or located in deep brain regions). However these simulated results cannot entirely predict the performance of the method in real conditions. More complex simulation scenarios, in particular regarding the number of sources, their spatial conguration and the level of synchronisation between neuronal populations as well as the comparison of our approach with other methods more specially designed for of λ 2 is not satisfying. Therefore, we show in gure B.10, for q = 1, 2, that proposition P2) is almost surely true for λ 2 much less than 1 (λ 2 = 0.12 for q = 1, λ 2 = 0.16 for q = 2) by representing the estimate of the probability P ( Θ2q-ExSo-MUSIC (λ) ∩ Θ = Θ) from the same 27 ExSo's Θ. Interestingly, this probability is equal to one for λ closer to λ 1 than 1. Probability of the estimated extended source Θ2q-ExSo-MUSIC (λ) (q = 1, 2) to be a subset of the actual extended source Θ (referred as P 1 ) and probability of Θ2q-ExSo-MUSIC (λ) to entirely recover Θ (referred as P 2 ).

Table 1 AuC for all four localization methods obtained in the context of three multipatch ExSo congurations.

Figure 1 :

 1 Figure 1: (a) Schematic representation of a cortex exhibiting a region of epileptic activity Θ (in gray).

Figure 2 :

 2 Figure 2: The ExSo-MUSIC principle. (a) Schematic representation of the ExSo-MUSIC principle. When the threshold λ increases, the number of concatenated circular sets θ increases, enlarging the surface of the estimated extended source. (b) Example of result using 4-ExSo-MUSIC method. The actual extended

Figure 3 :

 3 Figure 3: Simulation of scalp EEG spikes. (a) An extended source of 1000 triangles (≈ 10 cm 2 ) was manually dened on a mesh of the cortical surface using Paraview (Paraview, Kitware Inc., NY, US).

  a graph of the mathematical expectation of the True Positive Fraction, E[TPF(λ)], as a function of the mathematical expectation of the False Positive Fraction, E[FPF(λ)]. Estimations, TPF and FPF, of these mathematical expections are computed by averaging the results of 50 monte-carlo realizations. For the ℓ-th realization, TPF ℓ (λ) and FPF ℓ (λ) were expressed as

  order to study the performance of the methods as a function of the area of the source, an ExSo made of a single patch was placed in the left dorsolateral frontal cortex. The size of this patch varied from 50 to 2000 mm 2 (examples are illustrated on g. 5(a)). The AuC for all four methods are displayed on gure 5(b) along with with the corresponding MSBR values (as the patch area has a direct inuence on the MSBR in simulated EEGs).

Figure 4 :Figure 5 :Figure 6 :

 456 Figure 4: Performance of the source localization approaches with respect to the MSBR. (a) In this study, data were simulated from one single-patch ExSo of 1000 mm 2 placed in the left lateral occipital region. Six dierent values of MSBR (from 18.3 to 3.3 dB) are illustrated. (b) For the 2q-ExSo-MUSIC and 2q-MUSIC (q = 1, 2) methods, the Area under the ROC Curve (AuC) is plotted as a function of the MSBR.

Figure 7 :

 7 Figure 7: Examples of localization results of 4-ExSo-MUSIC, 2-ExSo-MUSIC, 4-MUSIC and 2-MUSIC. The threshold was choosen such that FPF ≈ 5.10 -3 (this corresponds to an area of about 350 mm 2 ). (a) 1000 mm 2 ExSo in the left pole. (b) 1000 mm 2 ExSo in the left frontodorsal gyrus. (c) 1000 mm 2 ExSo in the left parahippocampus gyrus. (*) Misestimated parts are located in the cingulate region.

Figure 8 :

 8 Figure 8: Localization results of the 2q-ExSo-MUSIC and 2q-MUSIC (q = 1, 2) methods in the context of three multi-patch congurations. (a) two single-patch ExSo's (i.e. with interpatch independent activity) located in the superior and middle temporal gyrus. (b) one double-patch ExSo (i.e. with interpatch synchronized activity). (c) one double-patch ExSo (i.e. with interpatch synchronized activity) and one single-patch ExSo (i.e. withe activity independent from the rst two) located in the middle frontal gyrus.

  52 mm) EEG signals were simulated from two single-patch ExSos of 1000 mm 2 located in the left anterior frontal and inferior parietal region. For this study inter-patch dipole activities were independent. The MSBR was xed to a high value (15 dB) in order to not disadvantage the SO method with respect to the HO methods. Thus only the inuence of modeling errors is involved. The localization results, expressed as a function of the mean electrode position error, are given in gure 9(b).

  Figure 9: Performance of the methods with respect to the mean electrode position error. (a) The distances between the original (red) and perturbed (yellow) positions were 0.00±0.00, 2.02±0.65, 6.92±2.60 and 11.58±5.52 mm. Data were simulated from two 1000 mm 2 single-patch ExSo's located in the frontal and parietal regions. (b) AuC of the four methods as a function of the mean electrode location error.

  Figure B.10: Probability of the estimated extended source Θ2q-ExSo-MUSIC (λ) (q = 1, 2) to be a subset of the actual extended source Θ (referred as P 1 ) and probability of Θ2q-ExSo-MUSIC (λ) to entirely recover Θ (referred as P 2 ).

  scheme. IEEE Trans Biomed Eng 55 (2), 490501.

  Figure 8 Localization results of the 2q-ExSo-MUSIC and 2q-MUSIC (q = 1, 2) methods in the context of three multi-patch congurations. (a) two singlepatch ExSo's (i.e. with interpatch independent activity) located in the superior and middle temporal gyrus. (b) one double-patch ExSo (i.e. with interpatch synchronized activity). (c) one double-patch ExSo (i.e. with interpatch synchronized activity) and one single-patch ExSo (i.e. with activity independent from the rst two) located in the middle frontal gyrus.

Figure 9

 9 Figure 9

Figure

  Figure B.10

Table 1 :

 1 AuC for all four localization methods obtained in the context of three multipatch

	ExSo congurations.

in a more scattered area. SO approaches displayed quasi-null performance (2-ExSo-MUSIC approach) or failed (2-MUSIC).

Acknowledgement

This work was supported by a fellowship from the Region Bretagne (Program ACOMB, project Solare).

Appendix A. Cumulants and statistical matrix

The cumulants are quantities that entirely describe the statistical behavior of a random variable. For more theoritical informations see (Mc-Cullagh, 1987). In practice, the cumulants are computed from K available samples of the data x(t) ordered in the (N × K) spatio-temporal matrix X.

The covariance and quadricovarince matrices of the data are estimated by

, respectively, where ⊘ is the Khatri-Rao product operator (also refered as column-wise Kronecker product), ⊙ is the Hadamard product operator and 1 N is a N -dimensional vector of ones. The Khatri-Rao product X ⊘ X, also 33 known as column-wise Kronecker product, is equal to [x 1 ⊗x 1 , x 2 ⊗x 2 , . . . , x K ⊗

x K ] where x k is the kth column of X. The Hadamard product is a term by term matricial product.

Appendix B. Validity of the ExSo principle

The validity of our approach can be proved in a probabilistic sense. It is based on both following propositions: P1. For every Θ, it exists at least one threshold λ 1 such that for every λ smaller than or equal to λ 1 , the estimate Θ2q-ExSo-MUSIC (λ) is a subset of Θ.

P2. For every Θ, it exists at least one threshold λ 2 such that for every λ bigger than or equal to λ 2 the estimate Θ2q-ExSo-MUSIC (λ) entirely recovers Θ.

Let's verify that proposition P1) is almost surely true in our EEG/MEG localization context. With this aim, we introduce the probability, as a function of λ, that Θ2q-ExSo-MUSIC (λ) is a subset of Θ, say P ( Θ2q-ExSo-MUSIC (λ) ⊆ Θ).

If it exists a threshold λ 1 such that for all λ ≤ λ 1 , P ( Θ2q-ExSo-MUSIC (λ) ⊆ Θ)

is zero, then proposition P1) will be almost surely true. The probability P ( Θ2q-ExSo-MUSIC (λ) ⊆ Θ) was estimated from 27 ExSo's Θ of 315, 500 and 1000 mm 2 and the result is shown in gure B.10 for q = 1, 2. It appears that this estimated probability reached 0 beyond λ = 0.03 for q = 1 and λ = 0.05 for q = 2, hence the result.

As far as proposition P2) is concerned, it is true because Θ2q-ExSo-MUSIC (1) recovers all the cortex and consequently entirely Θ. Obviously, this value