Virtual localization of the seizure onset zone: Using non-invasive MEG virtual electrodes at stereo-EEG electrode locations in refractory epilepsy patients.

In some patients with medically refractory epilepsy, EEG with intracerebrally placed electrodes (stereo-electroencephalography, SEEG) is needed to locate the seizure onset zone (SOZ) for successful epilepsy surgery. SEEG has limitations and entails risk of complications because of its invasive character. Non-invasive magnetoencephalography virtual electrodes (MEG-VEs) may overcome SEEG limitations and optimize electrode placement making SOZ localization safer. Our purpose was to assess whether interictal activity measured by MEG-VEs and SEEG at identical anatomical locations were comparable, and whether MEG-VEs activity properties could determine the location of a later resected brain area (RA) as an approximation of the SOZ. We analyzed data from nine patients who underwent MEG and SEEG evaluation, and surgery for medically refractory epilepsy. MEG activity was retrospectively reconstructed using beamforming to obtain VEs at the anatomical locations corresponding to those of SEEG electrodes. Spectral, functional connectivity and functional network properties were obtained for both, MEG-VEs and SEEG time series, and their correlation and reliability were established. Based on these properties, the approximation of the SOZ was characterized by the differences between RA and non-RA (NRA). We found significant positive correlation and reliability between MEG-VEs and SEEG spectral measures (particularly in delta [0.5-4 Hz], alpha2 [10-13 Hz], and beta [13-30 Hz] bands) and broadband functional connectivity. Both modalities showed significantly slower activity and a tendency towards increased broadband functional connectivity in the RA compared to the NRA. Our findings show that spectral and functional connectivity properties of non-invasively obtained MEG-VEs match those of invasive SEEG recordings, and can characterize the SOZ. This suggests that MEG-VEs might be used for optimal SEEG planning and fewer depth electrode implantations, making the localization of the SOZ safer and more successful.

Introduction

Epilepsy surgery is a treatment option in selected patients with medically refractory focal epilepsy. It aims at rendering the patient seizure-free by removing tissue that is responsible for seizure generation and/or propagation. This region is referred to as the epileptogenic zone (EZ), and is only know unequivocally after surgery if the patient becomes seizure-free (Luders et al., 2006). Preoperatively, surrogate markers of the EZ can be obtained by non-invasive techniques. A proportion of patients need additional invasive stereo-electroencephalography (SEEG) monitoring of seizures to substantiate the hypothesis about the location of the EZ by establishing the seizure onset zone (SOZ). In resection cases, the area that is to be resected (resection area, RA) usually includes the SOZ as it is a well-established surrogate marker of the EZ (Luders et al., 2006).

Successful SOZ identification by SEEG depends upon correct planning and placement of intracerebral electrodes. The hypothesis regarding the SOZ location used for electrode placement needs to be focused, since SEEG coverage of brain structures is limited and the risk of complication rises with increasing number of electrodes implanted. Although SEEG provides essential information in selected cases, some of its disadvantages – such as high costs, risks and patient burden involved in implantation of the electrodes and the long-term invasive monitoring– apply to all patients. In addition, SEEG fails to reveal the SOZ in some cases, despite using all conventional non-invasive studies, such as video-EEG and MRI. Additional information from non-invasive studies is needed to improve the SEEG placement and to limit the number of required electrodes.

By using magnetoencephalography (MEG) for assessment of epileptiform activity, as well as for more advanced source localization and network analysis, one might overcome some SEEG limitations. With spatial filtering (beamforming) it is possible to reconstruct time series of neuronal activation at a-priori defined target locations –so called virtual electrodes (MEG-VEs) (Hillebrand and Barnes, 2005; Hillebrand et al., 2005). MEG-VEs have provided information on epilepsy, e.g. interictal spikes location concordant to the SOZ (Mohamed et al., 2013). Moreover, MEG-VEs allow detection of epileptiform discharges in deep structures, like the hippocampus, and have been used to assess functional connectivity and brain network properties in epilepsy (Hillebrand et al., 2016; Nissen et al., 2017; Nissen et al., 2016). Simulated SEEG recordings are feasible, since the location and number of MEG-VEs can be varied and MEG-VEs can be placed in user-defined locations.

Our group and others have previously used resting-state functional brain network characteristics, established from pre- and postoperative fMRI and MEG to increase prediction of success of surgery in addition to conventional interictal measures such as interictal spikes and spectral power (Quraan et al., 2013; Stam, 2014; van Diessen et al., 2013a). Epilepsy-related disturbances include both increases and decreases in mean functional connectivity, decreased long distance efficiency, altered distribution of regions with high centrality (hubs), and presence of hubs in or near the EZ (Stam, 2014; van Diessen et al., 2013a; Wendling et al., 2010). Removal of these hubs correlates with favorable outcome after epilepsy surgery (Nissen et al., 2017; Ortega et al., 2008). These measures are essentially assumption-free, meaning they involve all brain regions without focusing on the specific regions that were thought to be involved in seizure generation and propagation. The value of subnetwork analysis involving only specific brain areas needs to be established further.

To aid the planning of SEEG placement, we investigated whether interictal activity measured by MEG and SEEG at identical anatomical locations corresponded with respect to spectral, functional connectivity and network properties in patients with medically refractory epilepsy. Furthermore, we determined whether the spatial distributions of these properties were indicative of the later RA.

Methods

Patient selection

We retrospectively included nine patients with refractory epilepsy who underwent resective surgery at the VU University Medical Center. The inclusion criteria included availability of MEG and SEEG data that had been recorded as part of their presurgical evaluation, postoperative MRI and CT (Computed Tomography) with implanted depth electrodes. All patients gave informed consent to use their data for research purposes and the local Medical Ethical Review Committee at VUmc (in accordance with the Declaration of Helsinki) approved the study. One year postoperative surgery outcome was recorded using the Engel score (Engel Jr et al., 1993).

SEEG recordings

Implantation of intracerebral electrodes with multiple contacts (Ad-Tech, Medical Instrument Corporation, USA, 10–15 contacts, electrode diameter: 1.12 mm, intercontact spacing 5 mm) was performed with a stereotactic procedure and planned individually based on the hypothesized SOZ from non-invasive pre-operative studies, including video-EEG, structural MRI, MEG, Positron Emission Tomography (PET) and in selected cases ictal Single Photon Emission Computed Tomography (iSPECT). The number of electrodes per patient varied between 10 and 15 and the total number of contacts between 64 and 121, resulting in SEEG with 64–121 channels (see supplementary Table S1 for details of individual patients). Recordings were performed with a sampling frequency of 1024 Hz and downsampled by a factor of 2, except for one patient (500 Hz). From the recorded data, 10 consecutive interictal broadband (0.5–70 Hz) epochs of 4096 samples (8–8.2 s each) of the first day of the recording with the patient awake, at rest and with eyes open were selected for further analysis. No changes to anti-epileptic medication was indicated at the time of this part of the SEEG recording, hence SEEG did not differ from MEG recordings with regards to medication.

MEG recordings

A whole-head MEG system with 306 channels (102 magnetometers and 204 gradiometers, Elekta Neuromag Oy, Helsinki, Finland), placed in a magnetically shielded room (VacuumSchmelze GmbH, Hanau, Germany), was used to obtain the MEG recordings. Patients were in a supine, no-task, eyes-closed resting-state condition. A sampling frequency of 1250 Hz, an anti-aliasing filter of 410 Hz, and a high-pass filter of 0.1 Hz were used online for the recordings. Four head position indicator (HPI) coils were used to determine the head's position during the recordings. The scalp outline and HPI coil positions were digitized using a 3-D digitizer (Fastrak, Polhemus, Colchester, VT, USA). The temporal extension of Signal Space Separation (tSSS) (Taulu and Simola, 2006) (subspace correlation limit of 0.9 and a sliding window of 10 s) using MaxFilter software (Elekta Neuromag Oy, v. 2.2.15) was used to remove artifacts. The MEG data were co-registered with the patient's anatomical preoperative MRI using surface-matching software. A single sphere was fitted to the outline of the scalp and used as a volume conductor model for the beamforming approach.

MEG-VE reconstruction at SEEG electrode locations

Broadband (0.5–48 Hz) MEG-VEs were reconstructed at the locations of the SEEG contact points using a scalar beamforming method (Elekta Neuromag Oy; beamformer; v. 2.2.10), with normalized broadband beamformer weights that were computed using the broadband data covariance, a unity noise covariance matrix, and an equivalent current dipole (ECD) as the source model. The coordinates of the SEEG contact points were obtained from the post-implantation CT scan (containing the SEEG electrodes) that had been fused (linear coregistration) with the preoperative MRI scan (iPlan Net 3.0.0, BrainLab Ag, Germany) (Hillebrand and Barnes, 2005; Hillebrand et al., 2005): Fig. 1 (upper left image) shows that the SEEG contacts have a high signal density relative to brain tissue. Signal density was also higher than other high density structures, such as blood vessels. By setting a threshold, all contacts were isolated and stored. Subsequently, the x-, y- and z-coordinates for these contacts were used as input locations for the manual placement of the MEG-VEs. MEG signals at each VE location was reconstructed with beamforming: The beamforming methodology is described in detail in Hillebrand et al., 2005 (Hillebrand and Barnes, 2005; Hillebrand et al., 2005). In summary, for a particular location (in our case, the locations of the SEEG contacts) and source orientation) a set of beamformer weights are computed. The weights are chosen such that the activity at that location is reconstructed without the contribution from noise (or signal from other locations). The beamformer output at a target location for a source is the weighted sum of the output of all signal channels, and is called virtual electrode:where VE is the beamformer output or virtual electrode, W is the 1 × N weight vector and B an NxM matrix containing the magnetic field at the N sensor locations at all M latencies (Hillebrand and Barnes, 2005; Hillebrand et al., 2005). The choice of weights determines how accurate the signal for a region is reconstructed and these are determined by the formula:where Cb is the data covariance matrix, L the lead field (signal produced by the source [with unity amplitude]) and T the matrix transpose (Hillebrand and Barnes, 2005; Hillebrand et al., 2005).

Fig. 1

Overview of analysis pipeline.

This was repeated for the locations of all SEEG contacts independently, resulting in a reconstructed MEG-VE time series at every SEEG location. These resulting MEG-VE time series were downsampled by a factor two to 625 Hz. Data of 100–140 consecutive epochs of 4096 samples each (6.55 s) were selected for further analysis.

MEG and SEEG time series analysis

Fig. 1 provides an overview of the time series analysis for both MEG-VE and SEEG data. Care was taken to select epochs without ictal activity and containing as few interictal epileptiform activity (spikes and spike and wave complexes) as possible. Due to the nature of the patients included (refractory epilepsy), interictal epileptifrom activity could not be avoided in some channels in some epochs. No seizure activity was present in the selected data. The selected epochs from both MEG-VE and SEEG, with individually matched number of electrodes and locations, were analyzed with BrainWave software (v. 0.9.152.4.1) available at http://home.kpn.nl/stam7883/brainwave.html. The discrete Fast Fourier Transform (FFT) was used to obtain the peak frequency (frequency with maximum power in the range 4–13 Hz), and relative power per channel in five frequency bands: delta (0.5–4 Hz), theta (4–8 Hz), alpha1 (8–10 Hz), alpha2 (10–13 Hz), beta (13–30 Hz) and gamma (30–48 Hz) for each epoch. FFT values were averaged over epochs for each patient individually.

We estimated functional connectivity (FC) which signifies the statistical relationship between the regions' time series, by using the Phase Lag Index (PLI) (Stam et al., 2007). PLI quantifies the asymmetry in distribution of phase differences between time series (Δφ = φ1 − φ2) (Stam et al., 2007). This asymmetry relates to the proportion of phase difference (Δφ) in the interval −π < Δφ < 0 compared to that in the interval 0 < Δφ < π. No phase relationship between signals will give a symmetric distribution of differences around zero, whereas asymmetry of the phase distribution reflects a systematic phase‑leading of one channel over the other (Stam et al., 2007). PLI is obtained for each time-series (tK) with the formula:

PLI ranges between 0 and 1, where 0 indicates no coupling, or coupling with zero lag, of signals and 1 indicates total coupling or synchronization. PLI is relatively insensitive to volume conduction/field spread since it discards zero-lag (modulus π) connectivity (Stam et al., 2007; Hillebrand et al., 2012). Broadband PLI was computed between all channel pairs, yielding a functional connectivity matrix for both MEG and SEEG. PLI values were averaged over channels for each epoch, and subsequently over epochs for each patient separately.

Functional network analysis was performed by reconstructing the minimum spanning tree (MST) based on the PLI matrices for MEG-VE and SEEG (Kruskal, 1956; Stam et al., 2014). Thus, the MEG-VE and the SEEG channels served as network nodes and the inverted PLI connectivity values as undirected weighted links in the networks. The MST construction avoids arbitrary thresholding and is unique provided that the entries of the PLI matrix are unique. In addition, the number of nodes and links is the same for all networks, independent of the underlying average PLI, which enables straightforward comparison of network topology between patients and modalities. From the MST graph, measures that indicate node importance (centrality) were computed for each epoch and averaged per patient. These measures included: degree: number of links connected to the node; betweenness centrality (BC): fraction of all (shortest) paths (sequence of connections connecting two nodes) on the tree that go through the node; eccentricity (Ecc): length of the (longest) path between any two nodes (Tewarie et al., 2015). We selected the five nodes with the highest values for degree and BC separately and designated those as ‘hubs’.

Lastly, to assess the influence of interictal epileptiform activity in the resting-state data on MEG, we scored the number of spikes and spike-and-wave discharges in each MEG-VE and SEEG channel and grouped them into three categories: ‘no spikes’, ‘< 10 spikes/80 s’,’≥ 10 spikes/80 s' (80 s being the length of the SEEG selection used for the analysis).

MEG virtual electrodes (green dots in bottom left panel) were reconstructed at SEEG electrode locations (intracranial bright dots in top left panel) using beamforming. Time series were obtained per technique (MEG and SEEG), and spectral measures (peak frequency and relative power) and functional connectivity (PLI) were computed. A minimum spanning tree graph was constructed using the PLI connectivity matrix and network analysis was performed. Correlations between techniques, per measure were obtained. Comparisons between measures in the resected areas and the non-resected areas were conducted. MEG-VEs: MEG virtual electrodes; MST: minimum spanning tree; NRA: non-resected area; PLI: Phase Lag Index; RA: resected area; SEEG: stereoEEG.

Correlation, reliability and agreement between MEG-VE and SEEG findings

All statistical analyses were performed in IBM SPSS software package version 22.0.0.0. Normality was checked with the Kolmogorov-Smirnov test, significance level was set at p < 0.05. False discovery rate (FDR) correction for multiple comparison was applied (Benjamini and Hochberg, 1995).

To analyze the relationship between MEG-VE and SEEG measures, we assessed the correlation, reliability and agreement. The correlation between all MEG-VE and SEEG channels, per patient and for all measures (peak frequency, relative power, PLI, and MST measures [BC, degree, Ecc]) was estimated with the Spearman's rank correlation coefficient, yielding correlation values (rho) and p-values (determined by a Student's t-test). The correlation between MEG and SEEG was classified as weak if rho values were 0.2–0.39, moderate if 0.40–0.59, strong if 0.6–0.79 and very strong if 0.8–1.0. Interictal epileptiform discharges may have influenced the correlation analysis. Therefore, for the most significant result, and the relative delta power (which may be related to interictal epileptiform activity), we re-computed the correlation analysis with the subgroups ‘no spikes’ and ‘spikes’ to assess the effect of epileptiform activity on the correlation.

The reliability of all measures was examined with the intraclass correlation coefficient (two-way mixed ICC). An ICC < 0.2 was considered poor, 0.21–0.4 fair, 0.41–0.6 moderate, 0.61–0.8 good and 0.81–1 excellent reliability. Bland-Altman plots were used to graphically show level of agreement between modalities (Bland and Altman, 1986) for the group mean across channels per measure. Bias, precision and 95% limits of agreement were calculated where the mean difference in measure values between both modalities is depicted as the “bias” and the standard deviation of the mean difference is depicted as the “precision”. Both modalities agreed and the bias was not considered relevant if the differences between MEG-VE and SEEG-based measures fell within the 95% limits of agreement.

Agreement between MEG-VE and SEEG with regard to occurrence of epileptiform activity was assessed by Kendall rank correlation coefficient.

Resected versus non-resected area comparisons

The resection area was manually segmented on the postoperative MRI, and the postoperative MRI and the post-implantation CT were linearly registered with the preoperative MRI (iPlan Net 3.0.0, BrainLab Ag, Germany). Each channel (in MEG-VE and SEEG) was labeled ‘resected area’ (RA) when the anatomical localization of the channel was located inside the resection area, and labeled ‘non-resected area’ (NRA) when the channel was located outside. For group comparisons, mean values were obtained per area (RA and NRA) for each measure (peak frequency, relative power, PLI, and MST measures [BC, degree, Ecc]). The mean values of RA and NRA were compared using a Wilcoxon signed rank test, because the data were not normally distributed. Within each patient, comparisons between areas (RA vs. NRA) were performed for the statistically significant results in the group analysis using an independent sample Mann-Whitney U test.

Additionally, to assess if hubs were located in the RA, the correspondence of any of the five highest BC or degree values with a RA channel was reported as a YES/NO answer.

Results

Patient characteristics are shown in Table 1.

Table 1

Patient characteristics. F: female; FCD: focal cortical dysplasia; M: male; MTS: mesial temporal sclerosis.

Patient	Sex	Age (Y)	MRI findings	Resection area	Engel Score
1	F	35	FCD, MTS	Left anterior temporal and amygdala/hippocampus	I A
2	M	48	MTS	Left anterior temporal and amygdala/hippocampus	I A
3	M	21	FCD	Right frontal	I A
4	F	49	Hamartoma	Left lateral temporal	I A
5	M	40	No lesion	Left anterior temporal and amygdala/hippocampus	II A
6	F	29	Multiple cavernoma	Left anterior temporal and amygdala/hippocampus	I A
7	F	31	FCD	Right frontal	I A
8	M	16	Polymicrogyria	Right frontal and parietal bordering the central and postcentral sulcus	IVB
9	F	39	Porencephalic cyst	Right temporal, hippocampus, amygdala	IA

Table 2 and supplementary Table S2 show the correlations and reliabilities between MEG-VE and SEEG. Peak frequency had a significant correlation (Spearman's rho) and reliability (fair to good ICC) over modalities in seven and eight out of nine subjects, respectively. Correlations and reliabilities between MEG-VE and SEEG relative powers differed per patient and per frequency band: Relative delta, alpha2 and beta power correlated significantly in five to seven out of nine patients. These bands also showed the highest reliability (Table 2). Functional connectivity correlated significantly and showed good reliability between MEG-VE and SEEG in six patients. Five out of the seven patients with significant correlation in peak frequency between MEG-VE and SEEG also showed significant PLI correlation. Of the two patients who continued to have seizures after surgery, patient 5 showed good correlation and reliability between modalities in PLI and beta relative power and patient 8 for peak frequency. Neither consistent correlation nor reliability was found for the network measures (degree, BC and Ecc).

Table 2

Correlation and consistency between MEG virtual electrodes and stereo-EEG.

Patient	CC [rho (p)]					ICC
	PF	PLI	Deltaa	Alpha2a	Betaa	PF	PLI	Deltaa	Alpha2a	Betaa
1	0.49b (<0.001)	0.42b (<0.001)	0.32 (0.001)	0.42b (<0.001)	0.54b (<0.001)	0.46b	0.55b	0.41b	0.74b	0.57b
2	0.32 (<0.001)	0.01 (0.34)	0.11 (0.29)	−0.002 (0.98)	−0.05 (0.59)	0.58b	0.11	0.39	0.21	0.42b
3	−0.10 (0.3)	0.02 (0.9)	0.19 (0.07)	−0.03 (0.78)	0.28 (0.004)	−0.13	−0.03	0.24	−0.25	0.23
4	0.67b (<0.001)	0.36 (0.003)	0.75b (<0.001)	0.30 (0.02)	0.62b (<0.001)	0.66b	0.58b	0.66b	0.27	0.65b
5	0.19 (0.06)	0.46b (<0.001)	0.16 (0.1)	0.21 (0.04)	0.57b (<0.001)	0.25	0.42b	0.15	0.17	0.68b
6	0.66b (<0.001)	0.44b (<0.001)	0.36 (<0.001)	0.72b (<0.001)	0.19 (0.06)	0.58b	0.47b	0.53b	0.74b	0.55b
7	0.20 (0.04)	0.28 (0.003)	0.66b (<0.001)	0.36 (<0.001)	0.79b (<0.001)	0.35	0.33	0.61b	0.34	0.84b
8	0.45b (<0.001)	0.18 (0.07)	0.11 (0.28)	0.07 (0.49)	−0.19 (0.05)	0.51b	0.07	0.22	0.32	−0.17
9	0.53b (<0.001)	0.54b (<0.001)	0.72b (<0.001)	0.40b (<0.001)	0.77b (<0.001)	0.52b	0.54b	0.72b	0.46b	0.75b

CC: correlation coefficient; ICC: intraclass correlation; MEG-VE: MEG virtual electrodes; PF: peak frequency; PLI: Phase Lag Index; SEEG: stereo-EEG.

Relative power.

Values indicating moderate to strong correlation (CC) or moderate to excellent reliability (ICC).

Table 3 shows the distribution of the MEG-VE and SEEG channels over the categories of the interictal epileptiform activity.

Table 3

Distribution of MEG virtual electrodes and stereo-EEG channels over the categories of interictal epileptiform activity.

SEEG spikes
		No spikes	<10 spikes/80 s	≥10 spikes/80 s	Total
MEG-VE spikes	No spikes	770	34	33	837
	<10 spikes/80 s	141	0	2	143
	≥10 spikes/80 s	72	3	0	75
Total		983	37	35	1055

As can be seen from the table, most channels did not contain spikes. Kendall's tau was −0.09 (p < 0.001), but this was largely due to the channels containing no spikes in MEG-VE as well as SEEG. When these channels were discarded from the analysis, the correlation was not significant anymore (tau = 0.08, p = 0.30). Because of this distribution, we merged the spike categories into one ‘spike group’ for further analysis.

We found no systematically increased correlations in the channels containing spikes compared to the channels without spikes for peak frequency and relative delta power: Two subjects did not have spikes in the selected epochs and therefore were not included in the analysis. For peak frequency, three out of seven subjects had stronger correlation in the ‘spike group’ than in the ‘no spike group’, for two patients the correlation was weaker, and for one patient there was no difference between the two groups of channels. For delta power, stronger correlation in the ‘spike group’ was found in four of seven subjects. Supplementary table S6 shows correlation values for peak frequency and relative delta power for ‘spike’ and ‘no spike’ groups separately.

Bland-Altman plots are shown in Fig. 2 for peak frequency, delta and alpha relative power and PLI (based on strong correlations and good reliabilities; data for the other measures can be found in supplementary Fig. S1). All spectral and connectivity measures (but not network measures) fell within the limits of agreement, except one PLI and relative theta power correlation, therefore the agreement between the two modalities is good. For relative beta power and PLI we found indications for a proportional error (larger mean values lead to larger differences between modalities), although for PLI the differences were small. We found a bias for peak frequency and relative power in the high frequency bands (alpha1, alpha2, beta and gamma) with MEG-VE being systematically higher than SEEG. For the relative power in the low frequency bands (delta and theta) and PLI, this was the opposite.

Fig. 2

Bland-Altman plots for the differences between SEEG and MEG-VE for peak frequency (A), relative delta and alpha2 power (B and C, respectively), and functional connectivity (D). Bland-Altman plots were used to graphically show level of agreement between modalities for the group mean across channels per measure. Agreement of both techniques and the bias was not considered relevant if the differences between MEG-VE and SEEG-based measures fell within the 95% limits of agreement. Note that all except one PLI difference fall within limits of agreement (D). Blue lines: limits of agreement mean ±1.96 SD; black line: mean difference between MEG and SEEG (mean); dotted black line: no difference; yellow dotted lines: 68% confidence interval (CI) of the mean. LA: limits of agreement.

Resected versus non-resected areas

Regarding group comparisons, significant differences in mean relative power between RA and NRA were found for both MEG-VE and SEEG, with relatively more slow activity (delta power) in the RA and relatively less fast activity (alpha1 and 2, beta and gamma power) than in the NRA (Fig. 3A, supplementary Table S3). The spectral differences were also reflected in differences in peak frequency between RA and NRA (Fig. 3B, supplementary Table S3), with significant lower peak frequency in the RA compared to the NRA for both MEG-VE (Z = −2.55, p = 0.01) and SEEG (Z = −2.55, p = 0.01). Functional connectivity was higher in the RA in MEG-VE and SEEG, but the differences with the NRA were not statistically significant after correction for multiple comparisons (MEG-VE: Z = −1.72, p = 0.05; SEEG: Z = −1.96, p = 0.09) (Fig. 3C, supplementary Table S3).

Fig. 3

Resected area (RA) and non-resected area (NRA) differences for MEG-VE and SEEG. MEG virtual electrode (MEG-VE; red) and stereo-EEG (SEEG; blue) relative power in RA and NRA for different frequency bands (A), peak frequency (B) and Phase Lag index (PLI; C).

*: significant difference after correction for multiple comparisons. Note that RA has higher relative delta power and lower alpha2 power (A) and peak frequency (B) for both MEG-VE and SEEG. Functional connectivity (PLI) tended to be higher in the resected area, although non-significant after multiple comparison correction (C). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Post-hoc analysis at the individual level was performed for delta and alpha2 power and peak frequency, because these measures were significantly different between RA and NRA for both MEG-VE and SEEG (other frequency bands were only significant for one modality). Table 4 shows the differences in peak frequency between RA and NRA for MEG-VE and SEEG for each patient separately. Differences were significant for six patients in MEG-VE and for five patients in SEEG. Results for relative delta and alpha2 power are given in supplementary Tables S4 and S5, showing significantly higher RA delta power in three and six patients (MEG-VE and SEEG, respectively), and with significantly lower RA alpha2 power in six and five patients (MEG-VE and SEEG, respectively). The two patients who were not seizure-free after operation showed significant lower peak frequency and higher delta power in SEEG RA channels, but not in MEG-VE RA channels. For alpha2 power, only one patient of those two showed decreased alpha2 power in RA channels either for MEG-VE or SEEG.

Table 4

Peak frequency differences between resected and non-resected areas per patient.

Patient	MEG-VE		SEEG
	RA	NRA	RA	NRA
1	6.70(0.28)	7.53⁎(1.00)	6.09(0.38)	7.19(1.26)
2	6.70(0.43)	7.12⁎(0.55)	5.36(0.57)	5.31(0.73)
3	5.27(0.45)	6.55⁎(0.45)	4.76(0.25)	4.83(0.37)
4	6.76(0.66)	8.04⁎(0.73)	6.64(1.02)	7.87⁎(1.45)
5	7.45(0.20)	7.60(0.40)	6.35(0.71)	6.96⁎(1.23)
6	5.98(0.44)	7.49⁎(0.71)	4.99(0.57)	6.53⁎(1.64)
7	7.25(0.30)	7.23(0.41)	5.08(0.42)	5.99(0.73)
8	6.57(0.99)	6.43(1.27)	4.84(0.35)	5.22⁎(0.58)
9	6.20(0.47)	6.59⁎(0.48)	5.86(0.80)	6.43⁎(0.94)

MEG-VE: MEG virtual electrodes; NRA: non-resected area; RA: resected area; SEEG: stereo-EEG. Values are mean peak frequency (SD).

p < 0.05, corrected for multiple comparisons.

Network measures did not differ significantly between RA and NRA for either modality (supplementary Table S3). The additional analysis of functional network hubs showed that in six of the nine patients (67%), an RA channel was among the five highest values for both centrality measures (BC and degree) for both modalities (except BC of MEG-VE that corresponded in seven patients [78%]). There was no significant difference in the mean Ecc between RA and NRA for either modality.

Discussion

The purpose of this study was to determine whether a new method with interictal MEG-VE reconstructed at the same location as SEEG electrodes might facilitate the epilepsy surgery workup. We correlated spectral, connectivity, and network measures from MEG-VE time series with those obtained from interictal SEEG recordings. In addition, we aimed at identifying markers for the resected areas, as such might assist SEEG electrode planning. We show that the invasively and non-invasively determined measures are correlated and reliable and that spectral analysis measures (relative power and peak frequency) are significantly different between the RA and the NRA for both modalities at group level, although not consistently at the individual level. Furthermore, we show that network hubs overlap with the RA, but the results do not reach statistical significance.

Correlation, reliability and agreement between modalities

In this study, we used per patient the same location for the MEG virtual electrodes as for the SEEG electrodes, thereby enabling direct comparison of the two modalities. Spectral measures and functional connectivity, but not network measures, showed significant correlations, reliability and agreement between MEG-VE and SEEG. These findings are valuable given that spectral analysis (and a non-significant trend using connectivity) was also shown to be sensitive in detecting changes between the RA and the NRA and that simultaneous recording is not necessary for this correlation. Taken together, this shows that significant differences between RA and NRA are picked up by measures that also show reliability over modalities, and over time, and that consistent RA versus NRA differences could not be detected by measures that were also not consistent over modalities. It seems that for network measures to be useful in identifying the RA, they need to be relatively modality-independent.

To our knowledge this is the first study comparing MEG-VE and SEEG in terms of spectral analysis and functional networks. Previous studies comparing the two modalities in terms of spike locations using non-simultaneous recordings also found good correspondence between the modalities (Bouet et al., 2012; Grova et al., 2016; Murakami et al., 2016). Bouet et al. showed good correlation between spikes on MEG and SEEG when using a beamformer method (Bouet et al., 2012). Additionally, it has been shown that concordance between MEG and SEEG spike localization was related to seizure freedom after surgery (Murakami et al., 2016). Overall, in these studies, MEG was compared to the “gold standard” SEEG. However, since SEEG electrode placement may have been dependent on the results of MEG spike localizations, the presented correlations may be artificially high as there could have been reduced SEEG electrode coverage in MEG areas without epileptiform abnormalities. Although in the current study SEEG placement was not based on the parameters under study (spectral, connectivity and graph measures), the data were not completely free of interictal epileptiform activity. Therefore, our results might also have been influenced by the presence of spikes and spike-and-wave discharges in the selected epochs. However, we found no correlation between the qualitatively assessed amount of epileptiform activity in MEG-VE and SEEG in the time series containing interictal epileptiform activity, possibly because we attempted to avoid such activity during the epoch selection. This was done since it was not the topic under study and we aimed at minimizing the effect of interictal epileptiform activity on the other analyses. This might also be the explanation for the incongruence with the study of Bouet et al., where spikes were investigated specifically. When assessed separately, we found no indications for a systematically higher correlation between MEG-VEs and SEEG for relative delta power or peak frequency in channels with spikes, compared to the spike-free channels. Overall, we therefore did not find indications for a large effect of interictal spikes on the identified correlations between MEG-VE and SEEG data.

Simultaneously recorded MEG and SEEG have repeatedly confirmed the relationship between spikes in the two modalities (Badier et al., 2017; Gavaret et al., 2016; Kakisaka et al., 2012). In our study, we showed that non-simultaneous and interictal activity can still confirm the correlation between MEG and SEEG, at least when sampled from the same locations.

Characterization of the resected area using interictal MEG-VE and SEEG

Spectral analysis showed that the interictal activity in the RA was significantly slower compared to the NRA. Previous studies have shown that dipole modeling of high density EEG and MEG delta activity is related to lesion location on MRI in the majority of patients with lesional epilepsy (Baayen et al., 2003). Asymmetric slowing in MEG was found to reliably lateralize the resected area, although localizing value was limited (Englot et al., 2016a; Ishibashi et al., 2002). Our observations from MEG-VE and SEEG are in line with these results, and in addition provide localizing information. In general, spectral contrast between brain areas that are included in SEEG exploration may be lower than between time series that would be used in a whole brain analysis because of their suspected common involvement in seizure generation. Nevertheless, our study indicated that, even within this subset of selected brain areas, significant slowing can be found in the areas that were later involved in the resection compared to the non-resected areas.

Functional connectivity was higher in the RA compared to the NRA in our study, although not significantly so after correction for multiple comparisons. Previous studies report mainly functional connectivity increases in the epileptogenic zone (Nissen et al., 2016; Bartolomei et al., 2017; Englot et al., 2016b), and one study relates higher functional connectivity in the resection area to better surgery outcome (Englot et al., 2015). Together, these results indicate that the epileptogenic zone is well-connected in the interictal state (Nissen et al., 2016; Bettus et al., 2008; Maccotta et al., 2013).

Graph theory, combined with modern network analysis, has shown value in epilepsy surgery related research (for a review see (Bartolomei et al., 2017)). In our study, the RA included highly central nodes in the majority of patients, suggesting that hubs are often located in the EZ. Focusing on MEG and SEEG networks, some studies have found that hubs were related to the EZ or resection area (Nissen et al., 2017; Dickten et al., 2016; Li et al., 2016; Varotto et al., 2012; Vlachos et al., 2016). However, results have not been consistent since other studies reported low centrality within the EZ (Nissen et al., 2016; van Diessen et al., 2013b). Different centrality measures and directed connectivity studies may add to understanding the relation between hubs and the EZ.

Limitations

Our study was retrospective in nature. The challenge will be to utilize MEG-VE prospectively and evaluate the feasibility of achieving the same results and accuracy in localizing the EZ as with SEEG. Other limitations are the small sample size and the variability within (different number of channels included in the RA vs. the NRA) and among patients, which limited the statistical power of our findings. As noted above, MEG recordings typically precede SEEG recordings by several months in our patient population, therefore differences between the modalities could have been caused by differences in the patient's clinical characteristics at the time of the recordings, as well as their state (e.g. changes in the level of arousal at recording times, circa- and ultradian rhythm variations, serum AED levels, occurrence of seizure clustering). However, despite these factors, we were still able to find RA related slowing in both modalities, as well as high correlation and reliability between MEG-VE and SEEG.

Two patients were not seizure free after surgery. Therefore, the resection area did, by definition, not include the epileptogenic zone in these patients and therefore we cannot extrapolate our group findings to the EZ. At the electrode-level, MEG-VE and SEEG correlated equally well in the seizure-free and the not seizure-free individuals. We therefore conclude that MEG-VE and SEEG analysis of interictal activity correspond to each other, irrespective of surgery outcome.

The correlation coefficients per measure have to be interpreted with care given that they were performed among all the channels (without distinction between RA and NRA channels), as the correlation could be inflated due to differences between RA and NRA channels. A difference in correlation compared to when taking channel grouping (RA and NRA) into consideration might imply that RA and NRA channels differed or that the correlation among all the channels was not entirely linear. As an additional test, ICCs were performed and the measures derived from the two modalities proved to be reliable, suggesting that the correlations were not only driven by differences between RA and NRA channels. Finally, the accepted limits for bias and precision between modalities were subjective, and might therefore vary across clinicians and centers.

Conclusion and future directions

Currently, non-invasive interictal MEG cannot replace invasive ictal SEEG. Interictal abnormalities, such as spikes and spike and wave complexes, are not sufficient to identify the SOZ and ictal recordings are still needed. Our results indicate that by using virtual MEG electrodes at the anatomical locations of the SEEG electrodes, invasive time series can be non-invasively approximated and thereby reduce the gap between interictal MEG and ictal SEEG. These results open up the possibility of using non-invasive MEG to “predict” later SEEG findings at various cerebrocortical locations, and to use this information for optimal SEEG electrode placement. Additionally, this study paves the way for further analyses on e.g. spike location, high frequency oscillations, and directed network centrality. Furthermore, since MEG provides whole head recordings, several VE strategies might be tested for alternative SEEG implantation schemes based on the anatomico-electro-clinical hypothesis regarding the seizure onset and propagation in an individual patient.