Accuracy and spatial properties of distributed magnetic source imaging techniques in the investigation of focal epilepsy patients

Abstract Source localization of interictal epileptiform discharges (IEDs) is clinically useful in the presurgical workup of epilepsy patients. We aimed to compare the performance of four different distributed magnetic source imaging (dMSI) approaches: Minimum norm estimate (MNE), dynamic statistical parametric mapping (dSPM), standardized low‐resolution electromagnetic tomography (sLORETA), and coherent maximum entropy on the mean (cMEM). We also evaluated whether a simple average of maps obtained from multiple inverse solutions (Ave) can improve localization accuracy. We analyzed dMSI of 206 IEDs derived from magnetoencephalography recordings in 28 focal epilepsy patients who had a well‐defined focus determined through intracranial EEG (iEEG), epileptogenic MRI lesions or surgical resection. dMSI accuracy and spatial properties were quantitatively estimated as: (a) distance from the epilepsy focus, (b) reproducibility, (c) spatial dispersion (SD), (d) map extension, and (e) effect of thresholding on map properties. Clinical performance was excellent for all methods (median distance from the focus MNE = 2.4 mm; sLORETA = 3.5 mm; cMEM = 3.5 mm; dSPM = 6.8 mm, Ave = 0 mm). Ave showed the lowest distance between the map maximum and epilepsy focus (Dmin lower than cMEM, MNE, and dSPM, p = .021, p = .008, p < .001, respectively). cMEM showed the best spatial features, with lowest SD outside the focus (SD lower than all other methods, p < .001 consistently) and high contrast between the generator and surrounding regions. The average map Ave provided the best localization accuracy, whereas cMEM exhibited the lowest amount of spurious distant activity. dMSI techniques have the potential to significantly improve identification of iEEG targets and to guide surgical planning, especially when multiple methods are combined.


Magnetoencephalography (MEG) is a valuable neuroimaging clinical and
research tool to localize the epileptic focus in patients affected by focal epilepsy. Source localization of interictal epileptiform discharges (IEDs) detected with MEG can guide intracranial EEG (iEEG) targeting and help neurosurgical planning for cortical resection in drug-resistant patients (Knowlton et al., 2006;Ryvlin, Cross, & Rheims, 2014;Stefan et al., 2003).
ECD-based localization has been widely used and validated in clinical settings in North America as it is the unique technique approved by clinical US guidelines and Food and Drug Administration (Bagic et al., 2011;Barth, Sutherling, Engel, & Beatty, 1982;Bast et al., 2004;Knowlton et al., 1997;Stefan et al., 2003).
On the other hand, dMSI techniques are widely used in research settings, while gaining popularity as additional or even the sole method in many tertiary neuroimaging facilities for clinical applications, especially in Europe (Mouthaan et al., 2016;Mouthaan et al., 2019). An important development of dMSI relies on the usage of realistic anatomical constraints, distributing a set of dipolar sources along the cortical surface, where dipole orientation is set to be normal to the underlying cortical mesh. This was made possible through the early use of advanced segmentation techniques on individual anatomical MRI data (Dale & Sereno, 1993). dMSI allows a better assessment of the spatial properties of the generator (Dale & Sereno, 1993;Darvas, Pantazis, Kucukaltun-Yildirim, & Leahy, 2004), offering robustness to low signal to noise ratio data (Tanaka & Stufflebeam, 2013), and providing maps of cortical activations that are clinically relevant. In general, dMSI approaches are less operator dependent and rely on the assumption that the generators of MEG-recorded activity are distributed over a spatial mesh and spatially extended (Tao, Baldwin, Hawes-Ebersole, & Ebersole, 2007;von Ellenrieder, Beltrachini, Muravchik, & Gotman, 2014;von Ellenrieder, Beltrachini, Perucca, & Gotman, 2014;Pellegrino et al., 2019;Raffin et al., 2015).
In the present study, we report MEG source localization performance on patients with a well-characterized focus based on noninvasive and invasive evaluation of their epilepsy. We quantitatively compared the performance of dMSI approaches and tested whether combining multiple inverse solutions can improve localization accuracy. We compared four dMSI techniques most commonly used in presurgical evaluation and which are freely available and easy to translate into clinical settings: (a) cMEM (Chowdhury et al., 2013); (b) minimum norm estimate (MNE) (Hämäläinen & Ilmoniemi, 1994); (c) standardized low-resolution electromagnetic tomography (sLORETA) (Pascual-Marqui, Esslen, Kochi, & Lehmann, 2002); and dynamical statistical parametric mapping (dSPM) (Dale et al., 2000).

| Patients
The study was approved by the Research Ethics Board of the Montreal Neurological Institute and Hospital-McGill University Health Center. All patients signed a written informed consent prior to enrollment. All experimental procedures were performed in agreement with the Declaration of Helsinki and its later amendments. We retrospectively reviewed all focal epilepsy patients scanned over a period of 9 years between 2006 and 2015 and included those who underwent presurgical assessment leading to accurate identification of the epilepsy focus. The identification of the epileptic focus was deemed accurate and the patient included when at least one of the following was satisfied: (a) patient underwent surgery and became seizure free, (b) invasive EEG capturing interictal and/or ictal activity, and (c) epileptogenic MRI lesion concordant with scalp EEG findings. We excluded patients whose MEG had magnetization artifact, patients with extensive brain lesions that could affect the estimation of the forward model and patients affected by mesial temporal lobe epilepsy, as this condition does not require a presurgical MEG assessment (Pellegrino, Hedrich, et al., 2018). We also excluded patients having less than five spikes per study to comply with American Clinical Magnetoencephalography Society guidelines and to avoid a too low signal to noise ratio. Clinical data are reported in Table 1  to 600 Hz (whenever necessary) (Pellegrino, Maran, et al., 2018;Pellegrino et al., 2019). MEG IEDs were visually identified and marked at their peak by trained neurophysiologists (G. P. and E. K.). MEG signal was then segmented into epochs of 2 s around IED peak (−1 to +1 s). IEDs were classified according to topographical distribution of the magnetic field and their morphology. IEDs occurring at the time of artifacts or EKG QRS complex were excluded. Epochs of the same type belonging to the same run were averaged. Each averaged IED group corresponded to a source imaging "study" (Heers et al., 2016;Pellegrino, Hedrich, et al., 2016;Pellegrino, Hedrich, et al., 2018).

| Forward model
Individual MRIs were processed using the FreeSurfer toolbox (Dale, Fischl, & Sereno, 1999). Skull surface was reconstructed from individual MRI in Brainstorm. MEG-MRI coregistration was ensured using a procedure of surface fitting with a rigid geometrical transformation (three rotations, three translations) between the head shape segmented from MRI and the head fiducials and the head points previously digitized. We selected the mesh of the "mid" layer equidistant from the white/gray matter interface and the pial matter as source space for source imaging and downsampled it to 8,000 cortical vertices. The distributed source model consisted in dipolar sources on every vertex of the mesh, oriented perpendicularly to the cortical surface. The forward model assessing the contribution of every dipolar source to MEG sensors was built with a 1-layer boundary element method (BEM) method (Kybic, Clerc, Faugeras, Keriven, & Papadopoulo, 2006) as implemented in the OpenMEEG toolbox (Gramfort, Papadopoulo, Olivi, & Clerc, 2010). This BEM model was computed using the inner-skull surface segmented by Brainstorm for every patient and cortical conductivity was set to 0.33 S/m.

| dMSI methods
Source maps were estimated at the time of the peak of the average IED. Noise-covariance was modeled using a 1 s baseline without any visually identified IEDs.
We computed the following inverse solutions applying default Brainstorm parameters: 1 MNE (Hämäläinen & Ilmoniemi, 1994): This is the "depth-weighted" linear L2-MNE current density (also known as wMNE) which imposes a constraint of minimum energy on the source map, which is equivalent to Tikhonov regularization to solve an underdetermined linear problem (i.e., MNE minimizes the L2-norm of the current distribution).
2 sLORETA (Pascual-Marqui et al., 2002) is a standardized version of MNE, which normalizes the MNE solution by its resolution matrix, allowing sLORETA to be sensitive to deeper generators and to provide no localization error in the absence of measurement noise.
3 cMEM (Amblard et al., 2004;Chowdhury et al., 2013): This is a nonlinear probabilistic Bayesian approach relying on maximizing relative entropy. cMEM allows "switching off" the cortical parcels that do not contribute to the solution, while preserving the ability to create a contrast of current intensities within the active parcels and is sensitive to the spatial extent of the generator (Chowdhury et al., 2013;Chowdhury et al., 2016;Grova et al., 2016;Hedrich, Pellegrino, Kobayashi, Lina, & Grova, 2017;Heers et al., 2016). 4 dSPM (Dale et al., 2000) is another noise-normalized version of MNE using a source noise covariance matrix and the whitened forward solution to obtain time-varying statistical maps of neuronal activity.
In order to assess whether the combination of multiple source imaging techniques can achieve better accuracy, we computed an   original maps were rescaled so that the amplitude at each cortical vertex ranged between 0 (no activity) and 1 (maximum activity).

The result of dMSI source imaging is a map (dMSI map)
where each vertex of the cortical surface is characterized by an intensity indicating its contribution to the signal recorded at sensor level.
Finally, we have also performed a standard analysis with the ECD technique (Bagic et al., 2011), which is reported in the supplementary material (Supplementary Material, Dipole Analysis), as the focus of this study was on dMSI.

| Definition of the epileptic focus
Prior to MEG data analysis, the epileptic focus was manually drawn along the cortical surface based on the available clinical information for each patient. This information is included in Table 1, and consisted of (in order of priority; not all factors were available for every patient): resected region, ictal and interictal invasive EEG findings, visible lesion in the MRI, ictal and interictal scalp EEG findings. The identification of the region was never performed on scalp EEG findings alone, and an epileptogenic MRI lesion or invasive recordings were always required. In patients who became seizure free after surgery, the epileptic region corresponded to the surgical cavity and an ad hoc analysis was performed for this small group. The location of the epileptic focus was then independently identified by two epileptologists (E. K. and G. P.) over the individual cortical surface, and a consensus was reached on its extension upon discussion, blind to source localization results, following the same methodology proposed in our previous studies (Pellegrino, Hedrich, et al., 2018;von Ellenrieder et al., 2016).

| Assessment of performance and spatial properties of dMSI methods
The following quantitative metrics of localization accuracy and spatial properties were estimated.
Measures of distance from the focus and across methods: Quantitative assessment of source maps spatial characteristics.
1 Spatial dispersion (SD): SD measures the spatial spread of a current source density map and is computed as follows: where Θ is the reference region (the epileptic focus in our case), andĵ i is the amplitude of the current density distribution estimated for the dipolar source i. The amplitude of all the p cortical sources is weighted by their minimum distance from all the dipolar sources belonging to Θ. min j∈Θ D ij À Á provides the minimum distance between the source i and the closest dipolar source of Θ. SD provides a measure in millimeters and a source localization method with either high spatial spread around Θ or high localization error will result in high SD. Further details on the computation of this quantitative measure can be found in Molins, Stufflebeam, Brown, and Hämäläinen (2008) and Pellegrino, Hedrich, et al. (2016). The effect of thresholding on SD was also estimated by thresholding the map J before estimating the resulting SD.
2 Map size: dMSI map size was measured as the number of active cortical vertices. This measure is sensitive to the threshold applied to dMSI maps. When no threshold is applied (denoted as 0% threshold), the entire cortical surface is active, and all cortical vertices have nonzero values. As the cortical surface was tessellated into 8,000 points, Map size is 8,000 when 0% threshold is applied. When the maximum threshold is applied (denoted as 100% threshold), all cortical vertices but the one with highest amplitude are set to zero and the Map size is 1. Map size can therefore range between 1 and 8,000, depending on the threshold and on the spatial properties of the dMSI map. This metric was proposed to investigate the spatial extent of dMSI maps as a function of the threshold.
3 Map_Dmin: At a given threshold, the resulting dMSI map was binarized and the minimum distance between this new binary map accounting for the extent of the generator and the epileptic focus is computed.
This metric is expressed in mm. Therefore, as compared to Dmin, which only relies on the position of the map maximum, Map_Dmin accounts for the spatial extent of the generator and is highly dependent on the threshold applied. If the dMSI map maximum is outside the epileptic focus, higher threshold reduces the map size and increases Map_Dmin; conversely, lower threshold increases map size and reduces Map_Dmin.
Ultimately, the differences in Map_Dmin across dMSI methods depend on the spatial properties of the dMSI maps.  F I G U R E 1 (a) Median Dmin was below 1 cm for all methods. The performance of standardized low-resolution electromagnetic tomography (sLORETA) was slightly but significantly better than the one of dynamic statistical parametric mapping (dSPM). The performance of Ave (median = 0 mm) was significantly better than coherent maximum entropy on the mean (cMEM), dynamic statistical parametric mapping (dSPM), and minimum norm estimate (MNE). (b) Reproducibility of Dmin measured as within-subject interquartile range was not significantly different across methods. Its median was below 1 cm for all methods. All graphs are boxplots depicting median and quartiles. The dash line is set to 1 cm. * denotes p < .05; ** denotes p < .001 sLORETA = 7.557 mm; cMEM = 7.352 mm; dSPM = 9.366 mm,  (Figure 2).

| DISCUSSION AND CONCLUSION
Several studies have demonstrated that dMSI has a very good accuracy when compared to ECD, with good overlap between iEEG activity and the reconstructed sources Kanamori et al., 2013;Tanaka et al., 2009;Tanaka et al., 2018).
We now demonstrate that MNE, dSPM, cMEM, and sLORETA have excellent performance with similar distance between map maximum and the epilepsy focus as a measure of accuracy, as well as similar within-subject reliability. The subject-based analysis further confirmed the overall excellent performance of dMSI techniques, with a median Dmin of 0 mm for all methods suggesting that source localizations of data from the same patients were often providing maps with maximum activity localized inside the epileptic focus.
It is important to notice that in clinical settings, there is no a priori knowledge of the true location and extension of the focus. Different inverse solutions applied to the same data can provide different results. In this study we demonstrated that a simple average of maps from different inverse solutions improves significantly the localization accuracy, as compared to any other technique alone (Figure 1a). The overall effect size of such improvement was however small and in the order of only few millimeters. It is possible that several dMSI approaches localize the maximum of the generator in the vicinity of the focus, with a different degree and direction of error. The average, F I G U R E 2 Inter_dMSI distance. This measure provides an estimate of the spatial concordance between one method and all others and is expressed in mm. Inter_dMSI distance was significantly different across inverse methods, but the median value was consistently around 5 cm. This analysis confirms that there is a remarkable variability in the position of the map maximum across inverse solution techniques and supports the attempt to combine inverse methods to achieve a common map owning higher accuracy. Boxplots depict median and quartiles. The dash line is set to 5 cm therefore, may tone down the localization error and take advantage of F I G U R E 3 Spatial properties of distributed magnetic source imaging (dMSI) methods. (a) Spatial dispersion (SD) for all methods, computed without applying any threshold. SD was significantly lower for coherent maximum entropy on the mean (cMEM) when compared to every other method. The boxplot depicts median and quartiles. * denotes p < .05, ** denotes p < .001. (b) SD expressed as function of the threshold value. As the source maximum was typically found in the vicinity of the focus, increasing the threshold reduced the amount of spurious activity for all methods (lower SD), and especially for Minimum Norm Estimate (MNE), dynamic statistical parametric mapping (dSPM), standardized lowresolution electromagnetic tomography (sLORETA), and Ave. The SD curve of cMEM remained rather "flat" and below other inverse solution techniques, suggesting that cMEM maps (a) result into a very high contrast between the generator and surrounding region, (b) most of the map activity is found within the borders of the focus, and (c) the amount of spurious activity is significantly lower when compared to other techniques.
x axis = threshold ranging between 0 and 100%. y axis = SD expressed in mm. The curves denote average values and the variability is expressed as SE over "studies." (c) Size of the source map expressed as function of the threshold value. For a given threshold, the average map size of cMEM was smaller than for other methods. With increasing thresholds, the size of cMEM maps reaches a plateau already at 10-20% threshold (see also zoomed view between 10 and 30% threshold indicated by the dashed box). The other methods exhibited a more regular decrease in extent of the map with increasing threshold, suggesting mainly the influence of the noise level rather than sensitivity to the underlying spatial extent. In other words, within a large range of threshold values, we observed very little influence on the size of the source map for cMEM. x axis = threshold ranging between 0 and 100%. y axis = size of the generator expressed as number of active vertices. y scale ranges between 1 (0% threshold-the entire cortical surface of 8,000 vertices is considered as active) and 8,000 (100% threshold, only the vertex exhibiting maximum amplitude is active). The curves denote average values and the variability are expressed as SE over "studies." (d) Distance from the epilepsy focus (Map Dmin) as function of the threshold. For a given threshold, the average distances of cMEM and minimum norm estimate (MNE) from the focus were larger when compared to the other methods that were more likely to overestimate the size of the generator, especially for lower thresholds. x axis = threshold ranging between 0 and 100%. y axis = distance between the thresholded dMSI map and the focus expressed in mm. The curves denote average values and the variability is expressed as SE over "studies" proposed an interesting theoretical framework to combine several source localization results through Bayesian model averaging (Trujillo-Barreto, Aubert-Vazquez, & Valdes-Sosa, 2004). Averaging multiple dMSI solutions is certainly not the only possible approach. Future studies will need to address how to handle the validity and reliability of the results produced by different methods applied on the same data and how to integrate this information for the benefit of the clinical assessment.
Our study focused on the spatial properties of the maps generated with different dMSI approaches. One of the best advantages of dMSI is its ability to provide an extended generator. Not all inverse solution techniques are able to recover the spatial extent of the generator: MNE, dSPM and sLORETA may reach a good accuracy in localizing the peak of the source, but typically provide large maps, often biased, especially for deep generators. (Chowdhury et al., 2013;Chowdhury et al., 2016;Ding, 2009;Lin et al., 2006). For instance, Ding has demonstrated that the minimum norm model fails to recover the continuous cortical distribution of extended sources (Ding, 2009).
In agreement with our previous report (Heers et al., 2016) The nonthresholded distributed magnetic source imaging (dMSI) maps for all methods. Cortical surface has been inflated to improve its visualization. The reconstruction of the epileptic focus based on clinical information is depicted in magenta. All the dMSI methods provide a good localization of the epilepsy focus. Coherent maximum entropy on the mean (cMEM) map shows high contrast, with the maximum centered in the right frontal operculum and very little activity spreading outside the epileptic focus. Standardized low-resolution electromagnetic tomography (sLORETA) shows a higher amount of activity outside the epileptic focus and some strong localization in the anterior insula. Dynamic statistical parametric mapping (dSPM) map shows very similar features to sLORETA. Minimum norm estimate (MNE) retrieved a less blurred map when compared to sLORETA and dSPM, with a maximum in the right frontal operculum, but also some activity spread over the right frontal and temporal pole. Ave shows a good localization, with less spatial spreading as compared to sLORETA, dSPM, and MNE (p < 0.001, consistently). Overall, the relationship between threshold and distance from the focus suggests that the map size and its distance from the epilepsy focus strongly depend on the threshold applied and therefore on the underlying noise level. Applying a method that is minimally affected by thresholding avoids an additional step that is user-dependent and that could negatively impact the results. Whereas Ave showed a better performance in localizing the map maximum, its spatial properties were similar to those of MNE, sLORETA, and dSPM. This was because Ave was the result of the average of four maps, three of which (MNE, sLORETA, and dSPM) showing significant degree of distant spurious activity (Figures 6   and 7).
The appropriate threshold of source imaging maps strongly depends on the specific source imaging approach and an unambiguous solution to this problem has not yet been found (Becker et al., 2014;Becker et al., 2017;Chowdhury et al., 2016;Jung et al., 2013;Sohrabpour et al., 2016;Zhu et al., 2014). In the case of the cMEM algorithm, the thresholding issue is largely overcome by the inverse solution itself. cMEM is able to switch off some cortical parcels that do not contribute to the inverse solution, while preserving the ability to create a contrast of current intensities within the remaining active parcels. This allowed us in the past to propose and validate an empirical 30% threshold which allows removing most or almost all the background noise (Heers et al., 2014;Heers et al., 2016;Papadelis et al., 2016;Pellegrino, Hedrich, et al., 2018;von Ellenrieder et al., 2016) ( Figures 4 and 5). These findings are in agreement with recent studies, demonstrating that, among the most common dMSI techniques, cMEM recovers the generator with the highest contrast while allowing a careful estimation of the spatial extent (Chowdhury et al., 2013;Chowdhury et al., 2015;Chowdhury et al., 2016;Hedrich et al., 2017;Heers et al., 2014).
The spatial properties described above explain the relationship between distance from the epilepsy focus and map threshold as observed in this study. Although Figure 3c illustrated a better performance profile for MNE, dSPM, sLORETA, and Ave, this was largely dependent on the size of the underlying map. Methods that, for the same threshold, provide more widespread maps, will more often Our cohort included neocortical epilepsy patients, for whom MEG source localization is more likely to contribute during presurgical evaluation (Genow et al., 2004;Knowlton et al., 1997;Knowlton et al., 2006;Mamelak, Lopez, Akhtari, & Sutherling, 2002;Mu et al., 2014;Ryvlin et al., 2014;Shiraishi et al., 2001;Stefan et al., 2003;Stefan et al., 2011).
Realistic simulations are possibly the best initial approach to quantitatively characterize the properties of inverse solutions and all the methods applied here have been previously carefully evaluated using realistic simulations (Chowdhury et al., 2013;Chowdhury et al., 2016;Grova et al., 2006;Hedrich et al., 2017). Biophysical computational models continue pushing further the level of realism of F I G U R E 6 Illustrative patient with focal seizures originating from the right frontal cortex (Patient 14). Each row corresponds to a study. From left to right, interictal epileptiform discharges (IEDs) average, topography at the peak, unthresholded source localization results obtained with Ave, coherent maximum entropy on the mean (cMEM), minimum norm estimate (MNE), dynamic statistical parametric mapping (dSPM), and standardized low-resolution electromagnetic tomography (sLORETA). Source localization was performed at the peak. The reconstruction of the presumed epileptic focus is depicted in magenta. This example shows that all distributed magnetic source imaging (dMSI) methods provide a good localization of the epilepsy focus. There is good reproducibility, but in this specific case MNE, dSPM, sLORETA perform slightly better than cMEM (see supplementary material for details on the distance from the focus and spatial dispersion (SD)) simulations, proposing advanced EEG/MEG generative models Cosandier-Rimélé, Merlet, Badier, Chauvel, & Wendling, 2008;Garnier, Vidal, & Benali, 2016). Simulations, however, are never truly realistic and validation in a clinical setting is necessary to progressively bring new techniques in clinical practice Grova et al., 2016;Heers et al., 2016).
Typically, accuracy of a source localization method is assessed considering the lobar or sublobar concordance with the presumed focus/epileptogenic zone/seizure onset zone (Agirre-Arrizubieta et al., 2009;Heers et al., 2016;Rampp et al., 2019). Here, we performed a personalized evaluation to further validate clinical utility of the applied methods, proposing quantitative metric going beyond sublobar qualitative judgment.
This study demonstrates that all the dMSI techniques that we tested provided excellent localization performance. The average map Ave provided the best localization accuracy, whereas cMEM exhibited the lowest amount of spurious distant activity and was less sensitive to map thresholding as compared to other methods.

CONFLICT OF INTEREST
The authors declare no potential conflict of interest.

DATA AVAILABILITY STATEMENT
Individual results are available as supplementary material. The original raw data supporting the findings of this study are available upon reasonable request to the corresponding authors.
F I G U R E 7 Illustrative patient with focal seizures originating from the right frontal parasagittal cortex (Patient 18). The first row corresponds to Study 1, the second row corresponds to Study 3 (see supplementary material for further details). From left to right, interictal epileptiform discharges (IEDs) average, topography at the peak, unthresholded source localization results obtained with Ave, coherent maximum entropy on the mean (cMEM), minimum norm estimate (MNE), dynamic statistical parametric mapping (dSPM), and standardized low-resolution electromagnetic tomography (sLORETA). Source localization was performed at the peak. The example shows that distributed magnetic source imaging (dMSI) is sometimes able to reconstruct recover an accurate source also when the signal to noise ratio is not ideal and the topographical distribution is complex (first raw). Conversely, source localization might be inaccurate, at times with the maximum localized in the opposite hemisphere even when the signal to noise ratio is high and the topographical distribution is dipolar (lower row). This often occurs when the source is close to the midline ORCID Giovanni Pellegrino https://orcid.org/0000-0002-1195-1421 Tanguy Hedrich https://orcid.org/0000-0002-6238-5415 Ümit Aydin https://orcid.org/0000-0002-6327-7811