Functional connectivity of the irritative zone identified by electrical source imaging, and EEG-correlated fMRI analyses.

OBJECTIVE: The irritative zone - the area generating epileptic spikes - can be studied non-invasively during the interictal period using Electrical Source Imaging (ESI) and simultaneous electroencephalography-functional magnetic resonance imaging (EEG-fMRI). Although the techniques yield results which may overlap spatially, differences in spatial localization of the irritative zone within the same patient are consistently observed. To investigate this discrepancy, we used Blood Oxygenation Level Dependent (BOLD) functional connectivity measures to examine the underlying relationship between ESI and EEG-fMRI findings.
METHODS: Fifteen patients (age 20-54), who underwent presurgical epilepsy investigation, were scanned using a single-session resting-state EEG-fMRI protocol. Structural MRI was used to obtain the electrode localisation of a high-density 64-channel EEG cap. Electrical generators of interictal epileptiform discharges were obtained using a distributed local autoregressive average (LAURA) algorithm as implemented in Cartool EEG software. BOLD activations were obtained using both spike-related and voltage-map EEG-fMRI analysis. The global maxima of each method were used to investigate the temporal relationship of BOLD time courses and to assess the spatial similarity using the Dice similarity index between functional connectivity maps.
RESULTS: ESI, voltage-map and spike-related EEG-fMRI methods identified peaks in 15 (100%), 13 (67%) and 8 (53%) of the 15 patients, respectively. For all methods, maxima were localised within the same lobe, but differed in sub-lobar localisation, with a median distance of 22.8 mm between the highest peak for each method. The functional connectivity analysis showed that the temporal correlation between maxima only explained 38% of the variance between the time course of the BOLD response at the maxima. The mean Dice similarity index between seed-voxel functional connectivity maps showed poor spatial agreement. SIGNIFICANCE: Non-invasive methods for the localisation of the irritative zone have distinct spatial and temporal sensitivity to different aspects of the local cortical network involved in the generation of interictal epileptiform discharges. Crown

Introduction

During presurgical evaluation and surgical planning in patients with refractory epilepsy, interictal epileptiform discharges (IEDs) are used to indicate tissue involved in epileptic networks. The irritative zone, the cortical substrate generating interictal epileptiform discharges, can be detected in a variety of ways (Rosenow, 2001). Simultaneous spike-related electroencephalography and functional magnetic resonance imaging (EEG-fMRI) maps cortical and subcortical blood flow changes associated with IEDs with high spatial resolution. However, up to 70% of EEG-fMRI acquisitions are inconclusive because of sparsity of IEDs during the examination (Salek-Haddadi et al., 2006, Al-Asmi et al., 2003). To counter this, voltage maps of IEDs recorded outside the MRI scanner have been correlated with voltage maps obtained during EEG-fMRI and the time course of the correlation between voltage maps used as the regressor in the fMRI analysis. The voltage-map approach can demonstrate fMRI activation in the absence of IEDs during the scanning session and has been shown to increase the sensitivity of EEG-fMRI for detecting the irritative zone (Grouiller et al., 2011). Electrical source imaging (ESI) identifies cortical sources of IEDs using a combination of high-density EEG recording (≥64 electrodes), accurate forward models and distributed inverse solutions that yield more accurate results than dipole localisation methods (Tanaka et al., 2009, Tanaka et al., 2014, Pellegrino et al., 2018, Tamilia et al., 2019)

Previous studies assessing the clinical contribution of individual non-invasive methods for identifying the irritative zone have observed good spatial concordance with intracranial EEG (Benar et al., 2006, Grouiller et al., 2011). However, comparisons between methods have revealed that identified foci in the brain are not identical (Elshoff et al., 2012, Vulliemoz et al., 2010, Centeno et al., 2017). This has implications for the clinical utility and interpretation of these methods. Moreover, it has been proposed that epilepsy involves several abnormal brain networks (Spencer, 2002, Bartolomei et al., 2008). We hypothesised that different methods are sensitive to different networks involved in the generation of IEDs. To test this hypothesis we investigated the functional and temporal relationship between foci of the irritative zone identified by each method. Our aim was to quantify the BOLD temporal time course correlation and the spatial similarity of functional connectivity maps. The spatial discrepancy at the lobar and sub-lobar level was quantified and the functional connectivity between peaks identified by each method evaluated. We compared the outcome with functional connectivity measures from seeds outside the irritative zone in patients undergoing presurgical evaluation.

Methods

Subjects

Seventeen patients (median age 26 years; range 18 to 54 years) with IEDs over at least 16 h of continuous video-EEG recording were initially recruited from the Royal Brisbane and Women’s Hospital (Brisbane, Australia). Two patients were excluded from the study; one was unable to complete the fMRI study because of a seizure during the scanning session and one had residual EEG artefact after following the pre-processing steps. The patients were scanned the week after completing the video-EEG investigation, while on their usual antiepileptic medications. Demographic and clinical information of the patients are listed in Table 1. The study was approved by the respective hospital and university Human Ethics Committees. All participants provided written informed consent.

Table 1

Patient demographics and pre-surgical clinical work-up information. SE: History of status epilepticus. AED: Antiepileptic medication at the time of the EEG-fMRI. LTEEG: Long-term EEG monitoring. FIA: Focal Impaired Awareness. GTCS: Generalised Tonic-Clonic Seizure. FA: Focal Aware Seizures. L: Left. R: Right. FLE: Frontal Lobe Epilepsy. TLE: Temporal Lobe Epilepsy. PLE: Parietal Lobe Epilepsy. LVT: Levetiracetam; LTG: Lamotrigine; LAC: Lacosamide; CBZ: Carbamazepine; CLB: Clobazam; ZON: Zonisomide; OXC: Oxcarbazepine; PER: Perampanel; PRG: Pregabalin; TPM: Topiramate.

#	Age/gender	Age onset	AED	Seizure type	Epilepsy type	Spike max. electrode	MRI	PET	SPECT
P01	43 f	35	PHT, LTG, CLB, LVT, LAC	FIA, GTCS	L TLE	F7, T3, A1	L HS	L TL, PL and OL hypometabolism	hyperperfusion L antero mesial TL
P02	33 f	5	OXC, TPM	GTCS	FLE	T5, T3	normal	L P hypometabolism	N/A
P03	25 f	18	CLB, LAC, LTG, ZON	FA, FIA, GTCS	L FLE, TLE	F7, T3	L anterior TL resection	Diffuse L posterior PL hypometabolism	Hyperperfusion in the inferolateral aspect of the L TL
P04	26 f	12	CBZ, PRG, LAC	FA, FIA, GTCS	R/L FLE, TLE	T3, T5, A1	normal	Mild hypometabolism R TL	N/A
P05	54 m	3	LVT, CBZ, TPM	FIA, GTCS	L PLE	Cz, C3, Pz, P3	L P dysplasia	N/A	Hyperperfusion on the lateral aspect of the L PL lesion
P06	19 f	17	LVT, CLB, VPA	FA, FIA, GTCS	R TLE	F8, A2	Dysplastic lesion R amygdala	N/A	N/A
P07	23 m	12	VPA, TPM	FIA, GTCS	R F	FP2	R F gliosis	R FL photopenia	N/A
P08	26 m	11	VPA, FYC, LVT, TPM, VIM	FA, FIA, GTCS	R multi-focal	T4, F8, A2	normal	R posterior F hypermetabolism	R posterior FL hotspot
P09	53 m	32	LVT, LAM, TPM	FA, FIA, GTCS	R TLE	T4, A2	R HS	R TL hypometabolism	N/A
P10	20 m	15	LVT, VIM, CBZ	FA, FIA, GTCS	R PE	P4	normal	R PL hypometabolism	N/A
P11	18 f	3	LVT, CBZ	FA, FIA	R TLE	F8, T4	R HS	normal	N/A
P12	47 f	3	VPA, CBZ, PER	FA, FIA, GTCS	R TLE	F8	R HS	Bitemporal hypometabolim > R hipp. and R Insula	R Insula hyperperfusion
P13	40 m	30	VAP	FA, FIA, GTCS	L TLE	T3, A1	normal	normal	N/A
P14	23 m	11	CBZ, LVT, PER, TPM	FA, FIA	L PLE	T5, T3	L PL lesion (DNET)	N/A	Perilesional hotspot in the L PL
P15	22 f	5	LTG, TPM	FA, FIA	L TLE	T3, F7	L HS	L TL hypometabolism	N/A
P16	31 f	26	LTG, BS, LVT, VPA	FA, FIA, GTCS	R TLE	F8, A2	R anterior T dysplasia	Diffuse hypometabolism R anterior and mid TL	R anterior TL hotspot
P17	47 f	39	LTG	FA, FIA	R TLE	A2, T4	normal	mild-hypometabolism in the R TL	N/A

EEG-fMRI acquisition

A 64-channel MRI-compatible EEG cap (Waveguard, ANT Neuro, The Netherlands) was placed according to the 10–20 system with an impedance below five kΩ for each electrode. An ECG electrode was placed on the left thoracic back to allow correction of the ballistocardiogram artefact (Allen et al., 1998). The EEG was recorded at 5 kHz using two 32-channel BrainAmp MRplus amplifiers (Brain Products GmbH) synchronised to the MRI system clock to allow the correction of MRI gradient artficats (Allen et al., 2000). Functional and structural images were obtained on a 3 T Magnetom Trio Scanner (Siemens, Erlangen, Germany) with a 32-channel head coil for patient P01 and P02, and on a 3 T Prisma fit using a 64-channel head coil (Siemens, Erlangen, Germany) for all other subjects.

For every patient, the fMRI data were obtained in 5 runs containing 36 slices per volume covering the full brain. In total, 750 fMRI volumes were recorded in a single 30 min EEG-fMRI session using a 2D single-shot T2*-weighted gradient echo-planar sequence (TR = 2.34 s, TE = 30 ms, flip angle = 90°, matrix = 64 × 64, FOV = 210 mm with an isotropic voxel size of 3.3 mm). The patient was instructed to remain still, relaxed with their eyes closed. A structural three-dimensional T1 self bias-field corrected Magnetisation Prepared 2 Rapid Acquisition Gradient Echoes (MP2RAGE) was also acquired (Marques et al., 2010), covering the entire brain with a resolution of 1 mm isotropic voxel size, 240 × 256 × 176 matrix size, GRAPPA = 3, TR = 4000 ms, TE = 2.89 ms with the inversion times of 700 and 2220 ms.

EEG-fMRI processing

Functional images were preprocessed using SPM12 (Wellcome Trust Centre for Neuroimaging, London, http://www.fil.ion.ucl.ac.uk/spm). The pre-processing steps included removing the first four fMRI volumes, realignment, slice timing correction, coregistration to the structural 3D T1 MP2RAGE image, spatial normalization, and smoothing using an isotropic Gaussian kernel of 8 mm full width at half-maximum. The pre-processed fMRI time series were analysed using a general linear model approach. Nuisance regressors in both analyses included covariates obtained from the six parameter realignment step and from the cerebrospinal fluid and white matter timecourses.

The EEG recorded inside of the scanner was preprocessed using BrainVision Analyzer (version 2.1, Brain Products GmbH). Time-locked gradient and ballistocardiogram artefacts were removed using the average template subtraction method (Allen et al., 1998). EEG data were then down-sampled to 500 Hz and a 1 Hz high-pass filter was applied. Residual temporal components of artefacts were removed by visual inspection of the topographic and temporal components obtained using Infomax Independent Component Analysis. Data were then low-pass filtered using an infinite response filter with a cut-off frequency of 35 Hz and segmented to discard the EEG recording from the first 4 MRI volumes. The preprocessed EEG data were then exported to European Data Format to perform the different types of analyses detailed in the following section.

Localisation of the irritative zone

Three EEG-based analysis methods were used for the localisation of the irritative zone at the lobar and sub-lobar level. As measurements of the irritative zone can vary depending on several factors including sampling period, the number of seizures at the time of the scan, state of wakefulness and level of antiepileptic medication (Rosenow, 2001), in this study all analyses were performed on simultaneously acquired datasets. For each method used in this investigation, a single localisation at the global maximum was selected for further analysis. We chose the maximal source strength and the BOLD global maximum to assess the spatial concordance between ESI and EEG-fMRI, respectively as these have been previously been shown to map to the irritative zone for both ESI (Brodbeck et al., 2011, Birot et al., 2014) and EEG-fMRI (Thornton et al., 2010, Khoo et al., 2018, An et al., 2013). The procedure for analysing each method investigated was:

Spike-related EEG-fMRI analysis

The peaks of in-scanner spikes were marked in consecutive order and were treated as delta functions which were convolved with the canonical hemodynamic response function and its temporal and dispersion derivatives as implemented in SPM12. In patients with multifocal IEDs, the most frequent IED reported during clinical long-term EEG investigation was selected for further analysis

Voltage-map EEG-fMRI analysis

We followed the pre-processing steps described by Grouiller et al. (2011). In brief, we calculated the spatial correlation between the patient-specific voltage map of IEDs obtained during video-EEG monitoring with each timepoint of the EEG recorded during the simultaneous EEG-fMRI acquisition (Grouiller et al., 2011). The resulting correlation vector was squared, convolved with the canonical HRF and used as the regressor of interest along with time and dispersion derivatives for the voltage map EEG-fMRI analysis. Video-EEG recordings were performed using 21 scalp electrodes placed according to the 10–20 electrode system recorded at 500 or 2000 Hz on an E-Sleep 443 or Neuvo amplifiers (Compumedics, Melbourne, Australia). IEDs were marked on artefact free periods and were averaged after a 35 Hz infinite response lowpass filter had been applied to match in-scanner EEG settings. Average voltage maps at the peak of the Global Field Power, when signal-to-noise ratio is highest and topography is most stable (Lehmann, 1984), were resampled using spline interpolation to match the 64 electrode system of the in-scanner EEG acquisition. The resulting voltage map was used to calculate spatial correlations with in-scanner voltage maps.

Both spike-related and voltage-map statistical parametric maps were thresholded at p < 0.05 with False Discovery Rate correction for multiple comparisons and a cluster threshold of 0.001 across the whole brain. Only positive BOLD responses were analysed as previous EEG-fMRI studies investigating the clinical utility of EEG-fMRI for epilepsy localisation have demonstrated a higher concordance between the epileptogenic zone and BOLD activation as opposed to deactivation (Federico et al., 2005, Zijlmans et al., 2007, Pittau et al., 2012).

Electrical source imaging

Source localisation maps were generated using IEDs acquired during the EEG-fMRI acquisition when these were available. For patients without IEDs during the fMRI scan, the averaged spike used for the voltage map analysis was used. In both cases, the inverse solution was calculated at 50% of the rising phase of the averaged spike. This time point was selected as it has been demonstrated to represent regions of origin rather than propagation of spikes (Lantz et al., 2003, Ray et al., 2007).

To estimate the localisation of the electrical source, we used a locally spherical head model with anatomical constraints (LSMAC; Brunet et al., 2011) based on the local radiuses of the scalp and skull in the patient’s structural MRI scan. The latter was segmented into gray and white matter and the solution space was confined to gray matter regions in which 5000 evenly distributed solution points were defined. The spatial coordinates of the scalp electrodes were obtained from the T1 weighted structural image. The inverse problem was solved using the distributed inverse solution local autoregressive average (LAURA; Grave de Peralta Menendez et al., 2001). The results of the inverse solution were mapped to the fMRI volume using cubic spline interpolation. The ESI analysis and export of each inverse solution was performed using Cartool software v3.7 (Brunet et al., 2011), following current best practice steps to estimate the distribution of neuronal sources derived from multichannel EEG recordings (Michel and Brunet, 2019).

Spatial concordance and distance between global maxima

Comparisons between peaks obtained from the three methods in different subjects were performed in Montreal Neurological Institute (MNI152) space following standard non-linear anatomical normalisation as implemented in SPM12. Automated Anatomical Labelling (AAL) (Tzourio-Mazoyer et al., 2002) and Brodmann Areas as implemented in WFU PickAtlas toolbox (Maldjian et al., 2003) were used for anatomical parcellation at the lobar (Frontal, Temporal-Insular, Parietal and Occipital cortices) and sub-lobar levels. Euclidean distances between peaks identified with ESI and spike-related and voltage-map EEG-fMRI were calculated in Matlab R2018b (The MathWorks Inc., MA, USA). Statistical differences between Euclidean distances were assessed using the Wilcoxon-rank sum test at a significance level of p < 0.05.

Seed-based functional connectivity analysis

The functional connectivity of peaks determined by each method was assessed using seed-based functional connectivity analysis, as implemented in the CONN Toolbox v18b (https://www.nitrc.org/projects/conn/). Seed regions of interest (ROIs) comprised spherical ROIs 9 mm in diameter centred on each peak. Standard first-level denoising steps as implemented in CONN, including linear detrending, despiking and smoothing were employed. Additionally, cerebrospinal fluid, white matter, motion parameters extracted from the SPM12 realignment step and their derivatives were considered as confounds and regressed out. Prior to the functional connectivity analysis, the residual time series was low-pass filtered at 0.008 Hz. Individual correlation (bivariate) maps were obtained from each seed to each voxel across the whole brain and transformed into individual Fisher-z maps. The significance level was set at a FDR < 0.05.

Spatial similarity analysis

The spatial overlap between connectivity maps was assessed using the Dice similarity index (Zou et al., 2004)where X and Y represent the number of voxels of interest in two image datasets and X ∩ Y denotes the number of voxels of interest in common between the two image datasets. The Dice index increases as spatial overlap increases. We chose a threshold of 0.7 to reflect high spatial similarity between two connectivity maps (Zou et al., 2004, Zijdenbos et al., 1994). Functional connectivity measures were calculated across the whole brain and Dice indices were calculated to allow pairwise comparisons between methods for all significant connectivity clusters in the brain as well as for the global connectivity peak and for connectivity clusters other than the global connectivity peak (i.e. ‘distant’ peaks). Fig. 1 summarises the analysis method.

Fig. 1

Diagram representing the steps for the spatial and temporal comparison between non-invasive irritative zone localisation methods. Steps 1 to 3 show the acquisition and pre-processing steps using a single simultaneously acquired dataset. Step 4: Localisation of the irritative zone using spike-related, voltage-map EEG-fMRI and electrical source data. Step 5: Spatial maps were normalised to MNI space and the discrepancy between peaks was quantified at the lobar and sub-lobar levels. Step 6: Using each global maximum as seeds to assess the temporal correlation between BOLD timecourses and voxel-based functional connectivity analysis to investigate the spatial similarity between each functional network. BOLD: Blood-oxygen level-dependent;

Temporal correlation analysis

To assess the temporal relationship between the peak BOLD signal time course at each peak, a correlation analysis was performed between extracted time courses. R2 values were statistically compared with a significance level set at p < 0.05.

Contralateral functional connectivity seeds

We placed additional ‘control’ seeds to compare the temporal and functional connectivity patterns found within the irritative zone against those in regions outside this zone. Seeds were placed in the contralateral homologous region by mirroring their location across the midsagittal plane using Fslswapdim (FSL v5.0).

Statistical analysis

A Shapiro-Wilk test was used to test for normality and as the data were found not to be normally distributed, non-parametric Kuskall-Wallis and Mann-Whitney-U tests were used to compare the variables. For categorical variables, Fisher’s exact test was applied. The level of significance was set at p < 0.05. Data are presented as median with their respective interquartile range.

Results

Seventeen patients were scanned in this study. Two patients (patients 6 and 11) were excluded because of incomplete acquisition due to seizure inside of scanner (patient 6) and the presence of residual EEG hardware artefact after pre-processing (patient 11). Of the remaining patients, 11 (73%) had IEDs during the simultaneous EEG-fMRI acquisition. In total, 1126 IEDs were identified with a median of 78 IEDs (interquartile range 57 – 127).

Localisation of global maxima using ESI and EEG-fMRI methods

For patients without IEDs during the EEG-fMRI acquisition, only ESI and voltage-map EEG-fMRI maps were obtained. Using the traditional spike-related EEG-fMRI approach, the localisation of the global maxima was found in 8 out of the 11 patients with intra-scanner IEDs, whereas for voltage-map EEG-fMRI and ESI a significant BOLD global maxima and source of maximal strength was found in 12/15 and 15/15 patients, respectively. The number of spikes per patient and those patients with significant peaks for each method are detailed in Table S1.

Spatial localisation of the irritative zone and distance between global maxima

Localisation of the global maximum for spike-related, voltage-map and ESI analyses at the lobar and sub-lobar level and the Euclidean distances between peaks are shown in Table 2. For all methods maxima were co-localised in the same hemisphere. At lobar level, there was concordance between all three methods in the 8 patients in which all methods were feasible. Similarly, global maxima from the voltage-map analysis and ESI were in the same lobe in all 12 patients studied with the 2 methods. At the sub-lobar level, spatial concordance was significantly lower compared to the lobar level. Co-localisation within the same Brodmann Area was seen in only in 2/12 patients with both voltage-map analysis and ESI (p < 0.001, Fisher exact test) and in 3/8 (p < 0.03) of patients with both spike-related analysis and ESI and 3/8 (p < 0.03) patients with both spike-related and voltage-map analyses.

Table 2

MNI coordinates (x, y, z) of the maximum peak using Electrical Source Imaging (ESI), voltage map (VM) and spike-related (SR) EEG-fMRI. The Euclidean distances in mm between all peaks are shown on the right. For patients with no spikes inside of the scanner only Electrical Source Imaging and voltage-map EEG-fMRI peaks (*) and not considered in the displayed mean. L: Left, R: Right, AAL: Automated Anatomical Labelling atlas, BA: Brodmann Area.

#	Peak	Electrical Source Imaging	Voltage map EEG-fMRI	Spike-related EEG-fMRI	Distance to peak (mm)
					ESI VM	ESI SR	VM SR
P01	MNI	−33, −1, −34
	Lobe	L. Temporal	–	–	–	–	–
	AAL | BA	Fusiform | BA20
P02	MNI	−57, −22, −30
	Lobe	L. Temporal	–	–		–	–
	AAL | BA	Inferior temporal | BA20
P03	MNI	−38, 27, −10	−57, 9, 3	−51, 21, 9
	Lobe	L. Frontal	L. Frontal	L. Frontal	29.2	23.8	14.7
	AAL | BA	Inferior frontal | BA47	Rolandic operculum | BA44	Inferior frontal | BA45
P04	MNI	−33, 5, −33	−33, 9, −3	−27, 3, −12
	Lobe	L. Temporal	L. Insula	L. Temporal	30.3	21.9	12.4
	AAL | BA	Inferior temporal | BA38	Putamen | BA13	Amygdala | BA49
P05	MNI	−32, −42, 48	−42, −30, 51	−12, −57, 66
	Lobe	L. Parietal	L. Parietal	L. Parietal	30.8	15.9	43.1
	AAL | BA	Inferior parietal | BA7	Postcentral | BA1	Precuneus | BA7
P07	MNI	45, 42, −9	6, 51, 9	–
	Lobe	R. Frontal	R. Frontal		43.9*	–	–
	AAL | BA	Inferior frontal | BA47	Anterior cingulate | BA10
P08	MNI	51, 18, 18	51, 18, 15	51, 18, 27
	Lobe	R. Frontal	R. Frontal	R. Frontal	3	9	12
	AAL | BA	Inferior frontal | BA44	Inferior frontal | BA44	Inferior frontal | BA44
P09	MNI	45, −3, −33	27, −12, −33	–
	Lobe	R. Temporal	R Temporal		20.1*	–	–
	AAL | BA	Inferior temporal | BA20	Fusiform | BA36
P10	MNI	45, −55, −19	–	–
	Lobe	R. Temporal			–	–	–
	AAL | BA	Fusiform | BA37
P12	MNI	36, −3, −36	36, 3, −6	36, −3, −3
	Lobe	R. Temporal	R. Insula	R. Insula	30.6	33	6.7
	AAL | BA	Fusiform | BA20	Putamen | BA13	Putamen | BA49
P13	MNI	−39, 18, −33	−30, −9, −36	–
	Lobe	L. Temporal	L. Temporal		28.6*	–	–
	AAL | BA	Mid temporal pole | BA38	Fusiform | BA36
P14	MNI	−48, −51, 30	−51, −42, 57	−54, −60, 42	28.6	16.2	23.6
	Lobe	L. Parietal	L. Parietal	L. Parietal
	AAL | BA	Angular | BA39	Inferior parietal | BA40	Inferior parietal | BA39
P15	MNI	−30, 0, −27	−36, 0, −9	−36, 6, 0
	Lobe	L. Temporal	L. Insula	L. Insula	19	28.3	10.8
	AAL | BA	Amygdala | BA36	Insula | BA13	Insula | BA13
P16	MNI	36, −18, −18	27, −18, −18	30, −18, −18
	Lobe	R. Temporal	R. Temporal	R. Temporal	26.3	25.8	3
	AAL | BA	Hippocampus | BA54	Hippocampus | BA54	Hippocampus | BA54
P17	MNI	30, −1, −19	27, −9, −24	–
	Lobe	R. Temporal	R. Temporal		9.9*	–	–
	AAL | BA	Amygdala | BA53	Hippocampus | BA54
				Mean	24.7	21.7	15.8
				SD	9.6	7.7	12.6

The mean Euclidean distance between pairs of global maxima from each method (Fig. 2) was 20.7 ± 10 mm. The distances for all methods were significantly greater than zero (p = 0.0078, Wilcoxon signed-rank test) confirming that the peaks were spatially discordant. Voltage map and spike-related EEG-fMRI showed the smallest inter-peak distance (median 12.2 mm, interquartile range 7.7 – 21.4 mm) whereas the largest distance was between source imaging and voltage map EEG-fMRI (median 28.9, interquartile range 20.8 – 30.5 mm) but the distance between peaks was not statistically different between methods (p < 0.09, Kruskal-Wallis test).

Fig. 2

Graph showing the median (symbol) and interquartile range (bars) of all pairwise Euclidean distances between global maxima. ESI: Electrical Source Imaging, SR: Spike-related and VM: voltage map EEG-fMRI. * p < 0.05, ** p ≤ 0.01 One-sample Wilcoxon test.

Functional connectivity between global maxima

The Pearson correlation between peak BOLD time courses was calculated and shown in Fig. 3a. BOLD time courses were significantly correlated, with the median R2 value across all methods being 0.35 (interquartile range 0.14–0.6; p = 0.0078, Wilcoxon test). The highest correlation was between spike-related and voltage-map EEG-fMRI (0.49; interquartile range: 0.22–0.74) and the lowest was between voltage-map EEG-fMRI and ESI (0.26; interquartile range 0.1–0.4). However the differences between R2 values did not differ significantly between methods (p = 0.42, Kruskal-Wallis test).

Fig. 3

Functional connectivity analysis a) Bar graph showing the median Dice similarity value and interquartile range (bars) for ipsilateral peaks at each level. b) Violin plot showing the Dice similarity value comparison between ipsilateral and contralateral peak seeds. c) Temporal correlation between functional connectivity fMRI peaks in the ipsilateral and in their respective ROIs on the contralateral hemisphere (darker colours). Triangles and bars representing the mean and standard deviation, respectively. We found low temporal correlation across all ipsilateral BOLD peaks with no significant differences between techniques. The temporal correlation between ipsilateral peaks could not explain the majority of the variance observed between all methods. Contralateral ROIs, used as a control, showed significantly lower correlation values than their respective ipsilateral seed. (* p < 0.05, ** p < 0.01; *** p < 0.001, One-sample Wilcoxon test, ESI: Electrical Source Imaging, SR: Spike-related EEG-fMRI, VM: Voltage map EEG-fMRI, ipsilat.: ipsilateral seed, contralat.: contralateral seed).

At each of the three levels for which conducted a Dice similarity analysis, we found no significant differences in the spatial distributions of connectivity maps with global maxima from each of the analysis methods (Whole brain: p = 0.32; Seeding peak cluster: p = 0.58; Distant connectivity clusters: p = 0.33, Kruskal-Wallis test).

No pairwise comparisons exceeded the Dice index threshold for good spatial agreement at any of the three levels analysed here. The mean Dice index was the highest for pairwise comparisons of peak connectivity clusters (median 0.65, interquartile range 0.27 – 0.75). This value was significantly higher than that for clusters distant to the peak (median 0.09, interquartile range 0.08 – 0.47, p < 0.02) (See Table S2). For the whole brain, the mean Dice index across all functional maps was 0.5 ± 0.26. (Fig. 3b). For each pairwise comparison, both the temporal correlation between BOLD time courses of the global maxima for each method and the Dice similarity index were inversely related with the Euclidean distance between peaks (Fig. S1).

Contralateral seeds

We compared the connectivity patterns and temporal correlation of ipsilateral peaks with seeds located outside the irritative zone, placed at the homologous region of the contralateral hemisphere. We observed a significantly lower pairwise BOLD time course R2 values of ESI/VM (p < 0.04), ESI/SR (p < 0.03) and VM/SR (p < 0.007) in contralateral versus ipsilateral seeds (See Fig. 3). Similarly, Dice similarity analysis of contralateral seeds showed significantly lower spatial concordance between functional connectivity maps at the three levels of analysis used in this investigation (whole brain: p < 0.001; peak cluster: p < 0.001; distant cluster: p < 0.03) (See Table S3).

Illustrative case of the relationship between inter-peak distances and functional connectivity measures

Patient 14 (left parietal lobe epilepsy) (Fig. 4) had a left parietal glioma and ipsilateral parietal lobe IEDs. Although all non-invasive peaks were co-localised within the left parietal lobe adjacent to the epileptogenic lesion, their peaks were spatially distributed within the lobe with a mean Euclidean distance of 22.8 ± 6.2 mm between peaks. Functional connectivity maps showed poor overall spatial concordance with a mean Dice index of 0.6 ± 0.1 for pairwise comparisons between methods. The mean R2 value for correlations between BOLD time courses at global maxima identified with each method explained <30% of the variance between BOLD time courses.

Fig. 4

Patient 14 showed an epileptogenic lesion in the left temporal lobe. Interictal spikes were confined to the left posterior-quadrant region. a) Electrical Source Imaging maxima (ESI), Spike-related (ST) and Voltage-map EEG-fMRI BOLD co-localise its respective global maxima in the left parietal lobe, all adjacent to the structural lesion. However, peaks were all distant to each other (b) with a mean Euclidean distance of 22.8 ± 6.24 mm. (c) Seed-voxel functional connectivity maps (p > 0.05 FDR corr.) derived from seed ROIs (black dots) of each respective non-invasive localisation method. Clusters containing the peak are highlighted by the black arrow. The patient showed low spatial correlation between networks (d) with a mean DICE index of 0.6 ± 0.11 between peak clusters and (f) low correlation between BOLD time courses with a mean R2 value of 0.29 ± 0.12 between all methods.

Discussion

We investigated the relationship between non-invasive localisation of the Irritative Zone using ESI and spike-related and voltage-map EEG-fMRI. Although global maxima for the three methods were in the same lobe, there were significant differences in localisation. The time courses of the resting state BOLD signal at global maxima identified by each method showed modest temporal correlation which did not account for most of the variance in the time courses. The functional connectivity of each of the maxima did not show good spatial agreement between methods. Spatial and temporal connectivity measures between peaks inside the irritative zone were higher than for seeds placed in the contralateral homologous regions. Taken together, these findings suggest that ESI and voltage-map and spike-related EEG-fMRI have systematic differences in the co-localisation of their respective maxima each of which reveals a distinct aspect of the underlying network associated with the generation of IEDs.

Spatial agreement between the non-invasive localisation methods facilitates the interpretation and clinical yield of the presurgical epilepsy workup. In this investigation we analysed the spatial relationship between two different EEG-based modalities, ESI and EEG-fMRI. Each modality measures fundamentally different aspects of IEDs, with ESI providing an estimate of the electrical source of IEDs recorded on the scalp and EEG-fMRI reflecting the changes in the proportion of oxygenated and deoxygenated blood associated with the generation of IEDs. Consequently, we did not expect exact concordance between ESI and EEG-fMRI findings (Nunez and Silberstein, 2000, Bagshaw et al., 2006). Nonetheless, our spike-related data were simultaneously acquired so that the observations relate to the same neuronal phenomenon. We observed a significant difference in the location of global maxima from ESI and the EEG-fMRI methods. This finding is in agreement with previous studies assessing the spatial concordance between source imaging and EEG-fMRI which have reported differences of 15 to 60 mm (Lemieux et al., 2001, Bagshaw et al., 2006, Benar et al., 2006, Vulliemoz et al., 2009, Elshoff et al., 2012, Centeno et al., 2017). To our knowledge, differences in the localisation of the global maxima of spike-related and voltage-map EEG fMRI has not been previously compared. The median 12 mm inter-peak distance that we observed cannot simply be explained by the spatial resolution of T2* BOLD contrast, which is in the range of 4–5 mm (Kim et al., 2004, Uğurbil et al., 2003).

We postulate that the differences in location of global maxima is related to the differences in temporal resolution of each of the techniques. This is most evident when comparing ESI to fMRI-based techniques. The high temporal resolution of ESI allows differentiation between regions responsible for the generation of IEDs, which are related to the rising phase of the discharge, from those involved in propagation, which are related to later timepoints of the IED (Birot et al., 2014, Lantz et al., 2003). The EEG-fMRI approaches also differ in temporal resolution. The spike-related EEG-fMRI approach reveals the neural generators of the fast component of the discharge as the regressor of interest follows the time course of the peak of the IED. In contrast, for voltage-map EEG-fMRI the regressor of interest also reflects the topography of components such as the slow wave and focal EEG slowing that may show similar scalp voltage topography as the spike (Grouiller et al., 2011). These low frequency events introduce periods of high correlation that sensitise the voltage map EEG-fMRI approach to a different time scale to both ESI and spike-related EEG-fMRI.

The correlation between BOLD time courses could only explain 38% of the variance between resting state BOLD time courses centres around the global maxima from ESI and the two EEG-fMRI anaylsis methods. This supports the hypothesis that the methods are identify different subsets of the temporally-correlated network associated with IEDs. Our analysis of spatial concordance leads to the same conclusion. Even though the temporal correlation and spatial overlap was small, it was significantly higher when compared against the respective functional connectivity measures using seeds placed outside the irritative zone in the contralateral homologous region. Overall, these observations provide evidence that each peak reveals a different aspect of a locally distributed network involved in the generation of interictal spikes. A combination of information from the three methods could provide complementary information that improves the delineation of the extent of the irritative zone.

Our findings are in agreement with recent neurophysiological observations on the spatial and temporal network dynamics during the generation of interictal discharges. IEDs have traditionally been ascribed to discrete cortical generators of a spontaneous burst of action potentials resulting in a large depolarisation followed by long hyperpolarisation (Matsumoto and Marsan, 1964). However, recent investigations using simultaneous single unit and local field potential recordings have demonstrated that discharges are not a discrete paroxysmal event of hypersynchronous excitatory activity but rather arise from heterogeneous neuronal firing patterns of multiple distinct neuronal subsets of a locally distributed neuronal network (Keller et al., 2010, Feldt Muldoon et al., 2013). Additionally, intracranial observations have revealed that the neurological cascade of events that gives rise to scalp-recorded interictal discharges spans a longer timescale, with periods of increased neural synchronisation occurring during the fast component of the discharge and modulation of firing patterns prior to and after the spike peak (Keller et al., 2010).

Limitations

Our sample size is small. Of the 15 patients analysed in this study, only 8 (53%) patients had IEDs inside the scanner. This proportion is similar to previous investigations (Salek-Haddadi et al., 2006) in which nearly half of the SR EEG-fMRI studies were inconclusive due to the lack of IEDs. VM EEG-fMRI was included in the study to address this limitation of SR EEG-fMRI (Grouiller et al., 2011).

The purpose of this investigation was not to assess which of the methods is better for the non-invasive localisation of the epileptogenic region. Further studies with post-operative outcome data will be necessary to assess whether combining information from different EEG-based methods confers benefits in presurgical localisation and surgical planning. Studies using the location of each peak as a target for electrode implantation could also elucidate the links between the different aspects of the network involved in spike generation and seizure onset.

As the localisation of the irritative zone has been shown to evolve and fluctuate throughout time (Rosenow, 2001), we performed our analyses on simultaneously acquired EEG and fMRI data. For ESI, we used an inverse solution that incorporates biophysical models of electrical propagation and a patient-specific head model. This has been previously shown to identify neural correlates of interictal activity with similar accuracy to more computationally intensive approaches such as the boundary element model (Birot et al., 2014). EEGs recorded outside the scanner employed 21-channels, which may capture the detailed spatial representation of the scalp electrical distribution of IEDs (Lantz et al., 2003). Future studies could investigate the effect of spatial sampling with a greater number of electrodes on voltage-map EEG-fMRI.

Conclusions

We provide evidence that ESI and spike-related and voltage-map EEG-fMRI provide information on different subsets of the network underlying interictal spikes. Spatial and temporal mismatch between identified parts of the network are likely to relate to differences in temporal resolution between the techniques.

CRediT authorship contribution statement

Javier Urriola: Conceptualization, Formal analysis, Investigation, Methodology, Project administration, Writing - original draft. Steffen Bollmann: Methodology, Writing - review & editing. Fred Tremayne: Investigation, Methodology, Writing - review & editing. Hana Burianová: Supervision, Writing - review & editing. Lars Marstaller: Supervision, Writing - review & editing. David Reutens: Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing - review & editing.