Connectivity strength, time lag structure and the epilepsy network in resting-state fMRI.

The relationship between the epilepsy network, intrinsic brain networks and hypersynchrony in epilepsy remains incompletely understood. To converge upon a synthesized understanding of these features, we studied two elements of functional connectivity in epilepsy: correlation and time lag structure using resting state fMRI data from both SEEG-defined epileptic brain regions and whole-brain fMRI analysis. Functional connectivity (FC) was analyzed in 15 patients with epilepsy and 36 controls. Correlation strength and time lag were selected to investigate the magnitude of and temporal interdependency across brain regions. Zone-based analysis was carried out investigating directed correlation strength and time lag between both SEEG-defined nodes of the epilepsy network and between the epileptogenic zone and all other brain regions. Findings were compared between patients and controls and against a functional atlas. FC analysis on the nodal and whole brain levels identifies consistent patterns of altered correlation strength and altered time lag architecture in epilepsy patients compared to controls. These patterns include 1) broadly distributed increased strength of correlation between the seizure onset node and the remainder of the brain, 2) decreased time lag within the seizure onset node, and 3) globally increased time lag throughout all regions of the brain not involved in seizure onset or propagation. Comparing the topographic distribution of findings against a functional atlas, all resting state networks were involved to a variable degree. These local and whole brain findings presented here lead us to propose the network steal hypothesis as a possible mechanistic explanation for the non-seizure clinical manifestations of epilepsy.

epileptogenic zone

false discovery rate

full scale intelligence quotient

propagation zone

standard deviation

stereo-electroencephalography

resting state functional MRI

resting state network

Introduction

Localization related epilepsies are complex disorders characterized by local and remote brain alterations. They are widely accepted as network disorders (reviews in Stefan and Lopes da Silva, 2013; Bartolomei et al., 2017) and are known to be associated with hypersynchrony at their pathologic foundation (Bartolomei et al., 2001; Bartolomei et al., 2004; Schevon et al., 2007; Warren et al., 2010). The precise relationship between the epilepsy network and the abnormal synchrony in epilepsy remains unclear however.

Due to the temporal resolution advantages of electrophysiological techniques, synchrony in epilepsy has historically been explored via analysis of clinical intracranial EEG data or single unit recordings. While intracranial EEG represents the gold standard for classifying regions in terms of their role in epilepsy, fMRI offers the possibility of direct comparisons with controls. The specificity afforded by invasive electrophysiological analysis has been used alongside fMRI in the investigation of focal network changes in epilepsy relative to healthy controls both in terms of proximity with voxelwise fMRI connectivity contrasts (Stufflebeam et al., 2011), and atlas-based regions of interest (ROIs) (Weaver et al., 2013). However, as of yet, investigation of connectivity in fMRI using ROIs derived directly from intracranial EEG have been compared only to the ‘internal’ control of ‘non-involved’ regions within patients (Bettus et al., 2011). Furthermore, the limitations of intracranial EEG's limited spatial coverage cannot take into account the varying natural time lags between remote brain regions which have been shown to be particularly prominent in some areas, whereas fMRI permits whole-brain comparisons. This may be particularly relevant given recent work by Mitra, and colleagues (2014) and their description of the physiologic lag architecture of the healthy human brain through their investigation of the temporal features within both rsfMRI activity and its modulation in task states. Importantly, the patterns in the latency structure of the brain identified by this work were highly reproducible over several large cohorts. This suggests an underlying dynamic framework that is conserved across individuals. Focal modulation by various tasks supports the functional relevance of this framework and implies that neural processes, rather than hemodynamic lag, are the major determinants of this latency structure. Furthermore, prominent loci identified by this approach correspond to regions with behavior-regulating functions in task states and that these constitute key nodes in resting state networks (RSNs), such as the default mode, salience and frontoparietal networks, suggests that this lag architecture represents a framework for integration of the intrinsic connectivity structure represented by RSNs.

Thus, the investigation of latency analysis in rsfMRI data could speak to both focal alterations within components of the epilepsy network as well as epilepsy's impact on underlying, functionally-relevant systems at the whole brain level. Specifically, using SEEG-defined nodes of the epilepsy network (including cortical and sub-cortical brain regions demonstrating both initial seizure onset (EZ) and seizure propagation activity (PZ) as defined by clinical and quantitative SEEG interpretation (Bartolomei et al., 2017; Guye et al., 2007)), we investigate patterns of correlation strength and time lag relative to these network nodes in a cohort of 15 epileptic patients compared to 36 healthy control subjects using rsfMRI data with the hypothesis being that enhanced directed interdependencies exist between nodes of the epilepsy network.

Materials and methods

Subjects

Fifteen patients (mean ± SD age: 33.4 ± 10.4 years, 8 females) with intractable epilepsy managed at the Epileptology department of La Timone Hospital in Marseille, France, were included in this study (Table 1). All patients underwent a pre-surgical evaluation including history and neurological examination, neuropsychological testing, MRI and scalp EEG. Following multi-disciplinary review of this non-invasive evaluation, all patients were deemed candidates for SEEG-guided exploration of a surgical hypothesis due to insufficient anatomo-electroclinical concordance of non-invasively acquired data (Isnard et al., 2018). A summary of electrode locations and EZ/PZ localization results following SEEG evaluation are presented in Table 2.

Table 1

Clinical features of the patients. L: left, R: right, B: bilateral, FSIQ: full-scale IQ.

ID	Gender	Age at onset (y)	Duration (y)	Side	Seizure Onset Localization	Etiology	MRI Findings	Engel class	Follow-up (y)	FSIQ
1	F	8	14	R	Frontal	Post-infectious	R temporo-occipital ulegyria	IV	4	54
2	F	15	7.9	B	Temporal	Cryptogenic	Normal	III	3	83
3	F	9	13	L	Occipital	Cryptogenic	Normal	IV	2	83
4	M	7	35	L	Operculo- insular	Cortical Dysplasia	FCD, L insula	III	3	83
5	F	21	15.5	B	Temporal	Hypothalamic Hamartoma	Hypothalamic Hamartoma	NA	NA	116
6	F	4	19	L	Frontal	Cortical Dysplasia	FCD, L frontal	Ia	1.5	81
7	M	0.2	56	R	Fronto-temporal	Cortical Dysplasia	FCD, R frontal	Ia	4	73
8	M	30	3	R	Fronto-temporal	Cortical Dysplasia	FCD, R amygdala	IV	4	100
9	M	18	12	R	Fronto-temporal	Arteriovenous Malformation	R frontal encephalomalacia	Ia	4	None
10	F	17	10	R	Temporal	Hippocampal Sclerosis	R hippocampal sclerosis	II	4	87
11	F	18	24	L	Frontal	Post-traumatic	L temporal encephalomalacia	NA	NA	85
12	F	23	26.8	R	Temporal	Post-traumatic	R fronto-temporal encephalomalacia	Ia	3	86
13	M	7	30	R	Temporal	Cavernous Malformation	R temporal cavernoma	Ia	2	110
14	F	11	14	B	Occipital	Cryptogenic	Dilation of R occipital horn	NA	NA	56
15	M	17	14	L	Frontal	Peri-natal injury	L frontal porencephalic cyst	Ia	2	None

Table 2

Location of implanted SEEG electrodes and localization of epileptogenic and propagation zones. EZ: epileptogenic zone, PZ: propagation zone.

ID	Location	EZ	PZ
1	Right orbitofrontal	1, 2	1–15
	Left orbitofrontal		1–15
	Right premotor		1–8
	Right anterior cingulate		1–12
	Right mid cingulate
	Right frontal operculum
	Right basal temporal
	Right temporopolar	1–4
2	Right parietal operculum	1–3
	Right anterior hippocampus	8–10	2–5
	Right posterior hippocampus	8–10	2–4
	Right Heschl's gyrus	3, 4, 7, 8
	Right temporopolar	1, 2	8–10
	Left amygdala	1–4	6, 7, 9, 10
3	Left parieto-occipital
	Left posterior parietal
	Left cuneus		2–4
	Left posterior fusiform gyrus
	Left isthmus of the cingulate gyrus
	Left anterior calcarine fissure		5–8
	Left lingual gyrus	4–6
	Left posterior parahippocampus
	Left occipitotemporal
	Left anterior hippocampus
	Left basal temporal
	Right posterior fusiform gyrus
	Right occipitotemporal
4	Left parietal operculum	2–5	6–8
	Left frontal operculum		1–3
	Left superior temporal gyrus
	Left temporopolar
	Left posterior hippocampus
	Right isthmus of the cingulate gyrus
5	Left anterior hippocampus
	Left posterior hippocampus	1–4	2, 3
	Left temporopolar	2, 3
	Left basal temporal		1, 2
	Left posterior parahippocampus		1, 2
	Left heschl's gyrus
	Hypothalamic hamartoma
	Right orbitofrontal		1–4
	Right basal temporal
6	Left premotor
	Left supplementary motor	1–4
	Left mid cingulate	1–3
	Left occipitotemporal
	Left central sulcus		6–9
	Left posterior cingulate		10–13
	Left anterior parietal
	Left parietal operculum
	Left parieto-occipital
7	Right orbitofrontal		1–2
	Right insula
	Right premotor		1–8
	Right anterior cingulate		1–12
	Right pars triangularis		1–8
	Right basal frontal		1–10
	Right amygdala	1–4
	Right temporopolar	1–4
	Right superior temporal gyrus
	Right anterior hippocampus		1–3
	Right basal temporal		1–3
	Right mid cingulate
	Left orbitofrontal
8	Right temporopolar
	Right amygdala	1–3
	Right posterior parahippocampus
	Right anterior hippocampus	1–3
	Right basal frontal		1–3
	Right superior temporal gyrus		1–3
	Right anterior cingulate
	Right dorsal frontal
	Right mid cingulate
	Right orbitofrontal
9	Right frontopolar	1–10
	Right basal frontal	1–3
	Right orbitofrontal		1–5
	Right frontal operculum
	Right premotor		1–8
	Right supplementary motor		1–8
	Right temporopolar
	Right left mid cingulate		5–10
	Left premotor
	Left supplementary motor
10	Right pars triangularis
	Right orbitofrontal
	Right temporopolar	1–3
	Right frontal operculum		1–5
	Right amygdala	1–3
	Right basal temporal	1–4
	Right insula		1–5
	Right parietal operculum
	Right Heschl's gyrus		1, 2
	Right basal frontal	1–4
	Right posterior hippocampus
11	Left anterior cingulate	1, 2, 5–11
	Left mid cingulate	11, 12
	Left posterior cingulate
	Left temporopolar	1–3
	Left amygdala		1–3
	Left anterior hippocampus		5, 6, 10–12
	Left basal temporal
	Left posterior parahippocampus
	Left orbitofrontal	1,2
	Left posterior hippocampus	9, 10	2, 3
	Left frontal operculum	7–10	5, 6
	Left premotor	7–9
	Right orbitofrontal		6–9
12	Right orbitofrontal
	Right frontal operculum
	Right parietal operculum
	Right Heschl's gyrus	4, 5	1–3, 6–8
	Right temporopolar	1–4	5–8
	Right amygdala
	Right anterior hippocampus		1, 2
	Right basal temporal		9–11
	Right posterior parahippocampus
13	Right basal temporal		1–3, 5–9
	Right amygdala		1–4
	Right parietal operculum
	Right anterior hippocampus		1–3
	Right parieto-occipital
	Right posterior hippocampus
	Right occipitotemporal	1–3
14	Right amygdala
	Right anterior hippocampus
	Right occipitotemporal		4–8
	Right lingual gyrus	1–3, 10–12
	Right cuneus	1–10
	Right anterior calcarine sulcus
	Right parietal operculum
	Right posterior isthmus of the cingulate gyrus
	Left lingual gyrus	1–3
15	Left Anterior cingulate
	Left mid cingulate		6–10
	Left premotor		4–9
	Left supplementary motor		1, 2
	Left central sulcus	5–10
	Left frontal operculum
	Left orbitofrontal
	Left parieto-occipital
	Left Anterior parietal
	Left parietal operculum	3, 4	5
	Right orbitofrontal

Neuropsychological evaluations were reviewed for consistency across patients. Local records were available for 13 of the 15 patients. Given significant variability in the tests administered to each patient, only the FSIQ was available for each of the 13 patients (mean ± SD FSIQ: 84.38 ± 17.76 range: 116–54).

Comparisons were performed against 36 age and gender matched healthy subjects (mean ± SD age: 31.2 ± 8.5 years, 19 females) with no history of neurological or psychiatric disease. This study was approved by the local ethics committee and informed written consent was obtained from all patients and healthy participants.

Image acquisition

MR images were acquired on a 3T Siemens Magnetom Verio MR scanner (Siemens Healthineers, Erlangen, Germany). The protocol comprised a rapid gradient echo anatomical T1-weighted image (3D-MPRAGE, TR = 1900 ms, TE = 2.19 ms, TI = 900 ms, voxel size = 1 × 1 × 1mm3) and a functional MRI using a BOLD-sensitized single-shot EPI T2*-weighted sequence (TR = 3.6 s, TE = 27 ms, FA = 90°, voxel size = 2 × 2 × 2.5mm3, 50 slices, 350 volumes) during which the subjects were awake with their eyes closed. This resulted in 20 min-long fMRI time series for analysis.

Structural image preprocessing

T1 weighted images were processed with Freesurfer (v5.3, https://surfer.nmr.mgh.harvard.edu/), which includes cortical surface reconstruction and registration to a common surface space using a non-rigid curvature alignment procedure (Fischl et al., 1999). Subcortical structures (accumbens nucleus, amygdala, caudate, hippocampus, pallidum, putamen, thalamus) were segmented using FIRST (FSL v5.0, https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/), then surface models were obtained with SPHARM-PDM allowing for subcortical mesh parameterization and registration to the corresponding surface template (Besson et al., 2017a.b; Fischl et al., 1999).

Functional image preprocessing

Functional MRI data were corrected for motion and slice timing using SPM12 (http://www.fil.ion.ucl.ac.uk/spm/software/spm12/) and aligned with corresponding T1 image using boundary-based registration. Further preprocessing steps were handled using the NIAK toolbox (https://www.nitrc.org/projects/niak/) and included removal of time series offset and linear trend, band pass filtering to retain frequencies between 0.01 and 0.08 Hz, linear regression using the six parameters of the motion correction, average ventricular signal, average white matter signal, average whole-brain signal, and their temporal derivatives (18 regressors in total) (Murphy and Fox, 2017; BT Thomas Yeo et al., 2011). Volume scrubbing was performed according to framewise displacement and DVARS to avoid spurious connectivity (Jonathan D. Power et al., 2012). After volume scrubbing, the remaining time did not differ significantly (p > 0.7, t-test with unequal variances) between healthy controls (16.98 min ± 2.64 min) and patients (16.66 min ± 4.48 min).

Finally, corrected fMRI volumes were projected onto corresponding cortical and subcortical surfaces. Cortical surface projection was performed by trilinear interpolation of fMRI frames along the medial cortical surface (mid distance between inner and outer cortical surfaces) discarding the midbrain portion of the cortical surface (Freesurfer's medial_wall label). A 10 mm FWHM surface Gaussian kernel blur was then applied in native space and the cortical surface frames were registered to the common cortical surface (Freesurfer's fsaverage5). The fMRI signal for each vertex of the subcortical surface was determined by a spatial averaging, weighted by a FWHM 6 mm Gaussian kernel centered at the vertex over the voxel within the subcortical region only. This resulted in measurement primarily of the superficial subcortical signal at each vertex. Subcortical signals were then registered to corresponding surface templates. Cortical and subcortical surface projection is summarized in Fig. 1.

Fig. 1

Cortical and subcortical surface projection. Top row demonstrates the mid-cortical surface alignment on the T1 and average EPI with a magnified panel depicting a close-up representation of the right hemispheric cortical surface. Bottom row demonstrates the thalamic surface projection on the T1 and average EPI with a magnified panel depicting the projection filter centered on the green dot.

Electrode definition & localization

Following clinical review and quantitative analysis of the SEEG data, all implanted SEEG electrodes were defined as belonging to one of three regions (Isnard et al., 2018): 1) epileptogenic zone (EZSEEG, defined as those electrodes demonstrating initial seizure activity), 2) propagation zone (PZSEEG, defined as those electrodes demonstrating subsequent seizure activity following initial seizure onset) and, 3) brain regions not primarily involved in seizure activity (NPISEEG, defined as those electrodes which were not primarily involved in initial seizure onset or initial seizure evolution) (Bartolomei et al., 2017; Guye et al., 2007). All patients underwent either a post-implantation volumetric CT scan or T1 weighted MRI for co-registration with the pre-implantation T1 weighted image for SEEG electrode localization.

Regions of interest

Three regions of interest were defined for connectivity and time lag analysis based on the clinical SEEG interpretation as previously described (Besson et al., 2017a.b). Briefly, the nodes of the epilepsy network EZ and PZ were defined to include all cortical and sub-cortical vertices within a 10 mm radius around an electrode contact found to be within the EZSEEG or PZSEEG, respectively (Supporting Fig. 1). Ten millimeters was selected due to the mapped spatial sensitivity of electrode contacts in clinical practice (Haglund et al., 1994; Rolston and Chang, 2017; Sanai et al., 2008) and experimental conditions (He et al., 2008; Roland, 2002; Roland and Zilles, 1998). Vertices not assigned to either EZ or PZ were grouped to form a third zone termed ‘not primarily involved’ (NPI). Therefore, NPI brain regions included both SEEG-sampled brain regions and all non-SEEG sampled brain regions. A known limitation of SEEG is the balance between electrode density and clinical feasibility. In our patient cohort, SEEG electrodes were placed based on clinical suspicion for epilepsy network involvement but it is possible that additional regions involved in the epilepsy network were not sampled by SEEG and therefore were included in NPI. Distances from electrode leads were mapped on native surfaces, registered to the common surface space and thresholded at 10 mm to define zone membership. This procedure ensured anatomical correspondence between patients and controls, allowing zones defined in each patient to be mapped in healthy controls for further comparison.

Correlation and time lag

Correlation and time lag between two time courses, T1 and T2, were calculated using the lagged Pearson's coefficient of correlation, ρ:Where τ is the time lag, μ1and μ2 are time average of T1 and T2 respectively, σ1and σ2 their standard deviations, and N the number of time points. The connectivity strength, C, between T1 and T2 was defined so as to maximize the absolute value of ρ over lags:

The time lag structure L between T1 and T2 was set as the absolute value of the time lag that maximizes C. Said another way, L is a value ≥ 0 with 0 representing maximal strength of correlation between two signals at time lag 0. L > 0 indicates that the two signals correlate maximally if one signal is shifted in time (regardless of direction, we are interested primarily in interdependency rather than causality). Thus, the higher L, the less the two signals fluctuate synchronously:

Using a similar approach, one could define directed correlation and time lag from a region of interest R1 containing n1 vertices, to another region of interest R2 composed of n2 vertices. This way, the correlation strength between these two regions is obtained by computing the correlation between the average time course of all vertices in R1 and each vertex in R2:Where is the average time course of R1 and V is the time course of the vertex v in R2. Similarly, oriented time lag was obtained:

Therefore, using this approach, directed, inter-regional correlation strength and time lag can be calculated whereas (C(R1, R2) ≠ C(R2, R1) and L(R1, R2) ≠ L(R2, R1)). These directed relationships represent the tendency of all the vertices comprising R2 to be correlated or lagged with the average signal of R1, respectively. Self-correlation and lag may be obtained if R2 = R1 and indicate the tendency of vertices in R1 to exhibit the same fluctuation patterns as the regional average. Similarly, correlation from a region to a single vertex may be computed if n2 = 1.

Epilepsy network node-based analysis

Group comparison

rsfMRI correlation strength and time lag were computed within and across the SEEG-defined nodes mapped into the common surface space. For each patient, a total of 18 directed relationships were obtained (Eqs. (1) and (2)). Nine were directed correlation relationships and nine were directed lag relationships. Among these 9 directed relationships, 3 were intra-regional (EZ/EZ, PZ/PZ, NPI/NPI) and 6 were inter-regional (EZ/PZ, EZ/NPI, PZ/EZ, PZ/NPI, NPI/EZ, NPI/PZ). Similarly, directed correlations and lags were extracted from each healthy control using the same regions of interest so that for each patient, a collection of 36 correlation and lag matrices were obtained for comparison. All matrices were then corrected for age and gender using linear regression.

Z-scores were calculated in patients with respect to the healthy subjects whereas z-scores were calculated for healthy controls using a leave-one-out procedure. Group comparison was handled using a two-tail block bootstrap test under the null hypothesis of equal means and unequal variances, with one million random samplings, accounting for repeated measurements from healthy controls. P-values were corrected using FDR. The block bootstrap test was selected in order to allow us to deal with repeated measures, as each healthy subject was sampled several times therefore violating the independence assumption, and it does not rely on any assumption about the distribution of the data. Therefore, it is well-suited to assess mean difference across two distributions, each composed of several different patient-dependent distributions.

Clinical correlation analysis

Relationships between FSIQ and directed regional correlation and time lag in patients were assessed using a bootstrap procedure employing 300,000 iterations in order to robustly construct the null distribution of the correlations and determine their significance. P-values were corrected for multiple comparisons using FDR.

EZ-seeded whole brain analysis

Calculation of correlation strength and time lag when seeding from EZ

Correlation and time lag relative to EZ were computed at every vertex of the cortical and subcortical surfaces. For a vertex v, C(EZ, v) and S(EZ, v) were obtained according to Eqs. (1) and 2 for each patient and every healthy control. This resulted in correlation strength and time lag relative to EZ mapped over the entire cortical and subcortical surfaces. These maps were corrected for age and gender, and z-scored with respect to healthy controls for patients, or using a leave-one-out approach for controls. Significant group differences were investigated using a two-tail block bootstrap test employing 300,000 random samplings and results were reported for p < 0.05, FDR corrected and cluster size larger than 3 vertices.

Neuropsychological correlation analysis

Relationships between FSIQ and correlation strength or time lag were investigated with correlation analysis using a bootstrap approach employing 300,000 iterations to construct null distributions of correlations with IQ to determine the significance at each vertex. P-values were FDR corrected.

Comparison against a functional atlas

The topographic distribution of whole-brain group findings and clinical correlation findings were compared against a previously published functional atlas of the cortex(B. T. Yeo et al., 2011). We defined the weighted proportion of network involvement (NI) to quantify the proportion of significant results within each network, weighted by each vertex's individual t-value and corrected for each network's area. This provided the information of “which functional networks are most affected” and was calculated as follows:Where v is a vertex with a significant corrected p-value belonging to the functional network F, AF is the area of F. The denominator is the sum of t-values of all significant vertices divided by the total atlas area and is used for normalization such that . The network originally defined as ‘limbic’ by Yeo et al., was excluded from network inclusion as well as the total surface area due to its poor signal to noise ratio and poor reproducibility (Gordon et al., 2016; Yeo et al., 2014).

Results

Epilepsy network node-based correlation analysis

No significant directed correlation relationships were identified in epilepsy patients compared to controls on node-based analysis. Similarly, correlation analysis did not identify any significant relationship between correlation strength and pre-surgical FSIQ.

Epilepsy network node-based lag analysis

Decreased time lag was identified within EZ in patients with epilepsy compared to controls (p = 0.041, FDR corrected). Significantly increased time lag was also identified within all non-epilepsy network brain regions (all brain regions not included in EZ or PZ, p = 0.0002, FDR corrected) as well as between all non-epilepsy network brain regions and PZ (p = 0.02, FDR corrected). While not surviving correction for multiple comparisons, a trend towards increased lag was also identified between EZ and all other non-epilepsy network brain regions as well as between PZ and all other non-epilepsy network brain regions (Fig. 2). Clinical correlation analysis did not identify a significant relationship between time lag and pre-surgical FSIQ.

Fig. 2

Results from epilepsy network node-based lag analysis. Connection widths reflect relative nodal synchronization and are proportional to t-values (indicated adjacent to each connection; increased lag corresponds to decreased synchronization while decreased lag corresponds to increased synchronization). Significant findings are indicated by an asterisk (p < 0.005, FDR corrected). .

EZ-seeded whole brain correlation analysis

When evaluating whole brain correlation relative to EZ at the group level, multiple foci of increased connectivity and no foci of decreased connectivity were identified in patients compared to controls (Fig. 3). Foci of increased correlation strength were distributed over both the cortical and subcortical surfaces. Table 3 demonstrates the weighted proportion of network involvement by comparing the foci of increased connectivity (weighted by their t-value) against a functional atlas(B. T. Yeo et al., 2011). All six functional networks (excluding Yeo's ‘limbic’ network) demonstrated increased strength of correlation with EZ in patients compared to healthy controls. The somatomotor and salience networks together accounted for over 75% of the findings of increased strength of correlation when seeding EZ. All other networks including the visual, default mode, dorsal attention and frontoparietal control networks also demonstrated findings of increased strength of correlation to a lesser degree.

Fig. 3

Results from whole brain functional connectivity analysis. Top panel depicts cortical regions of increased connectivity when selecting the seizure onset node, EZ, as the seed for whole brain functional connectivity analysis. Comparison was made at the group level across all patients compared to all control subjects. Colorbar demonstrates significant t-values for both the cortical and subcortical surfaces. Bottom panel depicts findings at the subcortical surfaces with magnified panel demonstrating findings within the right caudate and putamen. A: anterior, P: posterior, S: superior, R: right (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.).

Table 3

Resting-state functional connectivity findings in epilepsy patients vs. controls relative to functional networks.

Functional Network	Weighted proportion of network involvement: Increased FC
Somatomotor	39.77
Salience	35.67
Visual	11.57
Default mode	5.84
Dorsal attention	5.57
Frontoparietal control	1.58

EZ-seeded whole brain lag analysis

When evaluating whole brain time lag relative to EZ at the group level, multiple foci of both increased and decreased lags were identified in patients compared to controls (Fig. 4). Foci of both increased and decreased lags were broadly distributed over the cortical and subcortical surfaces. Table 4 demonstrates the weighted proportion of network involvement of time lag findings by comparing the areas of altered time lag, weighted by their t-value, against a functional atlas (B. T. Yeo et al., 2011) again excluding Yeo's ‘limbic’ network (Gordon et al., 2016; Yeo et al., 2014). Again, all six functional networks demonstrated findings of either increased or decreased time lag with EZ in patients compared to healthy controls with the frontoparietal control, salience and default mode networks demonstrating both increased and decreased time lag with EZ.

Fig. 4

Results from whole brain time lag analysis. Top panel depicts cortical regions of both increased and decreased time lag when selecting the seizure onset node (EZ) as the seed for whole brain time lag analysis. Comparison was made at the group level across all patients compared to all control subjects. Colorbar demonstrates significant t-values for both the cortical and subcortical surfaces. Bottom panel depicts findings at the subcortical surfaces. A: anterior, P: posterior, S: superior, L: left. .

Table 4

Resting-state lag findings in epilepsy patients vs. controls relative to functional networks .

Functional Network	Weighted proportion of network involvement: Increased lag (i.e., hyposynchrony w/ RZE)	Weighted proportion of network involvement: Decreased lag (i.e., hypersynchrony w/ RZE)
Frontoparietal control	51.13	24.39
Dorsal attention	0	46.21
Salience	40.55	1.58
Default mode	4.20	23.15
Visual	0	4.67
Somatomotor	4.12	0

Discussion

Altered synchrony architecture in the epileptic brain

Neuronal synchrony mediates functional integration of brain regions as in the case of visual pattern recognition (Gray et al., 1989), conscious awareness (Engel and Singer, 2001) and attention (Fries et al., 2001) and also forms the physiologic basis of functional connectivity (J. D. Power et al., 2011; B. T. Yeo et al., 2011) although these investigations have been carried out on differing time and spatial scales. Prior synchrony analysis in epilepsy has demonstrated a range of findings within interictal data. These include increased phase lag index across irritative regions when compared against contralateral homologous regions using magnetoencephalography data (Nissen et al., 2016), increased synchronization likelihood within the temporal lobes in TLE patients versus epileptic patients without TLE using SEEG data (Bartolomei et al., 2013), as well as improved resection outcomes after the removal of all hypersynchronous regions defined in terms of mean phase coherence (Schevon et al., 2007) and phase synchronization (Ortega et al., 2008) using intracranial EEG data. More recently, Lagarde et al., described synchrony findings in epilepsy using SEEG data with results exactly aligning with those presented here, both at the nodal level within the epilepsy network as well as throughout the remaining non-epileptogenic brain (Lagarde et al., 2018). Here, we extend this prior work by capitalizing on the advantages afforded by rsfMRI data by considering synchrony in rsfMRI data sets compared to healthy controls, both with regard to directed relationships (i.e., between nodes of the epilepsy network) as well as with respect to the whole brain and its intrinsic connectivity architecture in the form of RSNs.

We considered two aspects of directed relationships: maximal cross-correlation and associated time lags of rsfMRI. The use of BOLD-derived measures permits whole-brain analyses, whereas comparison to control subjects permits an accounting for the spatial variability of the epileptic network as well as the variable SEEG electrode locations across patients. Furthermore, recent evidence suggests that lag structure at the timescale of rsfMRI is not only organized, but also reflects an underlying framework for the integration of functional information within and between networks forming an intrinsic connectivity architecture (Mitra et al., 2014). This may seem counterintuitive, given the historically low temporal resolution of rsfMRI data however cross-correlation can be considered a common approach for solving Time-Delay Estimation (TDE) problems (Hero A, 1998). The lag value can then be defined by maximizing the cross-correlation of BOLD signal over a range of possible lag values. This assumes the following: 1) time lags across brain regions are constant throughout the rs-fMRI acquisition, 2) the quantity of noise in the signal (SNR) and the number of time points are sufficient to afford an estimation. If these assumptions are met, the sampling interval (TR) has no effect on the accuracy and therefore, lags can be determined with a precision much finer than the TR (Jacovitti G, 1993). This, however, holds true only if the signals are band limited in the frequency range [-W, W], with TR≤1/(2 W). In this study, time series were low pass filtered at 0.08 Hz and TR=3.6 s, fulfilling this condition. As such, analysis of the lag structure between epilepsy network nodes (epileptic (EZ), propagation (PZ) or non-primarily involved (NPI) zones), identified multiple signatures of altered directed relationships in epilepsy. Given the range of localization related epilepsy syndromes considered, and comparisons with controls, our findings suggest that generalized features of epilepsy include 1) focal hypersynchrony (in the form of reduced time lags) within epileptogenic regions, 2) reduced integration between the NPI zones and brain regions subject to propagated ictal activity (manifest as increased time lags), and 3) diffuse hyposynchrony (increased time lags) in the rest of the brain. The isolation of nodes within the epileptic network has been demonstrated in lag analysis using intracranial EEG data (Warren et al., 2010). While capturing data at a distinctly different timescale, reduced lag times within the EZ node of the epileptic network using rsfMRI is consistent with this previous work. This finding is also in keeping with the focal preservation of structural connectivity in the EZ despite widespread reductions when investigating structural connectivity using SEEG-defined regions (Besson et al., 2017a.b).

These altered directed node-based time lag relationships are particularly interesting in the context of dramatically increased time lag found within all brain regions outside the epilepsy network (Fig. 2). This suggests a fundamentally altered time lag architecture in the epileptic brain compared to controls characterized by 1) decreased time lags within regions of seizure onset, 2) increased time lags between nodes of the epilepsy network and the rest of the brain and 3) globally increased time lags within all other brain regions. Previous work in autism identifies not only differences in the resting state lag architecture in the autistic brain but also a relationship between the degree of these changes and behavioral measures relevant to autism (Mitra et al., 2017). This would suggest a possible mechanistic link between lag and synchrony architecture and clinical disease manifestation.

Network steal hypothesis: proposed mechanistic explanation for the non-seizure clinical manifestations of epilepsy

When directed correlation was analyzed between SEEG-defined nodes of the epilepsy network using rsfMRI, no differences were observed between patients and controls. This is in contrast to prior evidence supporting differences in fMRI-based analyses of non-directed functional connectivity in clinically defined regions (Bettus et al., 2011; Stufflebeam et al., 2011; Weaver et al., 2013). Despite this, our findings are consistent with prior directed analyses using intracranial EEG (Dickten et al., 2016; Wang et al., 2017).

A range of alterations were observed, however, when assessing the relationship between the epilepsy network and the rest of the brain. This is consistent with previous findings of significant influence of EZ on non-epileptic regions when analyzing the directionality index using intracranial EEG data (Bettus et al., 2011). In particular, we found increased strength of correlation between the epileptogenic node and cortical regions known to be involved in the somatomotor, salience, visual and default mode networks (Fig. 2, Table 3). These findings are of particular significance in the context of other reports of these brain regions demonstrating decreased intranetwork correlation in epilepsy (de Campos, Coan, Lin Yasuda, Casseb, and Cendes, 2016; Liao et al., 2011; Pittau et al., 2012; Zhang et al., 2010).

Synthesizing these two reciprocal profiles of altered connectivity in the setting of epilepsy, together with the altered lag architecture identified herein, suggests the possibility of an underlying “network steal” phenomenon in which one network's function is compromised because its intrinsic connectivity is “stolen” by a coincident pathology. There is evidence to support this notion of reciprocal connectivity changes related to both the physiologic process of learning and memory (Rao et al., 2016) as well as pathologic conditions including depression (Rao et al., 2016), schizophrenia (Park et al., 2017), autism spectrum disorders (Kleinhans et al., 2008) and epilepsy (Dansereau et al., 2014).

In the context of our findings reported here, this network shunting may occur due to priming between nodes of the epilepsy network and other brain regions which may result from the altered lag architecture in epilepsy, serving to favor connectivity both within the epileptogenic network and between the epileptogenic network and other RSNs at the expense of connectivity within each RSN. While altered lag structure has not been described in all of the other conditions noted above, it has been described within autism spectrum disorders and has been found to correlate with a variety of behavioral scores relevant to autism (Mitra et al., 2017) suggesting a possible mechanistic link between lag structure and neurocognitive performance. Alternatively, there may be separate or additional causes of shunting contributing to the connectivity changes associated with other conditions beyond epilepsy.

We hypothesize that the non-seizure clinical manifestations of epilepsy (including motor discoordination, attentional and cognitive deficits, tendency toward depression and anxiety, and poor executive functions), manifests from this shunted connectivity away from physiologic networks and toward the pathologic epilepsy network. These concomitant neuro-psychiatric conditions have long been considered co-morbidities with epilepsy until recently when they have been proposed to be fundamentally related to the presence of epilepsy (Josephson et al., 2017) providing further support for this reciprocal relationship. Exploring the interplay between these non-seizure clinical manifestations of epilepsy and the connectivity profile of the epileptic brain further may not only enhance our understanding of epilepsy but may also further our understanding of other neuropsychiatric conditions.

Conclusions

Using SEEG-defined seizure onset and propagation locations, we identify consistent patterns of altered time lag structure and increased strength of correlation relative to EZ in patients with epilepsy compared to healthy controls using rsfMRI data. These patterns include 1) pathological decreased lag involving nodes of the epilepsy network, 2) global increased lag throughout the remainder of the brain and 3) findings of broadly distributed increased strength of correlation between the seizure onset zone and the remainder of the brain. This work brings closer together two disparate bodies of work in epilepsy, correlation analysis and synchrony. Despite these broad, global findings, consistent patterns of increased directed functional connectivity between EZ and the somatomotor, visual and primary attentional networks of the brain including the salience and default mode networks are identified. This, in the context of prior reports of decreased connectivity within these individual networks in the context of epilepsy, leads us to propose the network steal hypothesis as a possible mechanistic explanation for the non-seizure clinical manifestations of epilepsy including the panoply of motor, cognitive and neuropsychiatric challenges patients with epilepsy manifest. Together, these findings challenge us to expand our understanding of epilepsy to include the functional connectivity dynamics inherent within it in order to develop novel treatment opportunities for the disorder.

Funding

This work was supported by the following funding sources: 7T-AMI ANR-11-EQPX-0001, A*MIDEX-EI-13-07-130,115-08.38-7T-AMISTART, A*MIDEX ANR-11-IDEX-0001-02, CNRS (Centre National de la Recherche Scientifique) and The American Association of Neurological Surgeons.