MEG Source Localization of Spatially Extended Generators of Epileptic Activity: Comparing Entropic and Hierarchical Bayesian Approaches

Localizing the generators of epileptic activity in the brain using Electro-EncephaloGraphy (EEG) or Magneto-EncephaloGraphy (MEG) signals is of particular interest during the pre-surgical investigation of epilepsy. Epileptic discharges can be detectable from background brain activity, provided they are associated with spatially extended generators. Using realistic simulations of epileptic activity, this study evaluates the ability of distributed source localization methods to accurately estimate the location of the generators and their sensitivity to the spatial extent of such generators when using MEG data. Source localization methods based on two types of realistic models have been investigated: (i) brain activity may be modeled using cortical parcels and (ii) brain activity is assumed to be locally smooth within each parcel. A Data Driven Parcellization (DDP) method was used to segment the cortical surface into non-overlapping parcels and diffusion-based spatial priors were used to model local spatial smoothness within parcels. These models were implemented within the Maximum Entropy on the Mean (MEM) and the Hierarchical Bayesian (HB) source localization frameworks. We proposed new methods in this context and compared them with other standard ones using Monte Carlo simulations of realistic MEG data involving sources of several spatial extents and depths. Detection accuracy of each method was quantified using Receiver Operating Characteristic (ROC) analysis and localization error metrics. Our results showed that methods implemented within the MEM framework were sensitive to all spatial extents of the sources ranging from 3 cm2 to 30 cm2, whatever were the number and size of the parcels defining the model. To reach a similar level of accuracy within the HB framework, a model using parcels larger than the size of the sources should be considered.


Introduction
Epilepsy is a neurological disorder characterized by the recurrence of clinical seizures. The state during which the seizure takes place is called the ictal state. In between the seizures, abnormal neuronal discharges, the so-called inter-ictal spikes may take place and usually occur more frequently than the seizures. They are generated by the brain without any clinical manifestations and originate partially from brain regions similar to the ones involved during the seizures, i.e., from the epileptogenic focus. Thus analysis of inter-ictal spikes is widely used as a marker of epilepsy [1][2][3]. The context of the present study is the identification and localization of the epileptogenic focus using these markers, which is crucial during pre-surgical evaluation of epilepsy surgery [4,5].
Epileptic activity originates from abnormal excitability and synchronization of neurons. The large pyramidal neurons of the cortical layer V, which are oriented perpendicularly to the cortical surface of the brain, are the main generators of brain electro-magnetic activity. Magneto-Encephalography (MEG) measures the magnetic fields generated by the neuronal currents, using a helmet of few hundred sensors uniformly distributed around the head [6,7]. This non-invasive modality is used to localize brain regions involved during the generation of epileptic discharges [3,[8][9][10].
The amplitude of MEG signals for physiological brain activity is expected to range from femto-Teslas to pico-Teslas. As mentioned by Huiskamp et al. [11], inter-ictal spikes are spontaneous signals that can have relatively large amplitude (,3 pT in MEG). This implies that epileptic MEG signals are likely to arise from large spatially extended regions of active cortex [12]. A study by Mikuni et al. [13] suggested that MEG can detect epileptiform activity when a cortical area greater than 4 cm 2 is synchronously involved. Comparing MEG spikes with Electro-CorticoGraphy (ECoG) spikes, studies performed by Oishi et al. [14] and Huiskamp et al. [11] showed that MEG sensitivity varies for different regions in the brain. As a result, not only the size of the generators matters, but their location and orientation affect the detection of the MEG epileptic activity [1,13].
The MEG inverse problem of source localization consists in inferring the location of the generators of brain activity from signals detected outside the head [15]. Following a previous study in which we proposed source localization techniques that are sensitive to the spatial extent of the generators of epileptic activity in EEG [16], the present study aims at evaluating the performance of similar methods when applied on MEG data in this context. MEG source localization did show excellent spatial accuracy when validated using invasive studies such as ECoG [14,17], depth electrode recordings [18,19] and post-operative follow-up [20]. While, MEG offers an excellent temporal resolution (few milliseconds), our main objective is to propose a source localization technique that is sensitive to the spatial extent of the underlying generators.
MEG source localization is an ill-posed problem, as it admits no unique solution unless additional information is used to regularize the problem. Such regularizations consist in adding some a priori knowledge or constraints to the problem. For instance, Dale and Sereno [21] introduced anatomical constraints that provided prior information about the sources by fixing the position of the sources along the cortical surface in a distributed source model. This type of constraint makes the inverse problem linear. However, the problem is still under-determined due to the use of few sensors (around 300) to estimate brain activity over a large number of sources (around 4000).
In order to obtain a unique solution, additional constraints in the form of a regularization scheme are required. Minimum Norm Estimate (MNE), which chooses the minimum energy solution [22], and Low Resolution Electromagnetic Tomography (LOR-ETA) [23], which chooses the solution with maximum spatial smoothness are among the first and still very popular regularization techniques proposed to solve this issue. In the present study, we compared two regularization schemes based on the following statistical frameworks: (1) the Maximum Entropy on the Mean (MEM) [16,24,25] and (2) the Hierarchical Bayesian (HB) framework [26,27], because of their flexibility in including prior information or constraint models of different natures.
Based on the rationale of obtaining realistic constraint models describing the generators of epileptic activity, two types of spatial models have been investigated. The first one is the idea that brain activity may be modeled as organized among cortical parcels, that can be active or not, when contributing to specific activity [16,24,28]. The second model is an extension of the spatial smoothness constraint originally proposed in LORETA [23] but locally constrained within cortical parcels as proposed by Trujillo-Barreto et al. [29]. Clustering of the brain activity into nonoverlapping cortical parcels is achieved using a Data Driven Parcellization (DDP) technique similar to the ones described in Amblard et al. (2004) [24], Lapalme et al. (2006) [28], and Grova et al. (2006) [16]. In this study we denoted P(s) the spatial clustering of the whole cortical surface at a spatial scale s, controlling the spatial extent and the total number of parcels.
In order to implement these above-mentioned spatial models, we proposed two new source localization methods within the MEM framework (MEM-s and CMEM-s) and one within the HB framework (COH-s). MEM-s refers to the MEM approach proposed in Grova et al. [16] at a specific clustering scale s, while CMEM-s refers to ''Coherent''-MEM-s, introducing local spatial smoothness within each cortical parcel. On the other hand, within the HB framework, we proposed the ''Coherent at scale s'' (COHs) localization method, modeling the covariance of the sources as a linear combination of source covariance components [27], where each component defines local spatial smoothness over a parcel of P(s). COH-s uses the same spatial model as CMEM-s and has been designed to compare MEM and HB frameworks in similar conditions. In order to assess the ability of these three methods to localize spatially extended epileptogenic generators, we evaluated them within a fully controlled environment using realistic simulations of MEG data. MEM-s, CMEM-s and COH-s were evaluated together with two HB methods, proposed in Friston et al. (2008) [27]: the independent and identically distributed sources (IID) model and the spatially coherent sources (COH) model, as implemented in the SPM8 software (http://www.fil.ion. ucl.ac.uk/spm/software/spm8). We assessed the detection accuracy of all the methods by simulating sources of several spatial extents s e ranging from ,3 cm 2 to 30 cm 2 , and at different cortical depths. Secondly, for the three methods (COH-s, MEM-s and CMEM-s) using the spatial model P(s), we assessed the influence of the spatial clustering scale s on their detection accuracy. We quantified the performance of each method using the area under the ROC curve (AUC) as an index of detection accuracy. We also considered the Mean Square Error (MSE) and minimum geodesic distance (Dmin) as localization error metrics [16].
After introducing the MEG inverse problem using a distributed source model, the definition of the two general spatial models considered in this study is provided: (i) the DDP and (ii) the local spatial smoothness. Then, the MEM framework and the implementation of MEM-s and CMEM-s methods are described, followed by the description of the HB framework and the corresponding methods (COH-s, COH and IID). The evaluation procedure of the source localization methods using realistic simulations is then introduced. Finally the results and a detailed discussion are presented.

MEG Inverse Solution Using Distributed Model
A distributed source model consists of a large number of dipolar sources distributed along the cortical surface. We considered the orientation of each dipole to be fixed perpendicular to the cortical surface. Using this anatomical constraint, the relationship between source amplitudes and MEG measurements is expressed by the following linear model [21]:

M~GJzE ð1Þ
where M is a q|t matrix of the MEG signal measured at q = 275 MEG sensors and t time samples. E models an additive measurement noise (q|t matrix). J is a p|t unknown matrix of the current density along the cortical surface (p,4000: unknown dipolar moment amplitudes). G indicates the q|p lead field matrix obtained by solving the forward problem, by estimating the contribution of each dipolar source on the sensors. However, the inverse problem is still an ill-posed problem as the forward matrix G is under-determined (p..q). There is no unique solution unless a priori model or assumptions regarding the distribution of the sources J are added to regularize the problem. To solve the ill-posed inverse problem, we investigated the relevance of two types of spatial models, DDP and local spatial smoothness, within two regularization frameworks (MEM and HB). In the next sections, we will describe these two spatial models before introducing their implementation within the MEM and the HB frameworks.

Definition of Realistic Spatial Models for Spatially Extended Generators
Data Driven Parcellization (DDP) of the cortical surface. We first assume that brain activity can be organized into functional cortical parcels. Characterizing brain activity, assuming functional homogeneity within brain parcels has proved to be an efficient approach to analyze neuroimaging data, either in EEG/MEG [16,[28][29][30][31], in fMRI [32][33][34] or in multimodal fusion [35,36].
In the present study, we proposed a Data Driven Parcellization (DDP) method performing full parceling of the tessellated cortical surface into non-overlapping parcels (see Figure 1). Such a partition at a specific spatial scale s is denoted by P(s). DDP consists in using partial information from the available data in order to guide this spatial clustering.
The key aspect of DDP lies in the pre-localization of the sources of brain activity using the Multivariate Source Pre-localization (MSP) method [30] followed by a region growing algorithm. MSP is a projection method that estimates a coefficient, which characterizes the possible contribution of each dipolar source to the data. A spatio-temporal extension of the MSP method is described in Appendix S1. From this extension, seed points were iteratively selected among the dipoles showing the highest MSP coefficients. Region growing around each seed points was then iterated until a given spatial neighborhood order s, resulting in a partition of the whole brain into K parcels. This way of choosing the seed points and parceling ensured dipoles contributing to the same underlying generator to be gathered within the same parcel, whereas dipoles contributing to distinct generators to be associated within distinct parcels. A brief description of this DDP technique is provided in Appendix S2.
Defining brain activity in terms of K parcels of functionally homogenous activity (K,,p) aims at better conditioning the under-determined inverse problem, while the inverse method will infer the local source intensity inside each parcel. Local spatial smoothness model. Spatial smoothness model assumes that nearby dipoles are more likely to have similar intensities. In this context, LORETA -originally proposed by Pascual-Marqui et al. [23] -used a discrete Laplacian operator to find the solution with maximum spatial smoothness over a 3D grid.
In order to introduce local spatial smoothness over a geodesic surface, we used the diffusion-based spatial prior proposed by Harrison et al. [37]. Diffusion-based spatial priors are actually constructed using the Green's function of the adjacency matrix defined over the geodesic cortical surface [37,38]. Let us denote A as the (p|p) adjacency matrix of the cortical surface, where A i,i'~Ai',i~1 , if the dipoles i and i' are distinct and directly connected on the mesh, 0 otherwise. The non-zero elements of A define a connection between dipolar sources in the immediate spatial neighborhood.
Let us defineÃ A, the discrete Laplacian over the geodesic surface at the first spatial neighborhood order as: Note that the non-null entries ofÃ A k represent spatial connections between dipolar sources within the k th order spatial neighborhood. We used the spatial smoothness model W introduced in Friston et al. (2008) [27], which is defined by: where s is a parameter that tunes the strength of spatial smoothness, varying between 0 and 1. In equation (3), the upper bound of summation was set to 8 as in Friston et al. (2008) [27]. Note that from equation (3), the spatial smoothness matrix W can be interpreted as a generalization of a discrete Laplacian over a large neighborhood order.

Regularization Techniques
Maximum Entropy on the Mean (MEM) framework. In the MEM framework, we consider the amplitude of the sources J to be estimated as a multivariate random variable of dimension p, with a probability distribution dp(j). The MEM principle aims at estimating the distribution dp p(j) that provides ''maximum uncertainty about missing information carried by the data'' [39], with respect to some reference model assumed on J [24]. Regularization in this framework is introduced by writing the solution in the form of dp(j)~f (j) dn(j), where the reference distribution dn expresses some assumptions on J and f (j) is a ndensity to be found such that it explains the data in average: Among all the distributions dp satisfying the above constraint, the MEM solution dp p~f f dn is the one with maximum n-entropy [16,24]. An interesting property of the MEM approach relies in its inherent flexibility for introducing constraints through the definition of the reference distribution dn. In this study, dn was defined using the parcellization model P(s) assuming brain activity to be described by K cortical parcels showing homogeneous activation state.
Each cortical parcel k is characterized by an activation state S k , describing if the parcel is active (S k~1 ) or not (S k~0 ). Assuming a collection of mutually independent parcels, the global dn was defined as a factorization of the joint probability distribution of the K parcels: with the following mixture model for each parcel: where a k~P rob S k~1 ð Þis the probability of the k th parcel to be active. Multivariate j k denotes the intensities of the p k sources in the k th parcel. d refers to the Dirac distribution allowing to ''shut down'' inactive parcels when S k = 0. N m k ,S k ð Þ is a Gaussian distribution describing the intensities of the k th parcel when active (S k = 1); where m k and S k represents respectively the mean and the covariance of the p k sources within the k th parcel. These parameters will be described in the next section.
The purpose of the present study was to evaluate different initialization of dn in the MEM framework, considering spatial modeling introduced using P(s) and local spatial smoothness within parcels. Once dn is initialized, the MEM solution dp p is obtained through the optimization of a convex function in a q-dimensional space (see Appendix S3 for details). Note that whereas MEM estimation was done iteratively at each time sample, the same clustering model P(s) was used over the whole time window of signal to localize.
Source localization methods within the MEM framework. When incorporating the parcels P(s) through dn in the MEM framework, the first step consists in the definition of the parameters (a k , m k , S k ) of dn k for each parcel (equation 6).
N A spatio-temporal Activation Probability Map (stAPM) was generated (see Appendix S4), by mapping the MSP coefficients of p k sources in the k th parcel (P k ) along time. Therefore, the probability of activation of the parcel P k was initialized at each time sample t as a k t ð Þ~Median i[Pk stAPM(i,t) ð Þ .
N In the reference model (equation 6), we assumed the Gaussian distribution of the active state to be a zero mean distribution; therefore, m k of each parcel was initialized to zero.
N The covariance matrix for each parcel is a time varying matrix S k (t) defined as follows: where g t ð Þ~0:05 1 2 is a scaling factor for the covariance of each parcel P k , estimated using the average of the square of the mean activity provided by the Minimum Norm solutionĴ J MN [22] within the k th parcel. This scaling factor was arbitrarily initialized as 5% of the energy within each parcel. In equation 7, W k s ð Þ is defined as the p k |p k matrix controlling local spatial coherence within the parcel, obtained by selecting the rows and columns of W (s) corresponding to the p k sources of the k th parcel..
Accordingly, under these assumptions, we propose the two following methods: (a) Maximum Entropy on the Mean at a specific cluster scale s (MEM-s) consists in setting s~0 leading to S k t ð Þ~g t ð ÞI p k , where I p k is a p k |p k identity matrix. (b) Coherent-MEM-s at a specific cluster scale s (CMEM-s) consists in setting s=0 leading to S k (t)~g(t)W k (s) T W k (s). We have used s~0:6, as suggested in Friston et al. (2008) [27], to introduce local spatial smoothness in each parcel.
For MEM-s and CMEM-s, we defined the reference distribution with mean m k = 0 with the hypothesis that we do not add much information a priori, since MEM provides inference on the mean of the distribution. On the other hand, we hypothesized that the initialization of the covariance matrix S k (t) for each parcel as 5% of the averaged energy of the Minimum Norm solutionĴ J MN will ensure a proper scale for the intensity of the reference distribution.
Hierarchical Bayesian (HB) framework. Solving the MEG inverse problem within the HB framework offers the advantage of accommodating multiple priors and proposes inference techniques to select the most likely combination of priors using model selection approaches [26,27,[40][41][42][43][44]. We thus chose HB as a key framework in which we could consider similar priors as the ones proposed for MEM thereby allowing an ideal comparison of the two approaches.
HB model allows integrating uncertainties at different levels, modeling the covariance in each level as linear combination of covariance components. The different levels are the sensor noise level and the source noise level.
At the sensor level (1 st level), the relationship between the MEG measurements (M) and the source amplitudes (J) is given by: At the source level (2 nd level), the prior distribution of the source amplitudes (J) is given by: where E 1 and E 2 represents additive random fluctuations in the sensor and source space respectively. The a priori distribution of these additive random noises is a zero mean Gaussian distribution with spatial covariances S sensor and S source and temporal covariance h t , such as: Here h t was modeled as the identity matrix. The sensor spatial covariance S sensor was modeled as: Where Q (1) refers to a spatial covariance component, i.e., the identity matrix here, and l (1) represents the corresponding hyperparameter.
The source spatial covariance S source was modeled as a linear combination of the form: where Q (2)~Q denotes the corresponding hyper-parameters, i[½1,m. The exponential term on l ensures the covariance model to be positive [27]. The hyperparameters were estimated using Restricted Maximum Likelihood (ReML) algorithm, selecting the most relevant linear combination of covariance components (see details in Friston et al. (2002) [42]).

Source
localization methods within the HB framework. In addition to two standard source reconstruction methods (IID and COH) implemented in SPM8 software package [27], we proposed a new method within this HB framework (COH-s).
(a) COH-s: Coherent at a specific cluster scale s: COH-s incorporates spatially smooth extended parcels within the HB framework, thus accounting for the same spatial priors as the ones considered in CMEM-s within the MEM framework. Three types of covariance components were considered: 1) Minimum norm component encoding independent sources Q (2) 1~I p , 2) global spatial smoothness Q (2) 2~W T s ð ÞW(s) and 3) K locally spatially coherent parcels of P(s) as independent covariance components denoted by (C 1 ,C 2 …C K ). C k is a (p|p) block matrix generated using the elements of W T (s)W(s), the block being extracted from the p k row and column indices of the k th parcel, and zero elsewhere. C k thus assumes local spatial smoothness over the whole k th parcel. To summarize, COH-s assumes the following spatial covariance model: where K represents the total number of parcels at a specific spatial clustering scale s. 1~I p (I p being a p dimension identity matrix). This method provides a minimum energy solution, similar to the one originally proposed by Hamalainen and IImoneimi (1994) [22].
(c) Spatially Coherent Sources (COH): This method provides a solution that is spatially smooth, similar to LORETA [23]. It consists in a model with two spatial components modeling respectively independent and spatially coherent sources: Note that COH-s is an extension of COH method using the concept of multiple parcels introduced in the Multiple Sparse Prior method proposed by Friston et al. (2008) [27]. Both COH-s and Multiple Sparse Prior methods are using several regional spatial covariance components. Whereas Multiple Sparse Prior models brain activity as small patches of coherent activity sparsely placed in the left and right hemispheres with a priori maximum variance at the center of the patch, COH-s incorporates spatially smooth extended parcels. In COH-s model, the non-zero terms of the diagonal of C k have a priori the same energy. Multiple Sparse Prior method was designed to localize focal ''sparse'' sources, and was proved to be efficient in cognitive studies [45]. As our objective was to localize spatially extended sources of epileptic activity, we proposed COH-s as a method in the HB framework to be compared with MEM-s and CMEM-s.

Evaluation Using Realistic Simulations
We evaluated the performance of the five above-mentioned source localization methods in their ability to localize spatial extended sources. To perform this validation, we proposed a fully controlled environment to generate realistic simulations of MEG data mimicking the generators of epileptic spikes with different spatial extents, similarly to the evaluation proposed for EEG source localization in Grova et al. (2006) [16]. This section describes the validation dataset and validation metric used for the evaluation.
Ethics statement. Realistic simulations were generated using MEG data obtained from a patient with focal epilepsy showing normal tracing with no epileptic activity. This patient participated as a research subject of the project entitled: ''Application of magnetoencephalography in the assessment of the epileptic focus'' (Dr. E. Kobayashi being the principal investigator for this project). Written informed consent for this study was obtained from the subject as approved by the Research Ethics Committee of the Montreal Neurological Institute and Hospital (MNI/H). At its full board meeting of June 14, 2011, the Research Ethics Board (REB) of the MNI/H has endorsed the review of this project and found this research to be acceptable for continuation at the McGill University Healthcare Centers. The REB of the MNI/H acts in conformity with standards set forth in the (US) Code of Federal Regulations governing human subjects' research and functioning in a manner consistent with internationally accepted principles of good clinical practice.
Validation dataset. The subject we selected to generate our realistic simulations had normal cortical surface segmented from his anatomical Magnetic Resonance Imaging data. This acquisition was done at the MEG center of Université de Montréal on a 275 channels CTF whole-head MEG system. The detection coils used in the system were first order radial gradiometers. The CTF system is equipped with reference sensors using a 3 rd order gradient correction to subtract background interferences. During the acquisition, the head position of the subject was tracked using localization coils placed on three fiducial points (nasion, left and right peri-auricular points).
A high resolution T1 weighted MRI was acquired on the same subject at the MRI center of the Montreal Neurological Institute. Co-registration between MEG sensors position and the anatomical T1-weighted MRI of the subject was obtained in three steps: (i) manual identification of the three fiducial points on the MRI, (ii) digitalization of the position of the fiducials on the head of the subject using a 3D Polhemus localizer and (iii) the rigid geometrical transformation between the MRI's space and the subject's space was obtained by fitting these points using Procrustes method [46]. A realistic head model was obtained by segmenting the surface of the brain from the subject's anatomical T1-weighted MRI [47]. The distributed source model was obtained by segmenting the white/gray matter interface from the MRI using Brainvisa software (BrainVISA: http://www.brainvisa.info). The source model consisted in a realistic 3D mesh of the cortical surface (4203 vertices, 7 mm mesh). The forward matrix G (in equation 1) was computed using the Boundary Element Method (BEM) proposed by Kybic et al. (2006) [48]. A 1-layer BEM model consisting of only the inner skull surface was considered and estimated using OpenMEEG software (OpenMEEG: http:// www-sop.inria.fr/odyssee/software/OpenMEEG/).
Simulation parameters. 100 simulation configurations involving one extended source were generated. The position of each source was selected by choosing a seed point randomly on the cortical surface mesh. The spatial extent of each source was obtained by region growing around the seed following the cortical surface using different spatial neighborhood orders ranging from a source spatial extent s e = 2 (,3 cm 2 ) to s e = 6 (,30 cm 2 ). The amplitude of each vertex of the simulated source was set to 9.5 nA.m, generating an overall maximum signal of 1.5 pT for MEG when all the sources of the cortex were set active. This value has been chosen to mimic realistic amplitude of a typical epileptic spike.
The time course of the simulated sources was the time course of an epileptic spike modeled with three Gamma functions, although only signal around the main peak of the spike was analyzed (about 21 samples around the peak with a sampling rate of 600 Hz). Let us refer to Jtheo as the simulated theoretical current distribution obtained from the spatial distribution of the simulated sources together with the corresponding time course. Noise-free MEG data were then simulated by applying the forward model G to the simulated current density (GJtheo). Realistic physiological noise was extracted from a three minutes segment of MEG background activity acquired on a patient with focal epilepsy showing normal traces without any epileptic discharge. MEG data acquired at 600 Hz were filtered between 0.3 Hz and 70 Hz. Periods with motion and eye blinks were excluded. Each noise-free simulated MEG signals were then corrupted by adding some real MEG background activity. In order to mimic MEG spikes averaging, 128 trials of 700 ms of MEG background activity were manually identified. For each simulation, 20 trials were randomly selected among the 128 trials, averaged and added to the simulated signal. The amplitude of all 128 trials was scaled to ensure a signal-tobackground ratio of 1 (0 dB) for most superficial sources when using reference source amplitude of 9.5 nA.m along a patch of 6 cm 2 . Consequently, the simulation of deep sources resulted in simulated signals with lower amplitude than the superficial sources. Therefore, the Signal-to-Noise Ratio (SNR, defined as the ratio of maximum activity at the peak to the standard deviation of the background activity) of the realistic simulated data varied depending upon the location and extent of the underlying sources.
In order to investigate the influence of the spatial clustering scale s of P(s) for MEM-s, CMEM-s and COH-s, we tested the performance of the methods when varying the spatial clustering scale s from s = 3 (K,200 parcels) to s = 6 (K,40 parcels), for each source spatial extent varying from s e = 2 (,3 cm 2 ) to s e = 6 (,30 cm 2 ), and for each of the 100 random source positions, leading to a total of 4(s)65(s e )6100(configurations)65 methods = 10,000 source localizations.
We also performed the following investigations: 1) to compare the performance of the methods that uses DDP model (MEM-s, CMEM-s and COH-s) when initializing parcels P(s) either with the data of interest or with some background MEG activity and 2) to compare the ability of the methods to localize single spike versus averaged spike data (average of 20 spikes). For these two tests we considered 50 source configurations of spatial extent s e = 3 (,7 cm 2 ) and the methods involving P(s) were localized using a spatial clustering scale of s = 5.
All the simulations were performed with Matlab (R2010a) using the simulation environment Pipeline System for Octave and Matlab (PSOM) [49].
Validation metric. In this section we describe the validation metrics used to evaluate the detection accuracy of the source localization methods presented in section: Regularization techniques. Note that the solution of the inverse problems was estimated and evaluated at one single time sample, at the peak of the spike.
(a) Area Under the ROC Curve (AUC): To assess the detection ability of the different localization methods, we used the Area Under the Receiver Operating Characteristic (ROC) curve [50], denoted by AUC, as a detection accuracy index assessing the sensitivity to the spatial extent of the sources. This metric was adapted by Grova et al. [16] to fit the context of a distributed source model. We chose the AUC index as the validation metric mainly because of the difficulty of providing a valid statistical threshold for all the proposed methods. The AUC index was estimated at the main peak t 0 of the simulated spike. We estimated the energyÊ E of the current density distribution at t 0 , for each localization. To compareÊ E with E ref (energy of the simulated current density distribution),Ê E and E ref were first normalized between 0 and 1 for each dipole i: . We quantified the specificity and sensitivity of the localization method by varying a threshold b between 0 and 1 and considering a dipole i to be active ifÊ E i,t 0 ð Þwb. ROC curves were then obtained by plotting sensitivity (b) against (1-specificity(b)). AUC was finally estimated to assess detection accuracy. In our study, a value of AUC.0.8 was considered to be good detection accuracy, suggesting 80% of detections were accurate.
However, to interpret the area under the ROC curve as a detection accuracy index, one should provide the same number of active and inactive sources to the ROC analysis [16]. Indeed, in the context of distributed source evaluation, the estimation of AUC is biased by the fact that among the p~4203 dipoles of the source model, only few dipoles (p a ) were actually active compared to the large number of inactive dipoles (p{p a ). A more accurate estimation of AUC was obtained by using as many inactive sources as active sources during the evaluation. This was done by randomly selecting p a inactive or fictive sources among the (p{p a ) available either within the immediate spatial neighborhood of the simulated sources (AUC close ) or within far local maxima of the source localization results (AUC far ). The final AUC index was computed as the mean of AUC close and AUC far , thus providing a metric assessing both the ability of the method to focalize the reconstructed activity and the eventual generation of spurious sources far from the simulated one (see [16] for more details).
(b) Mean Square Error (MSE): To assess the ability of the methods to accurately recover the amplitude of the simulated current density (Jtheo), we estimated the mean square error (MSE) [16] between the simulated current amplitude and the reconstructed one over the whole brain at the peak of the spike. Lower MSE values indicate that the method is able to recover the current amplitude with high accuracy.
(c) Minimal geodesic distance to the source (Dmin): To quantify source localization accuracy, we estimated the minimum geodesic distance between the dipolar source showing the global maximum of reconstructed activity source and the closest dipole belonging to the simulated source. This geodesic distance following the circumvolutions of the cortical surface was denoted by Dmin [16]. Solutions for which this global maximum was localized on the wrong hemisphere, Dmin could not be estimated since the surfaces of the two hemispheres were not connected geodesically. Therefore, Dmin was finally set at the largest Dmin value obtained over all source configurations in such cases. A value of Dmin close to 0 indicated that the maximum of reconstructed activity source was found within the simulated source.
In addition, AUC was measured as a function of eccentricity to check for the influence of the depth of the source on detection accuracy. The eccentricity of a simulated source was defined as the distance between the seed point of the spatially extended source to the center of the head, whereby the deepest source have a lower eccentricity value (10 mm) and the most superficial ones have a higher eccentricity value (90 mm). Sources with eccentricity ranging between 40 mm and 60 mm corresponded mainly to mesio-temporal sources and the ones with eccentricity less than 40 mm corresponded to the sub-cortical sources.

Qualitative Assessment
The purpose of this first section is to evaluate qualitatively the performance of three simulations together with the corresponding validation metrics AUC, MSE and Dmin. To visualize the results, we showed the absolute value of the reconstructed activity at the peak of the simulated spike, thresholded upon the level of background activity [51]. Figure 2 illustrates the ability of the five evaluated methods to localize a right occipito-parietal source with an extent of s e = 2(,3 cm 2 ) and an eccentricity of 79 mm (superficial source). Note the AUC values were in agreement with visual inspection. We observed that methods MEM-s and CMEM-s were the most accurate in detecting the spatial extent of the source (AUC.0.90, MSE < 0.70 and Dmin = 0 mm for MEM-s and CMEM-s at s = 3 and 5). IID and COH showed slightly less accurate localization (AUC = 0.88, MSE = 0.98, Dmin = 28.7 mm), probably due to the presence of low amplitude frontal spurious sources. Note that both IID and COH underestimated the spatial extent of the source equally and exhibited very similar solutions. For COH method, ReML model selection actually pushed forward the minimum energy prior over the spatial smoothness prior (cf. ReML estimates for COH, l 2 ð Þ 1~0 .057 and l 2 ð Þ 2~0 in equation (12)). This makes COH interesting when localizing focal sources, as it is able to choose between the minimum energy solution for more focal sources and the spatial smoothness solution for spatially extended sources. Finally, in this specific case, COH-s failed to find the simulated occipito-parietal source (Figure 2b, 2c), as it exhibited a spurious source (Dmin = 120.3 mm at s = 3 and at s = 5 the maximum activity was found on wrong hemisphere) in the deep fronto-mesial region, resulting in poor localization accuracy (AUC = 0.49 and MSE = 891 at s = 3 and AUC = 0.75 and MSE = 6.8 at s = 5).
We also illustrated the impact of the clustering scale s in P(s) (see Figure 2b and The profile of reconstructed intensities using MEM-s is actually similar to typical profiles observed in a minimum norm solution; however, it exhibits a larger contrast over the actual extent of the source. We also noticed a slight overestimation of the extent of the source when using MEM-s and CMEM-s at s = 5. On the other hand, increasing the clustering scale s from 3 to 5 had very little impact on the accuracy of COH-s localization in this example, probably because of the presence of spurious sources. Whereas COH-s did show a larger AUC value of 0.75 for s = 5 when compared to AUC = 0.49 at s = 3, an accurate low intensity source was found at s = 5, but did not pass Otsu's threshold [51], as the intensity of the spurious source was larger. Figure 3 illustrates the ability of the five methods to localize the same right occipito-parietal source but more spatially extended (extent s e = 5 (,15.7 cm 2 ), eccentricity = 79 mm). All methods were able to localize this source, but the spatial extent of the source has been slightly under-estimated. COH was the most accurate in reproducing the source spatial extent with AUC = 0.97, MSE = 0.78 and Dmin = 0 mm. In this example, COH favored the spatial smoothness solution (cf. ReML estimates for COH, l 2 ð Þ 1~0 and l 2 ð Þ 2~0 .28), hence, provided better localization than for the previous example. IID showed less accurate localizations (AUC = 0.84, MSE = 0.88) due to under-estimation of the spatial extent of the simulated source, whereas the maximum of activity was accurately localized (Dmin = 0 mm, Figure 3d). MEM-s and CMEM-s reproduced the source spatial extent with good accuracy at s = 3 and s = 5 (AUC = 0.88 and 0.85 for MEM-s, 0.90 and 0.83 for CMEM-s, with MSE,0.80 and Dmin = 0 mm). CMEM-s was able to detect the local maximum activity of the source following a smooth diffusion along the cortical surface; however, the source spatial extent was slightly under-estimated. MEM-s was able to localize the source but lacked smooth diffusion along the cortical surface. COH-s provided less accurate source localization for both s = 3 and 5 (AUC = 0.78, and 0.71). At s = 3 COH-s was able to detect the source and its spatial extent but also presented a higher intensity deep spurious source (cf. Dmin: maximum located in the wrong hemisphere, MSE = 0.91) (Figure 3b). At s = 5, COH-s detected the source with a low intensity and under-estimated its spatial extent (AUC = 0.71), whereas the maximum was accurately localized (MSE = 0.67, Dmin = 0 mm). Note that the size of the parcels used in COH-s (s = 3 and 5) are smaller than the source spatial extent. Figure 4 illustrates the ability of the five methods to localize a deeper and less extended left orbito-frontal mesial source (spatial extent s e = 3 (,9.   Table 1 reports the medians of AUC values obtained over 100 source configurations for all the five methods and all the source spatial extents s e = 2, 3, 4, 5, and 6. For this comparison, MEM-s, CMEM-s and COH-s were considered using a clustering scale of s = 5. For all the spatial extents, COH, MEM-s and CMEM-s exhibited median AUC values greater than 0.8, indicating overall good detection accuracy with these methods. COH-s showed median AUC values greater than 0.8 for all extents less than s e = 5. Note that for s e $5, COH-s provided low median AUC values (,0.8) when localizing using a clustering scale s = 5. IID showed median AUC values lower than 0.8 for all the extents.

Effect of the Spatial Extent of the Simulated Sources
The above results are presented with more details in Figure 5, which illustrates the distributions of AUC values for all the five methods using boxplot representations for s e = 2, 3, 4, 5 and 6. For this boxplot representation, MEM-s, CMEM-s and COH-s are considered using a clustering scale s = 5. We observed an overall very good accuracy for MEM-s and CMEM-s source localizations (median AUC.0.8 for all spatial extents). However, we noticed that their accuracy in localization decreased slightly when increasing the source spatial extent (Figure 5c and 5d), but still remained among the most accurate methods. For COH, we observed that accuracy in the localization of sources showed slight tendency to increase with the increase in source spatial extent (Figure 5b). For COH-s at s = 5, we observed poor localization accuracy for s e = 5 and 6, and excellent accuracy for smaller sources s e ,5 (Figure 5e). Overall IID was not really sensitive to the spatial extent of the sources, showing the smallest AUC values (Figure 5a). Table 2 reports the median and L1 dispersion (i.e., the average of the absolute deviations from the median) of MSE and Dmin metrics obtained over 100 source configurations for the five methods and all the source spatial extents s e = 2, 3, 4, 5, and 6. For all the spatial extents, MEM-s and CMEM-s provided very similar MSE values (,0.90) indicating an accurate recovery of the source amplitude for all the spatial extents. Similar MSE values were found for COH and IID for all the spatial extents. From all the five methods, only COH-s showed the highest MSE error with a large standard deviation. Overall, MSE was not very informative to compare localization methods in our context, except for detecting important false detection when using COH-s.
The minimum geodesic Dmin distance indicated that in most cases the maximum of reconstructed activity was located within the simulated sources (Dmin = 0), except for some large false detections observed with COH-s at smaller cluster scales s and for the more focal sources s e = 2. Keeping in mind that the average distance between two vertices of the cortical surface was 7 mm, Dmin results at s e = 2 illustrated that for MEM-s, CMEM-s and IID, whenever Dmin was not 0, the maximum activity was found within the 1 st or 2 nd spatial neighborhood order of the source, which is still quite close (Dmin,10 mm).
Whereas Table 2 showed that, except for clear mis-localizations (cf. COH-s for s = 3), all methods performed similarly well according to standard localization metrics (MSE and Dmin), only AUC results (Tables 1 and Figure 5) were able to illustrate that MEM-based methods were indeed sensitive to the spatial extent of the sources. In other words, most methods localized accurately the maximum intensity of the sources, but only methods using a DDP model were sensitive to their spatial extents.  Figure 6 shows the distribution of AUC values for these three methods, at the different clustering scales (s = 3, 4, 5, and 6) and for three spatial extents of the source (s e = 2, 4 and 6). Distributions of AUC values illustrated that the choice of the underlying clustering scale s had no impact on MEM-s and CMEM-s detection accuracy. On the other hand, accuracy of COH-s clearly increased when increasing clustering scale s. Figure 6, Figure 2b and Figure 3b demonstrate the poor performance of COH-s at s = 3. We observed that COH-s provided accurate localization for sources with spatial extent s e lower than the clustering scale s. However, this was not the case for the smallest spatial extent, s e = 2. In this case, even though the clustering scale s = 3 was greater than the spatial extent s e = 2, the overall detection accuracy for this spatial extent was accurate only when s .3 (Figure 2b and Figure 6). Table 2 reports the median and L1 dispersion of MSE and Dmin values for different cluster scales for MEM-s, CMEM-s and COH-s. Similarly to AUC findings, we found that MSE and Dmin remained unaffected by the size of the clusters in MEM-s and CMEM-s. For COH-s we found that MSE and Dmin indicated accurate localization as soon as the cluster scale was larger than 3

Effect of the Depth of the Simulated Sources
We assessed the effect of the depth of sources on detection accuracy by plotting for each method, AUC as a function of eccentricity of the source for source spatial extents s e = 2,3,4,5 and 6. We illustrate the results of only one source spatial extent s e = 4 in Figure 7. The solid lines in Figure 7 correspond to the local moving average of the AUC values for each method. For all the source spatial extents s e , COH-s (for s.s e ), MEM-s, CMEM-s and COH were able to localize most superficial sources (eccentricity .60 mm) with high accuracy (AUC .0.90). For sources with eccentricity ranging between 40 mm and 60 mm, which corresponded mainly to mesio-temporal sources, almost all the methods showed lower localization accuracy than for the superficial sources. Notably, COH-s at s = 6, MEM-s and CMEM-s still provided relatively good localization accuracy for these mesiotemporal sources (most AUC values .0.8). However, none of the methods (except COH-s at s = 6) were able to localize accurately deep sources (eccentricity,40 mm). COH-s at s = 6 performed better than other methods when localizing deep sources, although AUC was greater than 0.8 for s e = 4 only (Figure 7). IID exhibited poor localization accuracy for all the deep and mesio-temporal sources, although it provided better localization accuracy for superficial sources (AUC,0.8). Figure 7b clearly demonstrated that when increasing the clustering scale s, COH-s showed an increase in detection accuracy for all sources; this was observed for the other source spatial extents too. From figure 7c and 7d, we observed that MEM-s and CMEM-s showed no notable influence of the clustering scale s on detection accuracy, which was observed for the other source spatial extents too.
This analysis also confirmed that most of the low AUC values considered as outliers in boxplot distributions, presented in   Effect of Parcellization using Data of Interest or Background Activity Figure 8a illustrates the effect of initializing the parcellization P(s) (cf. Appendix S1) when using some MEG background activity or the data of interest, whereas only data of interest was considered to initialize the other parameters (a k , m k , S k ). This investigation was performed using 50 source configurations with source spatial extent s e = 3 (,7 cm 2 ) and the methods involving P(s) were localized using a spatial clustering scale of s = 5. As localization of deeper sources would provide lower AUC values (see section: Effect of the depth of the simulated sources), in order not to bias our comparisons, we decided to consider only the most superficial sources with an eccentricity greater than 70 mm for this section i.e., 22 out of the 50 simulated sources were thus considered.
We can see that there was no impact on the detection accuracy of MEM-s, CMEM-s and COH-s when using the data of interest or MEG background activity to initialize the DDP model (Figure 8a), suggesting the importance of using a parcellization, but the accuracy of such a parcellization was not an issue. Figure 8b illustrates the comparison between detection accuracy of all the five methods when localizing single spike versus the average of 20 spikes. For this investigation, we used the same superficial source configurations as for the previous evaluation (section: Effect of Parcellization using data of interest or background activity). As expected, results showed a higher accuracy for the averaged data than the single spike data. We can also notice that the percentage of single spike simulations that were localized with a good detection accuracy (i.e., AUC.0.8) were respectively: CMEM-s = 37%, MEM-s = 18%, COH-s = 9%, COH = 4%, IID = 4%. These results suggest the robustness of both MEM-s and CMEM-s methods to low SNR conditions.

Discussion
We have presented an evaluation of source localization methods for their ability to localize spatially extended sources of epileptic activity by incorporating realistic spatial models. We proposed three new source localization techniques that can detect the location of the sources with a good sensitivity to their spatial extent when using MEG data: MEM-s, CMEM-s and COH-s.
Standard localization metrics (Dmin and MSE) demonstrated overall good accuracy for most of the evaluated methods (except Table 2. Median (Med) and L1 Dispersion (Disp) of MSE and Dmin over 100 source configurations for all five methods, all five spatial extents s e = 2,3,4,5,6 and all four clustering scale s = 3,4,5,6 in MEM-s, CMEM-s and COH-s.

Methods
Metrics COH-s when s,s e and s e = 2), suggesting that the maximum of the activity was accurately localized in most cases. More importantly, only AUC metric was able to characterize better the ability of some methods to recover accurately the spatial extent of the sources, suggesting AUC to be a more appropriate metric in this context [16]. Whereas it is generally accepted that the minimum norm model (IID) is suitable for localizing the maximum of the activity [22], our results showed that IID was not sensitive to the spatial extent of the generators. Our findings are in agreement with the study of Ding (2009) [62], which demonstrated that a minimum norm model failed to recover the continuous cortical distribution of extended sources. In COH, with two variance components accounting for minimum norm and spatial smoothness models, ReML infers the most relevant combination of these two priors to fit the data. This method was able to localize accurately spatially extended sources. We also noticed that for focal sources, COH was able to choose the minimum energy solution over the spatial smoothness solution using ReML, which indeed is a very interesting property (see Figure 2d and 4d). COH-s seemed an appropriate method for estimating spatially extended sources; the full brain parcellization of the cortical surface into non-overlapping parcels with local smoothness proved indeed to be useful. This suggested that modeling source covariance as a linear combination of covariance components and inference using ReML is an interesting methodological approach, offering sufficient flexibility in the definition of the a priori model. However, COH-s showed some instabilities in the form of spurious sources mainly in cases where the clustering scale s was smaller than the spatial extent s e , as well as when dealing with very focal sources (s e = 2). The method was able to find accurately the sources in most conditions but the presence of some spurious sources far from the main generators led to lower AUC and larger MSE and Dmin values. From the three proposed methods using the DDP spatial model (MEM-s, CMEM-s and COH-s), MEM-s and CMEM-s were the most accurate and stable in localizing the spatially extended sources, especially since they provided accurate results whatever was the underlying clustering scale s.
All approaches, except COH, presented a loss of performance when increasing the spatial extent of the source. Although, a slight decrease in the localization accuracy of MEM-s and CMEM-s was noticed for more extended sources (Figure 5c and 5d), their overall accuracy remained better than for the other methods. This decrease in the accuracy could be explained by the fact that most of the spatially extended sources (s e = 5 and 6) involved multiple sulci and gyri including substantial radial and deep components leading to more cancellation effect of the MEG signals [63]. This resulted in low amplitude MEG signal for specific regions and high amplitude for other regions (superficial and tangential components). Indeed, a spatial extension of order 5 or 6 from a superficial seed will have more chance of spreading along some mesial or basal aspects of the cortical surface, thus involving generators that are more difficult to localize. Moreover, the spatial extent s e = 5 (,18 cm 2 ) and s e = 6 (,30 cm 2 ) are less realistic (or more rare) for an expected extended epileptiform activity, whereas the extent of 4 (,11 cm 2 ) or less are the more realistically expected extents according to Huiskamp et al. (2010) [11].

Localization of Deep Sources
Detection of deep sources is a difficult issue since deep generators will generate very low amplitude MEG data on the scalp, almost undetectable except under specific conditions. This usually requires lots of averaging to increase the SNR [64]. Within our simulation framework, we considered both deep and superficial sources, mimicking realistic SNR conditions in both cases. All the evaluated methods (except IID) accurately localized superficial sources (Figure 7), whereas they all demonstrated poor performance when localizing sub-cortical sources (eccentricity,40 mm) and some variability in accurate localization of mesio-temporal sources (40 cm,eccentricity,60 cm). Note that for eccentricity between 40 mm and 60 mm, ''mesio-temporal sources'' referred mainly to sources located on the mesial aspects of the temporal pole, rather than in the hippocampus per se. COH-s showed low accuracy when using small parcels (s,s e ), even for superficial sources. On the other hand, when using larger parcels (s.s e ) COH-s provided accurate results for superficial sources and most mesio-temporal sources. Overall, COH-s at s = 6 performed better than other methods when localizing sub-cortical sources in most cases of s e . It was not surprising that IID would give poor detection for deep sources because superficial dipolar sources with smaller magnitudes are favored by the minimum norm constraint. Similar trends towards most superficial sources were observed on all the evaluated methods. This bias towards the superficial sources has not been addressed in our study but few studies have been carried out to compensate such a bias [23,29,65]. We did not use depth weighting in any of the methods studied in this paper, in order to evaluate all the methods in the same context. However, it should be noted that even if no depth weighting was used, COH-s (for s.s e ), MEM-s and CMEM-s and were able to localize some deeper sources (Figure 7b, 7c and 7d), more likely because the underlying model was putting forward the involvement of these parcels in the solution, similar to the Bayesian Model Averaging (BMA) approach [29].

Spatial Models: Data Driven Parcellization and Local Smoothness Prior
The two realistic spatial models considered in this study, i.e. data driven parcellization P(s) and local spatial smoothness W, have been implemented within both the MEM and HB frameworks to model the activity of the underlying generators of epileptic discharges.
The first model assumed brain activity to be organized into several spatial clusters P(s) using DDP of the brain activity along the tessellated cortical surface. Few studies demonstrated how introducing parceling of the brain (obtained from some anatomical atlases) was quite useful to better condition the inverse problem [29,31,36,66]. Similarly, Lapalme et al. [28] used a data driven parceling technique to partition the whole cortical surface into functionally homogenous parcels. This approach is similar to the one used in the present study. The performance of MEM-s, CMEM-s and COH-s confirms the usefulness of DDP of the whole cortical surface in detecting extended sources. Although DDP was crucial to regularize the inverse problem, the accuracy of the underlying DDP was not required to obtain an accurate MEM-s or CMEM-s localizations ( Figure 6). In case of COH-s, while the size of the underlying DDP was crucial in the Hierarchical Bayesian approach, the type of data used for the parcellization did not alter the COH-s solutions. This was also confirmed with the assessment of the performance of MEM-s, CMEM-s and COH-s when the parcellization was obtained from MEG physiological background activity, instead of using the data of interest. We did not see any difference in their localization accuracy with the type of data used. We applied one method of parcellization here, but other methods could have been considered [29,31,36,66], with probably similar level of localization accuracy.
However, MEM-s was unable to accurately recover the spatial smoothness of the source along its extent. This led us to incorporate the second model W, i.e., local spatial smoothness within the parcels, giving rise to CMEM-s method. W models the local spatial smoothness of the distribution of the source activity along the cortical surface. The very first idea of spatial smoothness prior was LORETA [23] that allows the smooth reconstruction of cortical sources at a low spatial resolution but does not generally reflect the focal nature of most cortical activations. W used in this study introduces spatial smoothness following the geodesic surface [37], thus creating local spatial smoothness within each parcel. As a result, CMEM-s was able to accurately recover the spatial smoothness of activity within the extended sources. Similarly, COH-s, which also incorporated local spatial smoothness within the parcels, was also able to recover the spatial smoothness along the source extent as soon as the size of the clustering scale s was larger than the extent of the source s e .

Comparison between MEM and HB Frameworks
Two statistical regularization schemes, the MEM and the HB frameworks, were compared in this study.
COH-s method, in which we incorporated the parcellization model P(s) and the local spatial smoothness prior, was proposed to be the equivalent of CMEM-s method. It incorporates the same constraints as the CMEM-s method, but it uses the ReML algorithm within the HB framework to estimate the solution. Our findings show a good concordance between the MEM and HB frameworks when comparing the CMEM-s and COH-s for their detection accuracy, suggesting that both frameworks offer sufficient flexibility to build efficient source localization methods, especially in the context of MEG epileptic data.
In addition, we have tested the influence of the clustering scale s on these methods. We showed that clustering scale had no impact on detection accuracy for both MEM-s and CMEM-s, whereas COH-s method was indeed very sensitive to the clustering scale. The distribution of AUC values for COH-s shows a tendency of increase in detection accuracy when increasing the clustering scale. Indeed, we found that COH-s provided poor detection accuracy when using parcels smaller than the source spatial extent (i.e., s,s e ), whereas for s.s e it provided overall good detection accuracy (see Figure 6 and 7). We could consider the highest clustering scale s = 6 to accurately localize the range of sources' extents simulated in this study. However, our results still suggest that this scale parameter should be tuned from the data when localizing clinical data, as the underlying spatial extent of the generator could not be predicted.
MEM-s and CMEM-s provided very accurate results for any evaluated spatial clustering scale s. This is an important result, suggesting that MEM regularization is able to adapt the number of active parcels, whatever is the spatial scale of the clustering. In order to localize a spatially extended source as accurately as possible, MEM is able to ''switch on'' several parcels when using a lower clustering scale (small parcels) or only few of them when using a larger clustering scale (large parcels). Once the parcels have been identified as active, our results demonstrated that MEM inference is still able to create some local contrasts of dipole intensities within the active parcels, leading to the ability of localizing sources of different spatial extents. The regularization process was a bit different when using COH-s with ReML, as the hyper-parameters of the source covariance components (i.e., the parcels for COH-s) are first estimated through an Automatic Relevance Determination (ARD) scheme. Then, once the covariance model and its ''weights'' are estimated, sources are estimated using a regularized pseudo-inverse method. The diffusion-weighted prior will then push forward a spatially smooth solution over the selected parcels. When using smaller parcels with COH-s (s,s e ) we observed that some covariance components often associated with focal sources were then falsely enhanced, leading to the occurrence of spurious sources, while the main source was still found. In addition, it is possible that when using local spatial smoothness prior within the parcels COH-s was less efficient than CMEM-s in creating contrast of dipole intensities within the parcel. However, this was out of the scope of this study and will require further investigations. Since Litvak and Friston (2008) suggested that using a Greedy Search scheme instead of ARD scheme improved the performance of Multiple Sparse Prior method [67], we can expect that using a Greedy Search for COHs model could have similar impact on localization accuracy and stability. However, such implementation was out of the scope of the present study and will be considered for future publications. Overall, MEM framework and HB framework using ReML are two very efficient and sufficiently flexible regularization approach-es that can be adapted to localize spatially extended sources. An important aspect of ReML-based approaches is that quantitative model comparison using free energy [26] is possible and may be quite useful when testing the relevance of several models, as demonstrated in Henson et al. (2009) [45].
Note that, in this work, we compared the performances of the five methods using the same forward model for both source simulation and source localization. This is a standard approach when assessing the performances of source localization methods [16,27,29,31,58,68]. This is a best-case scenario and a decrease in performance should be expected when applied on real data. However, we investigated the detection accuracy of all the methods when using either the same forward model for both simulation and localization (1 layer BEM model) or a BEM model for simulation and an analytical spherical model for localization (results not shown). All the inverse solvers then outperformed in the same proportion when different models were used, without modifying the overall relative performances. We can conclude that our evaluation did not bias the results of any particular method, allowing appropriate comparison between methods, which was the purpose of this study.

Practical Application and Future Work
In distributed source modeling, the cortical surface constraint is defined from large cortical assemblies of pyramidal cells organized orthogonally to the grey-white matter interface. Most distributed methods adopt this constraint with either restricting the orientation to be perpendicular [21] or allowing some deviation from the surface normal [69].   [69] and Henson et al. (2009) [45] showed that the use of loose orientation dipolar sources increased the localization accuracy of minimum norm methods. In our distributed model, although we fixed the orientation normal to the surface, the loose orientation constraint could be easily adapted and would impact all the methods similarly; but this falls out of scope of the present study.
When dealing with real data, it is important to study the impact of the quality of segmentation and resolution of the cortical surface on the source reconstruction, and especially in pathological conditions. On the other hand, Henson et al. (2009) [45] showed that the nature of the cortical surface (obtained from individual MRI or warping the individual's MRI to a template MRI) had only minimal impact on MEG source reconstruction, as soon as the overall head model was transformed to take into account the actual head shape of the subject. They also showed that the resolution of the cortical surface (3000 sources or 7000 sources) had no reliable effect on the Minimum norm solution; whereas the 7000 sources surface resulted in better accuracy for their proposed Multiple Sparse Prior model. Whereas these are important issues to cope with, they will have similar influence on all the methods involved in the present evaluation study.
Our simulation paradigm was a spatial validation, studying detection accuracy only at the main peak of the simulated spike. In our future work, we plan to assess the spatio-temporal features of the sources by simulating different propagation patterns of epileptic discharges using models such as the extended source model developed by Cosandier-Rimélé et al. [53,70], which describes both the spatial distribution and the temporal dynamics of neuronal population. The code for MEM-s and CMEM-s methods has been recently implemented as a toolbox in Brainstorm Package [71] (http://neuroimage.usc.edu/ brainstorm/), to be released soon.

Conclusion
We have proposed three new methods (MEM-s, CMEM-s and COH-s) and evaluated their performance when localizing spatially extended generators of epileptic discharges using MEG data. We demonstrated that modeling brain activity using Data Driven Parcellization of the cortical surface and applying local smoothness within each parcel is particularly relevant to localize sources together with their spatial extent. Both MEM and HB frameworks are sufficiently flexible to allow the implementation of such spatial models. The present study is in agreement with the good performance of MEM we previously demonstrated on EEG data, although we added the evaluation of several other parameters such as a larger range of spatial extents and depths of the sources, as well as the scale of the spatial clustering.

Supporting Information
Appendix S1 Spatio-temporal extension of the Multivariate Source Pre-localization (MSP).