Magnetoencephalography resting state connectivity patterns as indicatives of surgical outcome in epilepsy patients

Objective. Focal epilepsy is a disorder affecting several brain networks; however, epilepsy surgery usually targets a restricted region, the so-called epileptic focus. There is a growing interest in embedding resting state (RS) connectivity analysis into pre-surgical workup. Approach. In this retrospective study, we analyzed Magnetoencephalography (MEG) long-range RS functional connectivity patterns in patients with drug-resistant focal epilepsy. MEG recorded prior to surgery from seven seizure-free (Engel Ia) and five non seizure-free (Engel III or IV) patients were analyzed (minimum 2-years post-surgical follow-up). MEG segments without any detectable epileptic activity were source localized using wavelet-based Maximum Entropy on the Mean method. Amplitude envelope correlation in the theta (4–8 Hz), alpha (8–13 Hz), and beta (13–26 Hz) bands were used for assessing connectivity. Main results. For seizure-free patients, we found an isolated epileptic network characterized by weaker connections between the brain region where interictal epileptic discharges (IED) are generated and the rest of the cortex, when compared to connectivity between the corresponding contralateral homologous region and the rest of the cortex. Contrarily, non seizure-free patients exhibited a widespread RS epileptic network characterized by stronger connectivity between the IED generator and the rest of the cortex, in comparison to the contralateral region and the cortex. Differences between the two seizure outcome groups concerned mainly distant long-range connections and were found in the alpha-band. Significance. Importantly, these connectivity patterns suggest specific mechanisms describing the underlying organization of the epileptic network and were detectable at the individual patient level, supporting the prospect use of MEG connectivity patterns in epilepsy to predict post-surgical seizure outcome.


Introduction
Presurgical evaluation of patients suffering from drug-resistant focal epilepsy aims at identifying the epileptogenic zone, i.e. 'the minimum amount of cortex that must be resected (inactivated or completely disconnected) to produce seizure freedom' [1]. The epileptogenic zone remains a theoretical concept first introduced by Talairach and Bancaud as the 'site of the beginning of the epileptic seizures and of their primary organization' [2]. In practice, there is no specific marker of this zone, and hence the goal of the clinician is to localize the seizure onset zone, i.e. the region of the cortex where seizures are initiated. Despite advancements in localizing the seizure onset zone, seizure freedom rates after epilepsy surgery range between 53%-84% for mesial temporal epilepsy and 36%-76% for neocortical epilepsy [3], indicating the need to identify more accurate predictors of surgical outcome.
The concept of focal epilepsy as a condition affecting a focal region of the brain has been recently challenged, and focal epilepsy could rather be a network disorder [4][5][6]. Several studies reported alterations of whole brain resting state functional connectivity patterns even in absence of any detectable epileptic discharges [7][8][9][10]. These functional connectivity studies, from our group and others, investigated slow functional magnetic resonance imaging (fMRI) hemodynamic fluctuations involving frequencies slower than 0.1 Hz [11][12][13][14][15][16][17][18]. The main findings consisted of network reorganizations characterized by reduced overall efficiency, i.e. a loss of distant connections and increased local connections in the vicinity of the epileptic focus. These reorganizations were found to depend on epilepsy duration and seizure frequency [19,20].
fMRI provides an indirect measurement of neuronal activity by measuring oxygen consumption and blood flow. Electroencephalography (EEG) and Magnetoencephalography (MEG) directly measure neuronal activity at the millisecond time scale offering together a complementary and unique approach to investigate the dynamic of resting state fluctuations. The resting state networks (RSNs) identified using EEG or MEG not only exhibit similar spatial patterns to the ones estimated using resting state fMRI, but EEG and MEG also provide a unique opportunity to explore much richer features in time, because of their excellent temporal resolution (at the millisecond level). Indeed, signals in different frequency bands allow studying different RSNs [21] as well as their dynamics [22]. Several resting state MEG or EEG analysis methods have already been considered in epilepsy, focusing on connectivity measures such as coherence analysis [23], lagged linear coherence [24], phase lag index [25,26], imaginary coherence [27,28], lagged phase synchronization [29] and weighted partial directed coherence [30], among others. Although most of these studies reported decreased long-range connectivity associated with epilepsy, findings were not as uniform as for the fMRI studies and demonstrate several discrepancies.
A promising application of functional connectivity patterns in epilepsy is to predict the probability of seizure freedom following resection of the epileptic focus, as suggested using MEG [27], fMRI [31][32][33] or a combination of fMRI and diffusion weighted imaging [34]. In this study, we calculated MEG orthogonalized amplitude envelope correlations (AEC) in theta (4)(5)(6)(7)(8), alpha (8)(9)(10)(11)(12)(13), and beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26) bands within MEG segments in which no epileptic discharges were visually identified. AEC was selected as the metric for connectivity analysis following the studies of Brookes et al [21] and Hipp et al [35], which demonstrated that AEC maps in alpha or beta band exhibited functional networks very similar to the ones obtained from resting state fMRI from healthy participants. The selection of theta band was based on the study of Hipp et al [35], which found hubs in medial temporal structures for this frequency band. MEG AEC has been widely used for studying RSNs in healthy brains [36,37], but has never been considered to characterize epileptic networks, and only one study applied EEG AEC in epilepsy [38]. Importantly, AEC of MEG signals for connectivity analysis does not just replicate the findings from fMRI studies but it also has the potential to provide additional information since AEC at different frequency bands could reveal different RSN patterns [21]. Therefore, using MEG AECs might help to differentiate fMRI connectivity patterns [21,22,[39][40][41], as for instance if epilepsy influences connectivity in the alpha band differently from the beta band, classical fMRI connectivity analysis would not be able separate these effects and would only provide a combined effect. Selecting alpha and beta frequency bands in absence of any detectable IEDs for our analysis, might not be directly measuring the epileptic network, but it rather evaluates the characteristics of the epileptic network by measuring how it interferes with other 'healthy' brain networks.
It is worth mentioning that several other markers have been suggested as potential predictors of surgical outcome with varying level of success. For instance, responses to single pulse electrical stimulation of the cortex [42], hippocampal atrophy and absence of generalized tonic-clonic seizures in temporal lobe epilepsy surgery [43], contralateral entorhinal cortex atrophy [44], reorganization of structural connectome [45], concordance between resected region and regions generating high frequency oscillations [46] or fMRI responses to interictal epileptic discharges (IEDs) [47] have been suggested as predictor of surgical outcome. Our study aims at assessing the possibility to consider MEG resting state connectivity patterns in this context.
As opposed to fMRI, EEG and MEG require solving an ill-posed inverse problem prior to characterizing brain networks along the cortex, from resting state data. A variety of source analysis methods have been suggested such as minimum norm estimate [48], dynamical imaging of coherent sources [49], and linearly constrained minimum variance beamformer [50]. We developed and carefully validated a nonlinear source localization method entitled the waveletbased Maximum Entropy on the Mean (wMEM) [51]. Benefitting from discrete wavelet transformation, wMEM is specifically designed to localize brain oscillations, as demonstrated in our recent studies aiming at localizing oscillatory patterns at seizure onset [52] or interictal bursts of high frequency oscillations [53]. Moreover, by studying the intrinsic spatial properties of the MEM operator, we recently demonstrated that MEM provides high spatial resolution and low spatial leakage [54] in comparison to other standard source localization approaches.
In the present study, we applied wMEM to localize MEG resting state fluctuations on cortical surface, within MEG segments in which no epileptic discharges were visually identified, and, therefore, considered as 'rest' data. Following the methodology proposed by Brookes et al [21] and Hipp et al [35], we assessed differences in functional connectivity measured by amplitude envelope correlations in theta (4)(5)(6)(7)(8), alpha (8)(9)(10)(11)(12)(13), and beta (13)(14)(15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26) bands for MEG data acquired prior to surgery. We then investigated seed based functional connectivity considering the following two regions of interest (ROIs): an ROI determined by localizing the generator of interictal epileptic discharges (ROI IED ), and a contralateral homologous ROI (ROI Hom ) serving, for comparison purposes, as a reference for each patient. We hypothesize that patients who exhibit an isolated epileptic network prior to surgery, characterized by low long-distance connectivity values (low AEC values) between the ROI IED when compared to connectivity with the ROI Hom , are more likely to become seizure-free after a focal resection than patients who exhibit a widespread epileptic network prior to surgery.

Patients cohort
This study was approved by the Montreal Neurological Institute (MNI) Research Ethics Board and all procedures were conducted in compliance with the Code of Ethics of the World Medical Association (1964 Helsinki declaration and its later amendments). All subjects signed a written informed consent prior to participation in the EEG/MEG acquisition.
We reviewed consecutive drug-resistant epilepsy patients that underwent MEG acquisition between 2008 and 2011. The inclusion criteria were: (1) surgical resection following MEG with at least 24 months follow-up (16 patients), (2) absence of large lesions or malformations or widespread atrophy, (3) focal epilepsy with strictly unilateral IEDs during MEG session (13 patients), (4) no magnetization artifacts during MEG acquisitions, (5) a minimum of 20 epochs of, each two-second-long, MEG intervals without artifacts and detectable IEDs (12 patients). Following these criteria, data from 12 patients (age, 32.7 ± 10.6 years) were considered. Seven patients were seizure-free (Engel I) after surgery and five continued to have seizures (Engel III or IV) (table 1).

MEG data acquisition and pre-processing
MEG scans were performed in supine position with a 275-gradiometer CTF system (MISL, Vancouver, Canada) and 11 out of 12 scans were acquired with simultaneous 54-channels EEG (10-20 system plus additional electrodes according to the 10-10 system; Easy-cap, Herrsching, Germany). Electrocardiogram and electrooculogram were acquired with additional electrodes. EEG/MEG scans of patients were acquired as a research protocol, during presurgical work-up, following our standard procedure already described [55,56]. At least six runs, each 6-minutes-long, were acquired for each patient with a sampling rate of 1200 Hz (400 Hz for one patient). Patients were instructed to relax and keep their eyes closed. The head position was continuously monitored with three localizing coils placed on Nasion, left and right preauricular points. Data was processed to 3rd order synthetic gradiometers correction to reduce the influence of background noise.

MRI data acquisition, analysis, and head modelling
In order to accurately model subject specific anatomy, for every subject a T1 weighted anatomical scan ( Boundary element method (BEM) was used to solve the EEG/MEG forward problem using Open-MEEG [57,58] as implemented in Brainstorm [59]. The BEM head models were constructed from individual T1-weighted MRIs and included brain, skull and scalp surfaces with conductivity values of 0.33, 0.0165, and 0.33 S m −1 , respectively [60]. Cortical reconstruction and volumetric segmentation were performed using Freesurfer image analysis suite [61] (http://surfer.nmr.mgh.harvard.edu/). Source space consisted of N Vert (8002) vertices, downsampled from the original Freesurfer mid cortex surface using Brainstorm software, and source orientations were fixed orthogonal to the surface. MRI-MEG coregistration was performed by fitting head shape points digitized using Polhemus system (Polhemus Inc., Colchester, VT, USA) to the skin surface derived from the individual MRI.

Definition of the ROIs: localization of interictal epileptic discharges
IEDs consist in transient and spontaneous epileptic discharges not associated with any clinical manifestation, likely to be generated in similar brain regions as those involved during epileptic seizures [62,63]. IEDs were marked visually by an experienced neurophysiologist (GP) in reading EEG and MEG recordings. The cortical generator of IEDs was localized using EEG/MEG fusion and a consensus map approach within coherent Maximum Entropy on the Mean (cMEM) framework, as previously described in [55]. Unlike standard EEG/MEG fusion strategies considering only Signal to Noise Ratio (SNR) normalization before inverting the problem, MEM fusion comprises additional steps, exploiting the complementarity of both recordings when defining the prior model for cMEM source reconstruction, notably when estimating an underlying data-driven cortex parcellation and initializing the probability each parcel to be active (see [64] for further details). Such fusion based source reconstruction was applied to every single IED marked. Moreover, as opposed to standard approaches consisting in localizing only one average IED and in order to provide more reliable and consistent source localization results, we proposed to compute an IED consensus map obtained by applying hierarchical clustering of all source maps obtained for every single IED. We previously demonstrated the excellent reliability to combine EEG-MEG fusion and consensus map to accurately identify the generator of IEDs along the cortical surface [55], while carefully considering variability between epileptic discharges. In the present study, these consensus maps were used to define an IED generator region-of-interest (ROI IED ) for each patient, as illustrated in figure 1(a). This ROI IED will then serve as seed for resting state connectivity analysis (see next section). The accuracy of cMEM localization and its sensitivity to estimate the spatial extent of the underlying generators were carefully demonstrated in our previous studies [54][55][56][64][65][66]. Another ROI with the same size as the ROI IED was drawn manually on the contralateral homologous region for every patient (ROI Hom , figure 1(b)). This ROI Hom was considered as a subject specific reference when comparing seedbased functional connectivity patterns, a reference allowing us to avoid spurious differences due to variable vigilance levels, medications, and ROI IED size.

Pre-processing and selection of the epochs and baseline
Cardiac and blinking artefacts were corrected using Signal Space Projection (SSP) as implemented in Brainstorm software [59]. Blink events were marked manually on EEG while automatic detection in Brainstorm was used to mark R-peaks from Electrocardiogram (ECG) channel. To isolate the artefacts from each other, R-peak events occurring within 1 s to blinks or IEDs were not included in the SSP process. After visual investigation, components corresponding to cardiac and the blink artefacts were selected separately for correction. In most cases, this corresponded to selecting the first SSP component for each artefact type. Patient specific MEG resting state functional connectivity was calculated from 20 epochs, each twosecond-long, exhibiting no detectable IEDs or artefacts. The rationale for selecting segments without epileptic discharges was to avoid connectivity analysis to be driven by excessive changes in connectivity occurring during IEDs. These 20 selected segments were later concatenated and filtered in corresponding band of interest: theta (4-8 Hz) or alpha (8)(9)(10)(11)(12)(13) or beta (13-26 Hz). We assessed that concatenation of time-series segments followed by filtering did not Figure 1. Analysis pipeline to calculate AEC laterality maps: For each seed vertex within the ROIIED (with magenta borders) (a) or ROIHom (bold black borders) (b) amplitude envelope correlations between the seed and all other vertices that were not inside the ROIs were calculated resulting in one connectivity map for each seed point. (c) Average connectivity map obtained for seeds in ROIHom was subtracted from the average connectivity maps for seeds in ROIIED resulting in an AEC difference map (d) Paired t-test between the connectivity maps was calculated for the seeds in ROIIED and ROIHom, and the vertices showing statistically significant differences (p < 0.01 after Bonferroni correction for the number of vertices on the cortex (NVert)) were thresholded to create a mask. Regions found to be statistically significant vertices were assigned to one (in red color in the figure), while the non-significant vertices were assigned zero and masked out. (e) The AEC difference map obtained in (c) was masked with (d) resulting in an AEC laterality map exhibiting only the statistically significant differences. The vector containing the AEC laterality map values of all vertices shown in (e) was denoted as LM, and LM corresponded to the scalar value denoting the corresponding empirical mean of the AEC laterality map values estimated over all vertices.
introduce any artificial off-set between segments that could lead to spurious power increase in certain frequency bands during the subsequent spectral analysis. We calculated the power spectrum of each twosecond-long segments and averaged the spectral results. This power spectrum was compared with the power spectrum calculated from the concatenated signal. Comparing these two power spectrums we did not observe any major effects of the concatenation process in the frequency bands of interest.
Before applying wMEM source localization in the band of interest, an appropriate noise covariance model had to be defined [51]. To do so, we calculated Hilbert transform of the MEG sensor data in the frequency band of interest and estimated the average absolute amplitude over all channels at each time instant. From this amplitude envelope signal in the sensor space, a two-second-long interval exhibiting the lowest absolute amplitude in the band of interest was used as baseline to estimate a diagonal noise covariance matrix. We ensured that this baseline segment did not overlap with any IEDs, artifacts or with the epochs selected for resting state analysis.

Source localization of resting state MEG data
Resting state MEG data were localized using the wavelet-based MEM (wMEM) [51]. This approach uses discrete wavelet transformation (Daubechies wavelets) in the sensor space to calculate a discrete time-frequency representation of the measurements. Discrete wavelet representation allows selecting and localizing only the time-frequency components of interest. wMEM inverse operator then allows estimating the corresponding time-frequency representation of the signal in the source space, i.e. along the cortical surface. Afterwards, the time courses of the sources were calculated from these estimated discrete wavelet coefficients, using an inverse wavelet transform [51]. wMEM does not only inherit the high spatial resolution and low spatial leakage properties of MEM [54], but also benefits from the efficiency of wavelet transformation to propose a sparse representation of the entire temporal dynamics of oscillations, using only a few specific time-scale components [67].  [21,35], we applied Hilbert transform to localized resting state source signals along the cortical surface. The purpose of Hilbert transform was to capture amplitude envelope oscillations within the frequency bands of interest (theta, alpha, and beta). Using these envelopes, we finally assessed MEG amplitude envelope correlations (AEC) between different brain regions as a measure of resting state functional connectivity. Hilbert transform was first calculated separately for consecutive frequency bands of 0.5 Hz widths. Resulting envelopes were then averaged to characterize the full frequency band of interest. For instance, when considering analysis in the beta band, Hilbert transform was first calculated within the following frequency bands 13-13.5, 13.5-14.0, …, 25.5-26.0 Hz separately and results were then averaged to calculate orthogonalized beta band amplitude correlations. This process was done to avoid calculating the Hilbert transform over wide frequency band signals.

Orthogonalized amplitude envelope correlations and AEC laterality maps
To address leakage due to volume conduction in MEG, we chose to remove source signals exhibiting zero phase connectivity with the seed, since most of them might be inherent to volume conduction issues and resolution of the inverse problem, therefore resulting in biased connectivity measures. To do so, we used the orthogonalization procedure proposed in Hipp et al [35]. After orthogonalization of Hilbert transforms, we produced seed based AEC connectivity maps between the seeds, i.e. ROI IED , and ROI Hom respectively, and the rest of the cortex.
In order to consider every subject as his/her own reference, we then generated an AEC laterality map that compared the connectivity differences between the seed-based AEC connectivity maps with ROI IED , and ROI Hom . The process to calculate this AEC laterality map is illustrated in figure 1. For each seed vertex within the ROI IED ( figure 1(a)), and the ROI Hom ( figure 1(b)), we computed orthogonalized AEC between the time courses of this seed vertex and all other vertices of the cortical surface that were outside these ROIs. Therefore, one connectivity map was estimated for every seed vertex belonging to each ROI. Afterwards, the connectivity maps obtained for all the seeds located in ROI IED (and ROI Hom respectively) were averaged. The averaged connectivity map for ROI Hom was subtracted from the averaged connectivity map estimated for ROI IED , resulting in an AEC difference map (figure 1(c)). Within this difference map, positive values corresponded to vertices exhibiting larger correlations with signals from the ROI IED when compared to signals from the ROI Hom . Negative values corresponded to vertices exhibiting smaller correlations with signals from the ROI IED when compared to signals from the ROI Hom . Considering all connectivity maps obtained from every seed in both ROIs, we performed two tailed t-tests for each vertex located outside the two ROIs. A mask of all vertices showing a t-value significantly different from zero was estimated, considering p-value < 0.01 Bonferroni corrected (corrected p-value threshold = 0.01/N Vert = 1.25 10 -6 ) (figure 1(d)). This Bonferroni corrected t-test procedure is a relatively conservative approach that preserves the most pronounced differences, especially considering signals at each vertex were not completely independent from each other due to volume conduction and inverse solution. Prior to t-test, Fisher transformation was also applied to all AEC values in order to reduce the skewness of the distribution [68]. On average, each ROI was composed of 57 ± 29 vertices, providing the corresponding range of degrees of freedom of the t-test. Finally, only significant differences were kept by multiplying the difference map (figure 1(c)) with this binary mask (figure 1(d)), resulting in an AEC laterality map reporting only statistically significant differences ( figure 1(e)). The vector containing the AEC laterality map values from all vertices reported in figure 1(e) is denoted as LM, and the scalar LM denotes the corresponding average of the AEC laterality values estimated over all vertices.

Influence of long-range connections in AEC connectivity
To assess and visualize the importance of long-range resting state connections for characterizing the epilepsy network, we performed a quantitative analysis of AEC patterns as a function of the distance to each ROI. We grouped the vertices according to their Euclidian distances from the ROI and calculated the average AEC value within each group. Inter-subject anatomical variability was addressed by using normalized distances, indicated as a percentage of the distance between the seed region and the most distant vertex on the cortex, for each individual. For example, we first calculated the mean AEC value over the vertices that were closer than 10% of the maximum distance to the seed ROI, then over the vertices exhibiting a distance ranging from 10% to 20% distance from the seed ROI and so on. These AEC values were then presented as bar plots in results section.
Finally, we also proposed a Laterality difference index (LDI) to quantify the change in laterality maps with respect to distance. We used a relative measure in percentage that normalizes the AEC laterality map differences, to provide an index that was not influenced by changes to the overall connectivity strengths for different patients and distances. The LDI was calculated as follows: where d is the center of the distance interval in which LDI is calculated (% of the maximum distance). ϕ d is the set of vertices that were located within a distance ranging between d − 15% to d + 15% to the corresponding seed ROI. LM ϕ d is the average of AEC laterality map values for the vertices in the set ϕ d , and |LM ϕ d | is the average of the absolute values of AEC laterality map values for the vertices in the set ϕ d .

Specificity analysis using healthy control data
We proposed a non-parametric testing scheme with random drawings to assess if differences in resting state patterns of seizure-free versus non seizure-free patients could be explained solely by the MEG sensitivity to the chosen anatomical locations of interest or the sizes of the ROIs. For this purpose, ROI locations from the 12 epilepsy patients and source localized from MEG resting state recordings of five healthy controls (age, 30.0 ± 7.6 years) were projected to the Montreal Neurological Institute ICBM152 template with nonlinear surface-based registration using Freesurfer spheres, as described in [69] (Brainstorm implementation). Following the same procedures as explained for the epilepsy patients, we calculated 60 AEC laterality maps from the healthy control MEG recordings, using anatomical ROIs extracted from patients' data. These maps were classified as 'seizurefree group' (35 AEC laterality maps, 35 = 7 ROIs × 5 subjects) when the corresponding anatomical ROI locations were obtained from seizure-free patients, and as 'non seizure-free group' (25 AEC laterality maps, 25 = 5 ROIs × 5 subjects) when the corresponding ROI locations were obtained from non seizure-free patients. Then, 12 out of 60 AEC laterality maps were randomly drawn, resulting in seven AEC laterality maps LM from the 'seizure-free group' and five from the 'non seizure-free group' . Drawn was repeated 20 000 times without replacement. Our null hypothesis was that seven mean AEC laterality map values (LM) from non seizure-free group are not larger than the five LM from the seizure-free group. We tested this hypothesis for each draw by calculating one tailed Wilcoxon rank-sum test. Finally, the histogram of these 20 000 Wilcoxon rank-sum test statistics (null hypothesis distribution) was used to assess the significance of the rank-sum value obtained from the patients' data. Detailed description of this procedure is provided in Supplementary Material.

Data availability
The data that support the findings of this study are available upon request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.

Results
In this section, we only report results obtained for fluctuations in the alpha band (8-13 Hz) considering it is the most prominent brain rhythm and it provided the most significant information related to seizure outcome. Findings for the theta and beta bands are provided in Supplementary Material.

Differences in MEG resting state connectivity patterns between seizure-free and non seizure-free patients
A complete analysis performed on MEG data obtained for a seizure-free patient (P5) with frontal lobe epilepsy is illustrated in figure 2. After marking eight epileptic spikes, source reconstruction of IEDs using EEG/MEG fusion with cMEM localized the ROI IED on the right orbito-frontal region (figures 2(a) and (b)). AEC of MEG resting state data exhibited a considerably more isolated resting state network for the IED generator ROI IED (figure 2(c)) when compared to a more extended (sometimes bilateral) network obtained from the homologous contralateral region ROI Hom (figure 2(d)) (non thresholded maps). Consequently, the corresponding AEC laterality map LM represented in figure 2(e) was dominated by negative values, indicating significantly lower AEC connectivity between ROI IED and bilateral parietal, occipital and posterior temporal regions when compared to ROI Hom . MRI investigation suggested a left parasagittal Focal cortical dysplasia (FCD), which was confirmed as type IIB after surgery. Patient was operated from the right orbito-frontal focus based on converging evidence from intracranial depth EEG recordings. Patient is seizure-free (Engel Ia) since the surgery (8 years). Findings for a non seizure-free patient (P9) with frontal lobe epilepsy are shown in figure 3. Nine left fronto-central spikes were marked and IED source reconstruction localized a left frontal region as IED generator at the peak of the spikes (figures 3(a) and (b)). AEC maps (non thresholded) for ROI IED (figure 3(c)) and ROI Hom (figure 3(d)) indicated a largely extended network connected with the IED generator, when compared to a more confined network associated with the ROI Hom . The corresponding laterality map LM (figure 3(e)) showed significantly higher AEC between ROI IED and bilateral orbitofrontal and some parietal regions when compared to ROI Hom . MRI investigation found a left parasagittal FCD, which was confirmed as type IIB after surgery. Patient went through resection of left parasagittal frontal region based on converging evidence from seizure semiology, scalp EEG and stereo EEG, however, continues to have seizures (Engel IV).
Findings for a non seizure-free patient (P11) with temporal lobe epilepsy are shown in figure 4. 51 left temporal spikes were marked and IED source reconstruction localized a left temporal superior IED generator at the peak of the spikes (figures 4(a) and (b)). AEC maps (non thresholded) for ROI IED (figure 4(c)) and ROI Hom (figure 4(d)) indicated a largely extended network connected with the IED generator, when compared to a more confined network associated with the ROI Hom . The corresponding laterality map LM (figure 4(e)) showed significantly higher AEC between ROI IED and bilateral occipital and frontal poles when compared to ROI Hom . MRI investigation showed left hippocampal rotation. Stereo-EEG showed neocortical and mesio-temporal widespread activity, therefore the lateral ROI we used in this study was considered as part of the epileptic network. Patient was operated (Amygdalohippocampectomy) based on converging evidence from seizure semiology, scalp EEG and stereo EEG, however, continues having seizures after surgery (Engel III).

Average AEC map values and LM for seizure-free and non seizure-free patients
We first computed average AEC map values for ROI IED and ROI Hom , and average AEC laterality map values (LM) for each of the 12 patients (averages calculated over all vertices of the cortex excluding those within the ROIs). No trend was observed between seizure free and non seizure free patients when considering average AEC map values, whereas the LM values were found significantly different for the two groups with negative LM for 6 out of 7 seizure-freepatients and positive LM for 4 out of 5 non seizurefree patients (see supplementary figure 2 is available online at (stacks.iop.org/JNE/17/035007/mmedia)). A negative LM indicates that, considering all the vertices outside the ROIs on the whole cortex, the average connectivity with the IED generator was significantly lower in comparison to the connectivity with the contralateral homologous area. A positive LM indicates that, considering all the vertices outside the ROIs on the whole cortex, the average connectivity with the IED generator was significantly higher in comparison to the connectivity with the contralateral homologous area.
Investigating clinical information reported in table 1, all three patients with mesial temporal sclerosis indicators in MRI were seizure-free after surgery. Other than that, groups were similar in terms of MRI findings and there were patients with FCD in both surgical outcome groups. However, for seizure free patients, histology reports one FCD type IIA, two FCDs type IIB and one FCD type IIIA. On the other hand, for non-seizure free patients, among four patients being diagnosed with FCD, three were of type IIA and one of type IIB. It is worth mentioning that the only seizure-free patient with FCD type IIA was P4 and he was also the only one showing a positive LM in the seizure-free group.

Influence of long-range connections in AEC connectivity
Histograms of AEC values as a function of the distance from the ROI are presented in figure 5, providing a more detailed view on the variation of AEC values as a function of the distance to the ROI. Investigating AEC values measured for vertices located close to the ROIs, we did not observe any clear difference between whole brain AEC values for seeds in ROI IED and ROI Hom . Even if orthogonalization was performed, these vertices close to ROIs correspond to regions that were under considerable influence of volume conduction effects. However, when considering only vertices located at a distance larger than 50% of the maximum distance, we found overall smaller AEC values associated with ROI IED in comparison to ROI Hom for six out of seven seizure-free patients. On the other hand, we found that AEC values associated with ROI IED for distant vertices were larger than AEC values associated with ROI Hom for all five non seizure-free patients. Therefore, as presented in figure 5, the differences found in AEC connectivity patterns between seizure free and non seizure-free patients were mainly due to differences in long-range connections. To further characterize and test the significance of the differences in long-range connections, we assessed the evolution of normalized laterality difference indices (LDI) as a function of the distance from the seed, as presented in figure 6. The LDI values were calculated for vertices that were within a certain distance. The horizontal axis represents the distance range (as percentage of the most distant vertex) for which the LDI was calculated. LDI values for the seizure-free group were significantly smaller than LDI values for the non seizure-free group, for all distances starting within the 40%-70% distance range and above (Wilcoxon rank-sum test: p = 0.015 for 40%-70% distance range, p = 0.009 for 50%-80% and 60%-90%, p = 0.024 for 70%-100%). Differences between seizure-free and non seizure-free groups were therefore larger when considering long range distant connections only, instead of all connections. LM values, considering all vertices, correctly identified 10 out of 12 patients with respect to their seizure outcome (supplementary figure 2), whereas LDI values for long range distances (>40%) correctly classified 11 out of 12 patients ( figure 6). These results are summarized in figure 7, where AEC laterality maps and LDI values estimated for 60%-90% distance range are provided for every patient, illustrating clear differences between seizure-free and non seizure-free patients. Figure 7 also illustrates the overall agreement between the laterality map and the long-range LDI values, suggesting the possibility to identify a Figure 5. The AEC values with respect to the distance from the ROI: Bar plots showing the mean AECs as a function of the distance from the seed ROI for seven seizure-free and five non seizure-free patients. Thick magenta bars were calculated from the AEC maps for the seeds in ROIIED and the distances to the ROIIED were shown. Similarly, the thin black bars were calculated from the correlation maps for the seeds in ROIHom and the distances to the ROIHom were plotted. The distances were given as a percentage of the maximum distance between the seed region and the furthest away vertex on the cortex for each individual. Figure 6. LDI values as a function of the distance to the seeds: The distance between the seed region and the furthest away vertex on the cortex for each individual was used to normalize the distances given in the horizontal axis. The LDI values were calculated for vertices that were within a certain distance as explained in methods section. The horizontal axis represents the distance range (as percentage of the most distant vertex) over which LDI was calculated. Therefore, the first distance interval corresponded to the LDI value for vertices within a distance range of 0%-30%, while the second interval indicated a distance range from 10% to 40%, and so on. Magenta triangles and black filled circles illustrate the non seizure-free and seizure-free patients, respectively. Seizure-free patient P4 exhibiting a different behavior than others in seizure-free group is shown in gray. One tailed Wilcoxon rank-sum test was used to check for our hypothesis of non seizure-free patients exhibiting larger LDI values than seizure-free patients for different intervals and statistically significant results were indicated with blue stars ( * p < 0.05, * * p < 0.01).
biomarker at the individual level. We provided the figure showing LDI as a function of the distance from the seed for theta band in supplementary figure 4. Analysis of theta band oscillations exhibited overall similar patterns to alpha band oscillation, i.e. seizure free patients presenting smaller LDI values when compared to non seizure-free patients (Wilcoxon rank-sum test: p < 0.05 for 20%-50% distance range, p < 0.01 for 30%-60%). However, this theta-band analysis presented overall less classification accuracy Figure 7. Laterality maps and LDI values for all patients: Patients 1-7 were seizure-free after and 8-12 were non-seizure-free. A positive value (red) refers to significantly larger correlation with the ROIIED, whereas a negative value (blue) refers to significantly larger AEC with ROIHom. Laterality difference index (LDI) given shows the relative connectivity difference, normalized to the total connectivity, for long range connections (60%-90%) of the maximum distance from the ROI).
considering LDI value for 60%-90% distance range, since P1, P2, and P8 were misclassified (P1 and P2 being seizure free patients showing positive LDI values and P8 a non seizure-free patient presenting a negative LDI value), whereas in this band P4 was correctly classified. These preliminary results on a too small database could not allow any further conclusions for this analysis, except that analysis in different bands should be carefully considered in future studies involving larger populations.

Specificity analysis using healthy control data
We applied a non-parametric hypothesis test to estimate the level of chance to observe laterality map differences similar to the epilepsy patients, when considering resting state MEG data from healthy controls, using the anatomical ROIs defined in our patient group. The null distribution of Wilcoxon ranksum test statistic was computed using 20 000 random drawings from 60 laterality maps, calculated from MEG recordings of healthy controls (see corresponding section in methods and the supplementary material). From this empirical null distribution, a Wilcoxon ranksum test statistic value equal to or smaller than the one calculated for our patient group could be found only in 1.5% of the cases (p < 0.015). We can conclude that the difference obtained in LM in seizure free versus non-seizure free patients was not driven solely by the anatomical locations of the generators or the size of ROIs in the two groups of patients.

Discussion
We investigated pre-operative MEG resting state functional connectivity patterns in patients with focal refractory epilepsy who underwent epilepsy surgery. Only patients with a minimum of two years follow-up were selected (average follow up = 65 months, range 32-97 months). Our main objective was to identify MEG resting state patterns that would predict surgical outcome from data acquired prior to surgery. For this purpose, MEG intervals not exhibiting any detectable epileptic activity were selected and MEG resting state activity was localized on individual cortex surface using wavelet-based Maximum Entropy on the Mean (wMEM). We calculated orthogonalized amplitude envelope correlations (AEC) in theta (4-8 Hz), alpha (8)(9)(10)(11)(12)(13), and beta (13-26 Hz) bands and used a seed-based approach for connectivity analysis. AEC was selected as the metric for connectivity analysis following the studies of Brookes et al [21] and Hipp et al [35], which demonstrated that AEC maps in alpha or beta band exhibited functional networks very similar to the ones obtained from resting state fMRI from healthy participants. Importantly, changes in resting state fMRI connectivity have been reported in several studies in epilepsy [13,[16][17][18], building the bases for our rationale of selecting AEC based metric for studying MEG resting state connectivity. Since our analysis was done in absence of any detectable IEDs, we did not directly measure the epileptic network, but we rather evaluated its characteristics by measuring how it interfered with other 'healthy' brain networks. AEC Laterality map and laterality difference index as a function of the distance to the IED generator were used to compare connectivity of the IED generator and the rest of the cortex against the connectivity between the contralateral homologous region and the rest of the cortex.
We found distinct connectivity patterns in the alpha band that allows differentiating between seizure-free and non seizure-free patients, using only resting state fluctuations in absence of any detectable epileptic discharges. We did not observe any significant difference in beta band fluctuations.
For alpha band fluctuations, seizure-free patients exhibited weaker connections between the IED generator and the rest of the cortex when compared to contralateral homologous region, suggesting a more isolated resting state network associated with the IED generator. On the other hand, non seizurefree patients were exhibiting stronger connections between the IED generator and the rest of the cortex, in comparison to the homologous region, suggesting a more widespread network. Detailed analyses revealed that these differences in resting state network structure were mainly present in long-range distant connections to the seed. Most importantly, these distinct connectivity patterns were not only observed at the group level, but they were also present and reliable at the individual level for 11 out of 12 patients studied. Connections to the contralateral side might play a big role in the long-distance connections we observed. In order to test this, we calculated the Laterality map (LM) values only for the vertices contralateral to the IED generator and plotted these values as boxplot distributions in supplementary figure 5. We found that negative median LM values contralateral to the IED generator, for six out of seven seizure free-patients, and zero or positive median LM values contralateral to the IED generator, for four out of five non seizurefree patients. Based on these results we could state that contralateral connections seem to play an important role on seizure-freedom and most but not all of the differences we observed between seizure-free and non seizure-free patients were due to these contralateral connections.
Unlike the predominantly negative AEC laterality maps observed for seizure-free patient, we found predominantly positive AEC laterality maps in non seizure-free patients, indicating a higher connectivity between the IED generator and the rest of the cortex in comparison to the connectivity between the contralateral homologous non-epileptic side and the rest of the cortex. This specific pattern could, therefore, reflect an underlying widespread epileptic network, involving long range distant regions connected to the main focus, as suggested in EEG-fMRI investigations [70]. Another possibility is that those patients exhibit less or insufficient compensatory mechanisms of the other 'healthy' brain networks, the contralateral ones notably, as opposed to the main trends observed in the seizure free group. Although not studied in the context of surgical outcome, increased contralateral fMRI connectivity and decreased connectivity within the epileptogenic networks were reported for patients with mesial temporal lobe epilepsy in comparison to healthy controls [71]. The increased compensatory contralateral activity hypothesis also agrees with the brain network model proposed by Fornito et al [72], suggesting anatomical and functional connectivity increase in regions distant to the pathological area.
Decreased distant connectivity in epilepsy, resulting in the evolution from a small world network configuration, which is characteristic of the healthy condition, to a more regularized network isolating the epileptic focus from the rest of the brain, was suggested by several functional and structural connectivity studies [15,25,73]. Using an fMRI data-driven sparse modeling approach, our group further demonstrated that such network reorganizations were characterized by significant disruptions of connector hubs within the mesial temporal lobe and the default mode network [13]. In a case report, Jackson et al [74] showed that the global fMRI functional network structure returned back to normal, i.e. restoring the small world network configuration, after successful surgery, even if the surgery was only involving a very small lesion.
A more lateralized fMRI functional connectivity has been reported for seizure-free in comparison to non seizure-free patients [32]. On the other hand, He and colleagues [31], also using resting state fMRI connectivity, showed an increase in number of connections between thalami and contralateral hemisphere in non-seizure-free patients with temporal lobe epilepsy, a pattern not found in seizure-free patients or healthy controls. This finding is in agreement with the more widespread epileptic networks found in cases of surgery failure in our study. Higher functional connectivity between the propagation zone and brain regions not-involved during the seizure was shown to be associated with worse seizure outcome using intracranial stereo-EEG as well [75].
However, using imaginary part of the coherence on MEG resting state signals, Englot and colleagues [27] showed higher connectivity between the resected region and the whole brain (in comparison to the homologous region and the whole brain) as an indicator of good surgical outcome for epilepsy patients. Whereas this seems conflicting to our findings, it is important to note that there is nothing that could prevent an increase in alpha-band AEC to be concomitant with a decrease in alpha band coherence. Coherence is indeed measuring spectral coupling (with constant phase shift), independently from the respective amplitude power. Therefore, these two metrics might refer to two different underlying mechanisms. Still, the neuronal origins of such resting state connectivity links are poorly understood and should be further investigated along with careful comparison between connectivity metrics, which falls out of the scope of this study [76]. Unlike fMRI, in which most studies use correlation as the connectivity metric, MEG studies provide access to a wide range of connectivity metrics, each focusing on different aspect of the signals, therefore studying different underlying mechanisms, to quote a few: coherence analysis [23], lagged linear coherence [24], phase lag index [25,26], imaginary coherence [27,28], lagged phase synchronization [29] and weighted partial directed coherence [30]. Therefore, seemingly high variability of the results reported for MEG connectivity studies could be explained by the high number of connectivity metrics, each of them focussing on different aspects of resting state dynamics. Furthermore, functional connectivity may be dynamically modulated in the epileptic zone [77]: increased functional connectivity during the IED, whereas during quiescent periods perhaps the functional connectivity is modulated downwards.
The ROIs considered to build AEC maps were obtained from MEG/EEG source localization of IEDs [55]. Therefore, for temporal lobe epilepsy patients, those ROIs were temporal lateral and not including the mesial aspects of the temporal lobe. However, it is important to mention that the resulting laterality maps for all six TLE cases were in agreement with our hypotheses. In many mesial TLE cases the signals at the peak of EEG/MEG IEDs are generated by lateral regions, where the activity propagates from the mesial structures [78][79][80]. Since in these cases the lateral and the mesial structures are part of the same epileptic network as suggested by the propagation of the spikes, we expect to capture the characteristics of the epileptic network even when using seed ROIs in temporal lateral structures.
The patient cohort in this study also consisted of seven FCD type II patients. It has been suggested that the balloon cells in FCD type IIB might disrupt connectivity to surrounding networks, which might result in a more isolated epileptic network and negative AEC laterality maps, as we observed for seizurefree patients in our study [81,82]. It is worth mentioning that the only seizure-free patient that showed an 'unexpected' positive laterality map was P4, who actually was diagnosed with an FCD of type IIA. All other three patients diagnosed with FCD type IIA were in the non seizure-free group. The absence of balloon cells in FCD type IIA might explain the positive laterality maps observed for these patients. However, a study with a larger cohort would be necessary to study this distinction in connectivity resting state patterns between these two FCD subtypes.
We focused on alpha and beta band AEC maps based on the study of Brookes et al [21] showing that reconstructing MEG resting state networks in these frequency bands provided network maps very similar to the ones reconstructed with resting state fMRI in healthy controls. Furthermore, highest correlation between global electrophysiological and hemodynamic brain networks was found to be within the alpha and beta bands, partly due to higher signal strength at these frequency bands in comparison to gamma band [35]. Similarly, using high density EEG resting state data in healthy subjects, Liu et al [83] reported that, when considering alpha (8-13 Hz) and beta (13-30 Hz) bands, they were able to reconstruct 13 out of 14 fMRI resting state networks, using independent component analysis on power envelopes. The default mode network (DMN) was fully reconstructed when considering alpha band oscillations, but only partially within the beta band. In another study, also involving healthy subjects, the DMN was reconstructed only in MEG alpha band, whereas other physiological resting state networks, such as sensorimotor or visual networks, were reconstructed from beta band oscillations [21]. On the other hand, our results are showing distinct connectivity differences between seizure-free and non seizure-free patients mainly in the alpha band, suggesting that epilepsy specific reorganizations could be associated mainly to distant connectivity changes within the DMN. For instance, patients with mesial temporal lobe epilepsy showed decreased fMRI connectivity from amygdala and hippocampus to DMN [18], and demonstrated significant disruptions of connector hubs within the mesial temporal lobe and DMN [13]. In accordance with these fMRI studies, DeSalvo et al [84] reported decreased long range structural connectivity as well as increased local connectivity within the DMN using diffusion weighted imaging. A study from Van Dellen et al [85] showed a positive correlation between lower alpha band connectivity and seizure frequency, which suggests that our findings specific to alpha band might be partly related to seizure frequency. Indeed, high seizure frequency prior to surgery has been associated with poor seizure outcome [86]. Considering the studies of Van Dellen et al [85] and Foldvary et al [86] together a link between alpha band connectivity and seizure outcome seems possible, although testing this hypothesis would deserve further investigations that were out of the scope of this study. It is also important to note that we are studying MEG oscillation within the 8-13 Hz, therefore the effects of mu-rhythms could also be a reason for observed connectivity in this band [87].
A major problem with functional connectivity studies investigating differences between epilepsy patients and healthy controls is to control for vigilance states and changes due to anti-epileptic drugs (AEDs) [88]. Investigating the effects of AEDs in fMRI connectivity, Hermans et al [89] reported increased connectivity after AED withdrawal questioning the cause and effect relationship for studies reporting decreased connectivity for epilepsy patients. As opposed to group comparison studies, which are prone to limitations due to differences in vigilance and medication, here we propose a metric in which the contralateral hemisphere serves as reference for the analysis to cope with these limitations. Therefore, by using laterality map and laterality difference index values, we aimed at minimizing the effects of several confounds that could influence MEG background activity. We believe that choosing such an approach was important to provide a marker of the surgical outcome at the individual level.
We acknowledge several methodological limitations, which warrants caution to interpret our findings and will deserve further validations. Even if we reported specific connectivity patterns in agreement with the postsurgical outcome for 11 of 12 patients, the small patient cohort limits generalization of our results. Furthermore, the studied patient population is heterogenous, involving six frontal lobe and six temporal lobe epilepsy patients. However, it is important that we found similar connectivity patterns that are related to seizure freedom even within this heterogenous group. It also requires caution when using our method for patients with bilateral epileptic activity, due to the comparison step of the IED generator and the contralateral homologous region. A high-density EEG study reported abnormal functional connectivity to be associated with cognitive impairments in epilepsy patients [30]. In this direction, the neuropsychological report for the seizure free patient P4 was indicating bilateral frontal deficits, which might be another reason for the very distinct connectivity pattern for this patient. Further investigations using other modalities (fMRI connectivity) or post-surgical recordings would be needed to clarify this issue. It is also important to note that we used the IED generator as the region of interest and not the resected region. One main reason for this choice was to ensure the analyses we suggested could also be employed in future prospective studies, before the surgery. It is important to mention that in this study, the localized IED generator was concordant with the resected regions for all extratemporal patients and it was on the main propagation pathway for all patients with selective amygdala hippocampectomy, in line with the reliability of MEM source localization of IEDs demonstrated in our previous studies [55,56]. It should be noted that orthogonalized AEC procedure removes zero lag signals to avoid spurious connectivity due to volume conduction (spatial leakage). It has been shown that real zero lag signals exist [90]; thus there is a risk in using orthogonalization. However, considering the dominant effects of spurious zero lag connections [36], we have decided to use orthogonalized signals in this study. Finally, we employed individual three-layer boundary element head models for solving the forward problem. Cho et al [91] showed that head models and thus the forward problem could affect connectivity analysis. Along with other studies showing the influence of forward problem in source localization [92], it would be interesting to employ more advanced head models [93] especially for EEG or simultaneous EEG/MEG source connectivity analysis.
It is important to note that in our proposed study, as in any functional connectivity analysis, our results are likely to depend on specific decisions related to the use of a specific source localization method and the choice of a specific functional connectivity metric. Such inherent variability is mainly due to the ill-posed nature of the inverse MEG problem. Each source localization algorithm imposes different constraints to obtain a unique solution and the results will be in line with these assumptions and constrains. In Hedrich et al [54], we compared the resolution matrices to characterize the intrinsic spatial resolution properties of Minimum Norm Estimate (MNE), dynamic Statistical Parametric Mapping (dSPM), standardized Low-Resolution Electromagnetic Tomography (sLORETA) and Maximum Entropy on the Mean (MEM) methods. Whereas localization error was similar for all methods, we found an overall better performance of MEM in terms of spatial resolution (spatial dispersion) and more importantly in terms of spatial leakage, which are important features for connectivity analyses. However, it is important to note that these intrinsic properties of the resolution matrix were evaluated in noise free conditions, therefore additional careful evaluations on real data would still be needed. Another important aspect is the metric chosen for connectivity analysis. In this direction, Colclough et al [36] tested 12 methods to estimate resting state network properties in MEG. Their results showed that AEC and partial correlation measures provided higher test-retest reliability than other metrics such as phase lag index and imaginary part of the coherence, which performed poorly. Whereas careful comparison between source localization methods and connectivity metrics was out of the scope of our proposed study, we assessed test-retest reliability of our findings using a bootstrap resampling strategy. Our results are suggesting an overall good test/retest reliability of our proposed methods considering wMEM for source localization and AEC for connectivity metric (Please see 'Test-Retest Reliability Analysis' section in supplementary material).
In addition to increasing the sample size, critical for further clinical validation, another important step to further improve the confidence of our results could be by considering multimodal studies. Morgan and colleagues [34] reported for TLE patients that considering both structural (diffusion MRI) and functional (fMRI) connectivity networks could provide a better prediction on the seizure outcome in comparison to relying on a single modality. However, considering that most of the information was brought by fMRI in their study, a more promising candidate for multimodal integration could be combining two modalities used for assessing functional connectivity, as for instance MEG and fMRI. Our group has demonstrated that the connector hubness calculated from resting state fMRI could bring important information on the structure of functional networks [13]. By combining fMRI, an indirect measure of neuronal activity through hemodynamic coupling, and MEG, a more direct measure of neuronal activity, as well as different functional connectivity metrics (AEC versus hubness) should provide converging evidence of our findings, as suggested by our preliminary data [94].
We demonstrated that specific resting state connectivity patterns calculated from MEG recorded prior to surgery together with accurate source localization of the IED generator from EEG/MEG fusion, can provide crucial information related to surgical outcome. Resting state networks characterized by weaker connections between the IED generator and the rest of the cortex, when compared to connectivity patterns with seeds in the corresponding contralateral homologous region, indicates an isolated epileptic network, and its resection is associated with a seizurefree outcome. On the other hand, a more widespread epileptic network, characterized by stronger connectivity between the IED generator and the rest of the cortex, is a strong indicator of an unfavourable surgical outcome. These differences were mainly manifested as changes in distant connections and not in close vicinity to the seed, i.e. the IED generator or the contralateral homologous region. Connectivity analysis might bring important information towards predicting surgical outcome in patients with refractory epilepsies, especially considering that these typical MEG connectivity patterns were found at the group and individual level, in absence of any detectable epileptic discharges in scalp recordings.