BOLD-fMRI activity informed by network variation of scalp EEG in juvenile myoclonic epilepsy.

Epilepsy is marked by hypersynchronous bursts of neuronal activity, and seizures can propagate variably to any and all areas, leading to brain network dynamic organization. However, the relationship between the network characteristics of scalp EEG and blood oxygenation level-dependent (BOLD) responses in epilepsy patients is still not well known. In this study, simultaneous EEG and fMRI data were acquired in 18 juvenile myoclonic epilepsy (JME) patients. Then, the adapted directed transfer function (ADTF) values between EEG electrodes were calculated to define the time-varying network. The variation of network information flow within sliding windows was used as a temporal regressor in fMRI analysis to predict the BOLD response. To investigate the EEG-dependent functional coupling among the responding regions, modulatory interactions were analyzed for network variation of scalp EEG and BOLD time courses. The results showed that BOLD activations associated with high network variation were mainly located in the thalamus, cerebellum, precuneus, inferior temporal lobe and sensorimotor-related areas, including the middle cingulate cortex (MCC), supplemental motor area (SMA), and paracentral lobule. BOLD deactivations associated with medium network variation were found in the frontal, parietal, and occipital areas. In addition, modulatory interaction analysis demonstrated predominantly directional negative modulation effects among the thalamus, cerebellum, frontal and sensorimotor-related areas. This study described a novel method to link BOLD response with simultaneous functional network organization of scalp EEG. These findings suggested the validity of predicting epileptic activity using functional connectivity variation between electrodes. The functional coupling among the thalamus, frontal regions, cerebellum and sensorimotor-related regions may be characteristically involved in epilepsy generation and propagation, which provides new insight into the pathophysiological mechanisms and intervene targets for JME.

Introduction

Juvenile myoclonic epilepsy (JME) (Spencer, 2002) is traditionally considered as an idiopathic generalized epilepsy, characterized by myoclonic jerks, tonic-clonic seizures and infrequent typical absences (Woermann et al., 1999; Genton et al., 2013), and stereotyped bilaterally synchronous spike-wave discharges (GSWDs) or polyspike-wave activity (Panayiotopoulos et al., 1994; Koepp et al., 2014). In recent studies, JME is no longer seen as a homogeneous disease but a system disorder of the brain with ictogenesis in distributed bilateral networks involving primarily the thalamus and other cortical areas (Baykan and Wolf, 2017). Altered neuronal interactions and cortico-cortical network re-arrangement were regarded to play a crucial role in triggering the epileptic activity (Clemens et al., 2013; Jiruska et al., 2013). Therefore, examining the abnormal dynamic network organization may provide new insight enabling a further understanding of the epileptic brain.

For the initiation and propagation of generalized discharges, cumulative evidences showed that aberrant thalamo-frontal circuit is the key element that contributes to GSWDs generation (Blumenfeld, 2003; Moeller et al., 2008; O'Muircheartaigh et al., 2011; Jiang et al., 2018). Moreover, it was suggested that widespread hyperconnectivity in sensorimotor and frontal cortex acts as the excitatory driver in discharge propagation (Vollmar et al., 2011; Clemens et al., 2013; Lee et al., 2017). Conversely, epileptic discharges, i.e., the sudden, transient disturbances of brain activity, can lead to abnormal organization of the diffuse brain networks (Engel et al., 2013). That is, the whole brain network is responsible for the electrographic phenomena associated with epilepsy seizures (Spencer, 2002); meanwhile, the discharges may be driven by the abnormal patterns of brain networks (Terry et al., 2012). Previously, dynamic functional connectivity has been used to accurately locate the seizure onset based on the intracranial or scalp EEG (van Mierlo et al., 2011; Staljanssens et al., 2017). In our recent study, according to the time-varying network patterns during different stages of the interictal discharges, it was showed that the electrode node close to the epileptogenic zone was initially involved with the strongest out flows before the interictal discharges (Zhang et al., 2017). In addition, the highest temporal variation of network synchronization has shown good performance as predictors of epilepsy-related fluctuations (Abreu et al., 2018). The dynamic network reorganization was proved to relate to the presence of epileptic activity through synchronous and asynchronous neurophysiologic process (Jiruska et al., 2013; Khambhati et al., 2015). Therefore, clinically, the knowledge of dynamic network-related brain fluctuations and functional couplings is essential to understand the spatiotemporal dynamics underlying epilepsy generation and propagation.

Technically, the combination of EEG and fMRI has been a powerful approach to investigate the BOLD fluctuations associated with the epileptic discharges. The common method is based on the visual inspection of the discharges of interest (i.e., the interictal epileptic discharges (IEDs) and the seizures), with the purpose of obtaining a predictor of discharge-related BOLD responses. EEG-fMRI fusion has provided a unique advantage in identifying epilepsy localization (Salek-Haddadi et al., 2006; Tyvaert et al., 2008). Also, in generalized epilepsy, spatial involvements associated with the generalized discharges have been found, including BOLD activation in thalamus, mesial frontal areas and deactivation in the parietal and anterior frontal regions (Aghakhani et al., 2004; Gotman et al., 2005). However, traditional EEG-fMRI fusion in the field of epilepsy has not fully explored the temporal characteristics of EEG and the architecture of functional network between electrodes that may provide important spatiotemporal dynamics of the brain. Therefore, the BOLD activity associated with the network dynamics of scalp EEG may extend the traditional discharge-informed fMRI analysis. In addition, previous studies also showed that the functional couplings between the regions in cortico-subcortical circuits may also interact with the cortical dynamics, thus resulting in low-frequency epileptic oscillation (Destexhe et al., 1999; Blumenfeld and McCormick, 2000). Therefore, we hypothesize that functional couplings between BOLD responding regions can be modulated by the cortical dynamic network organization, and the modulatory effect may have important implications for understanding the directional circuits in the dynamics of rapid spread of epileptic discharges (Blumenfeld, 2003; Sinha et al., 2019).

The aims of the current study included: (1) determining the fMRI activity associated with EEG dynamic network organization in JME; (2) investigating the directional functional couplings among target regions modulated by EEG dynamic network variation. Therefore, simultaneous EEG-fMRI recordings were conducted, and time-varying scalp network was generated based on the adapted directed transfer function (ADTF) between EEG electrodes. Different levels of the network variation were selected to predict the BOLD response. Furthermore, the EEG-dependent functional couplings among BOLD responding regions were analyzed using an extended psychophysiological interaction (PPI) analysis, which was a linear regression framework to identify the contribution of one region to another and incorporated the modulation of a third region or experimental factors (Friston, 1997). Finally, we evaluated the clinical relevance of the modulatory effects in JME, and key nodes in cortical and subcortical circuits were revealed.

Materials and methods

Participants

Eighteen JME patients (12 females; mean age: 21 years; standard deviation: 7 years; age range: 15–34 years) were recruited in this study. All patients were diagnosed by neurologists based on the clinical and seizure semiology information consistent with International League Against Epilepsy (ILAE) guidelines (Engel, 2001). No structural abnormalities were observed by routine examinations of CT and MRI scanning. Additionally, 3–6 Hz generalized spike-wave or polyspike-wave discharges were found by 24-h scalp EEG recordings. All participants gave written informed consent in accordance with the Declaration of Helsinki. This study was performed according to the guidelines approved by the Ethics Committee of the University of Electronic Science and Technology of China (UESTC). Data from all participants were previously included in a prior article in which authors demonstrated the linear and nonlinear discharge-affecting networks using IEDs onset information in eigenspace maximal information canonical correlation analysis (Dong et al., 2016), whereas the present study focused on the relationship between the network dynamic organization of EEG and voxel-level BOLD responses, as well as the network-level functional coupling.

Simultaneous EEG-fMRI recording

In this study, all participants had simultaneous EEG-fMRI recording in UESTC. FMRI data were collected using a 3-T MRI scanner (Discovery MR750, GE) in the Center for Information in Medicine of UESTC. All participants were instructed simply to keep their eyes closed and remain still. High-resolution T1-weighted images were acquired using a 3- dimensional fast spoiled gradient echo (T1-3D FSPGR) sequence (TR/TE = 5.936 ms/1.956 ms, flip angle [FA] = 9°, matrix = 256 × 256, field of view [FOV] = 25.6 × 25.6 cm2, slice thickness = 1 mm, no gap, 152 slices). Resting-state functional MRI data were acquired using gradient-echo echo planar imaging sequences (TR/TE = 2000 ms / 30 ms, FA = 90°, matrix = 64 × 64, FOV = 24 × 24 cm2, slice thickness/gap = 4 mm/0.4 mm), with an eight channel-phased array head coil. A 510-s resting-state scan was collected for each run, and five repeated runs were conducted for each JME patient.

Simultaneous EEG data of the JME patients were recorded using a 64-channel MR compatible EEG cap (Neuroscan, Charlotte, NC) according to the 10–20 standard system with a reference at the Fcz position. The amplifier (Neuroscan, synAmps2) was placed outside the scanning room, and the sampling rate was set at 5000 Hz. Electrode impedances were lowered to below 10 kΩ prior to recording. The EEG recording was synchronized with the MR scanner's internal clock to ensure the removal of the gradient artifacts in the EEG analyses. During the recording, all participants were instructed to close their eyes and relax without falling asleep.

EEG-fMRI data preprocessing

All fMRI data were preprocessed using SPM12 (Statistical Parametric Mapping, http://www.fil.ion.ucl.ac.uk/spm/) and NIT (http://www.neuro.uestc.edu.cn/NIT.html) toolboxes (Dong et al., 2018). The first five volumes were discarded for the magnetization equilibrium from all fMRI scans. The remaining images were slice-time corrected, and spatially realigned to the first volume. All subjects have <1 mm for head movement and 1°for head rotation during MRI scanning. The individual T1 images were coregistered to the functional images, and then segmented and normalized to the Montreal Neurologic Institute (MNI) space. Then, functional images were spatially normalized using T1-based transformation, resampled to 3 × 3 × 3 mm3 voxels, and spatially smoothed with a 6 mm full-width half maximum (FWHM) Gaussian kernel.

EEG data were preprocessed using Curry 7 software (Compumedics Neuroscan). MR gradient artifacts were removed by subtracting the averaged scanner artifact template from the continuous EEG recordings based on the scanner markups (Allen et al., 2000); then, the EEG data were bandpass filtered (1–45 Hz) and down-sampled to 250 Hz. The ballistocardiogram (BCG) artifacts were corrected using the OBS-based BCG correction using the ECG channel (Niazy et al., 2005). Then, we used ICA to manually reject movement-related and residual ballistocardiac artifacts (Luo et al., 2010). After the data were marked into periods of 2 s (one fMRI volume acquisition), the sampling points with obvious motion-related artifacts in each period were visually inspected and excluded for further analysis. Finally, the preprocessed EEG were re-referenced to the neutral reference using the reference electrode standardization technique (REST, www.neuro.uestc.edu.cn) (Yao, 2001; Dong et al., 2017). Onsets of IEDs for each patient were marked by two experienced neurologists. The run with the most global discharges was selected for the following studies. 19 channels approximating the standard electrode locations in the 10–20 system with the labels FP1, FP2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T7, T8, P7, P8, Fz, Pz, Cz were selected as nodes for the following network calculations.

Overview of EEG-fMRI analysis based on the network variation of scalp EEG

The proposed scalp network variation-informed fMRI analysis can be summarized as follows (Fig. 1). The time-varying scalp network was first constructed from the preprocessed EEG using ADTF. Then, the variation of the information flow of the ADTF networks over 1–20 Hz in each 2 s window was extracted to construct the scalp network variation time series. The values of the network variation time series above the mean, as well as the values above one standard deviation from the mean were used as the regressors in GLM-based EEG-fMRI analysis, respectively. Thus, statistical contrast maps associated with different levels of the network variation measure were obtained from the above GLM analysis.

Fig. 1

Overview of the proposed EEG-fMRI analysis based on the network variation of scalp EEG. (A): The time-varying scalp network was constructed from the preprocessed EEG using ADTF. (B): Variation of the information flow of the ADTF between electrodes in each 2 s time window was extracted to construct the network variation time series. (C): The values of network variation time series above the one standard deviation and the mean were used as the regressors in GLM-based EEG-informed fMRI analysis, respectively. (D): Statistical maps were obtained from the above GLM analysis.

Time-varying network between electrodes using ADTF

The directed transfer function (DTF) was used to explore functional connectivity between EEG channels. A multivariate autoregressive model (MVAR) was constructed to represent the multivariate dataset as a linear combination of its own past plus additional uncorrelated white noise. The MVAR was characterized by the following:where X(t) is the data vector over time, A(i, t) are the matrices of the time-varying model coefficients established by the Kalman filter algorithm (Arnold et al., 1998), E(t) is multivariate independent white noise and p is the model order that was chosen by the Schwarz Bayesian Criterion (Schwarz, 1978). The observation and state equations of this algorithm were solved by the recursive least squares (RLS) algorithm with forgetting factor (Campi, 1994).

As the function A(t) has been acquired, the DTF function, H(f), can be obtained from the MVAR model by transforming Eq. (1) into the frequency domain. The DTF is calculated by the following:

The elements of the transfer matrix, H(f), represent the directional causal interaction from the jth to the ith element at frequency f. The normalized ADTF, which describes the directional information flow from jth to ith, is defined by the elements of the transfer matrix in the spectral domain as follows:

To evaluate the total information flow from one single node, the so-called integrated ADTF was defined by summing the ADTF values over the frequency bands of interest, and normalized to be between (0,1).where the [f1, f2] is the frequency interval. In this paper, the frequency band 1–20 Hz was selected according to the spectral distribution of the IED trials (2 s before and after the IEDs onset for one trial), including the main IEDs covering the spike waves and sharp waves (Supplemental Material, Fig. S1).

Time-varying network variation extraction

After the time-varying network of scalp EEG was obtained, we calculated the variation of information flow of each node (F) in successive 2 s windows, which was consistent with the fMRI temporal resolution. F was calculated by the following:where Θ2(t) is the information flow from node j to node i, and m is the number of time points in each window (i.e. 500 in this study). Then, the sum of F of all nodes (∑ ∑ F) in the network was used to generate the new network variation time series (Supplemental Material, Fig. S2). In this time series, larger value means higher network variation. Therefore, in order to reduce temporal resolution of the network variation series, we calculated the upper envelope of signals using linear interpolation between the neighboring peaks. Thus, network variation time series of scalp EEG were generated for each subject and used for the following GLM-based EEG-fMRI analysis. In addition, to further diminish head motion influence, we checked the FD (frame-wise displacement) values during fMRI preprocessing, and time points with large FD (FD > 0.5) in the network variation signal were further processed through interpolation of the neighboring points (only 0.15% time points were included in this step). FD is a measurement to evaluate the head motion characteristics during fMRI scanning. Here FD was calculated as suggested by Power et al. (Power et al., 2012). The FD values of each subject were demonstrated in Supplemental Material, Fig. S3.

EEG-fMRI analysis informed by network variation of scalp EEG

After the network variation time series of scalp EEG were constructed, two measurements denoting the variation and medium variation were extracted from the time series (Supplemental Material, Fig. S2). To identify the BOLD response when the scalp EEG network had high temporal variation, the values above one standard deviation (i.e., higher than the mean plus the standard deviation) in the generated EEG time series were selected and taken as the regressors in the GLM-based EEG-informed fMRI analysis for each subject. In addition, the values above the mean value, which denoted the medium widespread variation of the time-varying network, were also extracted from the time series. The EEG-informed fMRI analysis was conducted using the Neuroscience Information Toolbox (NIT v1.2, http://www.neuro.uestc.edu.cn/NIT.html). The regressors related to different levels of scalp network variation were included after convolving with 4 canonical hemodynamic response functions (HRFs) peaking at 3, 5, 7, and 9 s. Additionally, nuisance signals (six motion parameters, linear drift signal, as well as the mean white matter and cerebrospinal fluid signals) were regressed out through the design matrix. Slow signal drifts with a period longer than 128 s were also removed (i.e., frequencies below 0.008 Hz). A one-sample t-test was then performed using the statistical images resulting from single-subject contrasts to examine the effect on the BOLD signal at the population level for the JME group. Statistical parametric maps were thresholded at an uncorrected voxel height P < 0.001.

In addition, traditional GLM-based discharge-informed fMRI analysis was conducted using the manually determined IEDs onsets as the regressors. This analysis was performed to elucidate whether the scalp EEG network variation-informed fMRI method provides extra information compared to the traditional EEG-fMRI analysis. Similarly, nuisance signals (six motion parameters, and the linear drift signal, as well as the mean white matter and cerebrospinal fluid signals) were regressed out through the design matrix. One-sample t-test was also conducted on the BOLD responses in the group with the threshold at P < 0.001.

Functional coupling modulated by EEG variation

After the EEG-informed fMRI analysis, BOLD responses correlated with the high level and medium level of dynamic EEG variation were obtained. Then, the regions with significant high BOLD activation and deactivation in the one-sample t-test were selected as the regions of interest (ROI, 3*3*3 voxels), and the mean BOLD signals within the ROIs were extracted. The EEG-dependent functional coupling for all pairs of ROIs was examined based on the extended PPI analysis (Friston, 1997; Di et al., 2017), which was used to examine the modulatory interaction for the EEG-network variation and the functional connectivity between regions. PPI aims to identify the contribution of one region to another with the modulation of other factors, and includes correlation PPI terms and modulatory interaction term. In present study, the correlation PPI terms were defined using the BOLD time series of the ROIs, and the EEG-network induced signals generated by convolving the high EEG-network variation signals with the canonical HRF. The modulatory interaction term was the point-by-point multiplication between the BOLD signal in one region and the EEG-network induced signals. In this study, the BOLD-related terms were calculated without deconvolution but with modulation signal centered (Di et al., 2017). The modulatory interaction generation model was as follows:where xand y represented the BOLD time series of the two ROIs, and x was the high EEG-network induced signal. β3 values indicated the modulation effect resulting from x on the functional coupling from region1 to region 2. Therefore, positive modulation effects denoted increased functional coupling between regions modulated by EEG, while negative modulation effects denoted decreased functional coupling between regions because of EEG modulation.

All pairs of ROIs extracted from the EEG-informed fMRI analysis were involved in the modulatory interaction analysis. Nuisance signals (six motion parameters, and linear drift signal as well as the mean white matter and cerebrospinal fluid signals) were regressed out. Then, group statistical analysis using one-sample t-test was conducted for β3 values. Before statistical analysis, the same regions in bilateral hemispheres, as well as all cerebellum regions were grouped as new ROIs and the averaged modulatory effects in one ROI was calculated. The statistical threshold was set to P < 0.05, corrected for family-wise errors (FWE), to show the significant EEG-based modulation effect on the functional coupling between the new ROIs.

Relationship between modulatory interaction and clinical features

Pearson correlation or Spearman correlation was conducted between the modulatory functional couplings and the clinical features according to the data distribution. The clinical features included the age of epilepsy onset, the duration, as well as the seizure frequency (times per month) of the patients (Table 1). The EEG-fMRI modulation effect included the input and output effect of all ROIs in the above modulatory interaction analysis. The input effect for a particular ROI was the sum of the modulation effects coming from the other ROIs to this ROI, while the output effect of one ROI was the sum of the modulation effect from this ROI to the other ROIs.

Table 1

Detailed demographic information and clinical characteristics of JME patients.

No.	Gender	Age (year)	Frequency of GSWDs (Hz)	No. of volumes with global GSWDs in selected run	Age at seizure onset (year)	Family history	AEDs	Seizure frequency (times per month)
1	F	17	2 Hz	1	10	–	VPA	1
2	F	17	3 Hz	20	14	–	LTG	1
3	F	33	6 Hz	6	20	–	VPM	3
4	M	22	4 Hz	3	8	–	VPM	2
5	F	19	3 Hz	6	12	Uncle with	MgV	1
						GTCS
6	F	20	2 Hz	111	6	–	VPM/LTG	6
7	M	15	3–3.5 Hz	12	5	–	–	4
8	F	22	3 Hz	2	14	–	VPA	2
9	F	17	3–3.5 Hz	2	3	Sister with JME	VPA	0.5
10	F	17	2 Hz	8	13	Sister with JME	VPA	0.5
11	F	29	3–3.5 Hz	2	10	–	TCM/VPM	5
12	M	18	5 Hz	2	14	–	VPM	1
13	F	27	4 Hz	3	16	Daughter with GTCS	–	3
14	F	21	2 Hz	117	11	–	VPM	0.5
15	M	10	4 Hz	2	5	Brother with JME	VPM	1
16	M	13	3–3.5 Hz	4	9	–	VPA	4
17	M	34	3.5–4 Hz	4	14	–	TCM/VPM	5
18	F	34	3 Hz	6	18	–	–	1

Note:

GSWDs: generalized spike-wave discharges; GTCS: generalized tonic-clonic seizures; AEDs: Antiepileptic drugs; VPA: valproic acid; LTG: lamotrigine; VPM: valpromide;

MgV: magnesium valproate; TCM, traditional Chinese medicine.

Results

BOLD response informed by temporal variation of scalp EEG network

As the time-varying network was acquired by the ADTF analysis between EEG electrodes, two types of temporal network variation, i.e., high variation and medium variation, were extracted to perform the GLM-based variation-informed BOLD response estimation. Group level analysis showed that positive BOLD responses highly correlated with the high EEG-network variation were mainly located in the sensory and motor regions, i.e. middle cingulate cortex (MCC), supplemental motor area (SMA), paracentral lobule, as well as the cerebellum, precuneus, and inferior temporal lobe (Fig. 2, Table 2). At the subcortical level, positive BOLD responses in the thalamus, and fewer voxels in the caudate and bilateral insula were also found. In contrast, negative BOLD responses associated with the medium EEG-network variation were found in the bilateral frontal lobe, anterior cingulate cortex (ACC), the postcentral, middle temporal, inferior parietal and occipital areas (Fig. 3, Table 3). These positive and negative BOLD responses were concordant with most of the previous IGE discharge-informed EEG-fMRI reports (Aghakhani et al., 2004; Gotman et al., 2005; Hamandi et al., 2006), but they were detected with different levels of scalp network organization. Furthermore, traditional GLM-based EEG-fMRI analysis using IEDs was also conducted, and the discharge-related BOLD activations in the cortical and subcortical regions, including the MCC, ACC, thalamus, putamen, and inferior temporal areas were found (Supplemental Material, Fig. S4).

Fig. 2

The BOLD response correlated with high EEG-network variation.

Table 2

Regions with significant BOLD response correlated with high EEG-network variation in JME (P < 0.001).

	MNI coordinates
Brain regions	x	y	z	Peak T-value	Cluster voxels
Cingulum_Mid_R	−1	18	33	7.2692	82
Cingulum_Mid_L	−3	15	33	8.4792	74
Supp_Motor_Area_L	-1	5	72	5.3883.	38
Supp_Motor_Area_R	1	3	71	5.2105	35
Paracentral_Lobule_L	-3	−21	78	5.8027	31
Paracentral_Lobule_R	1	−42	70	5.2676	23
Precuneus_R	9	−48	78	6.0834	21
Precuneus_L	−2	−78	47	5.1229	21
Thalamus_L	−14	−7	4	5.0042	35
Caudate_L	−10	3	12	4.1016	18
Insula_R	35	−12	4	5.002	18
Insula_L	−33	−18	9	5.6433	18
Temporal_Inf_R	57	−60	−18	5.3501	34
Cerebelum_4_5_L	−16	−45	−16	4.6384	71
Vermis_6	3	−67	−9	7.0837	67
Vermis_4_5	3	57	−9	8.8713	72
Cerebelum_4_5_R	12	−50	−16	5.585	47
Cerebelum_8_R	21	−62	−45	4.9414	55
Cerebelum_6_L	−9	−64	−15	4.3198	43
Cerebelum_6_R	8	−64	−14	4.9946	36

Fig. 3

The BOLD response correlated with medium EEG-network variation.

Table 3

Regions with significant BOLD response correlated with medium EEG-network variation in JME (P < 0.001).

	MNI coordinates
AAL regions	x	y	z	Peak T-value	Cluster voxels
Frontal_Sup_Medial_L	−4	52	5	−6.4968	69
Frontal_Sup_Medial_R	12	51	15	−6.9962	64
Cingulum_Ant_R	5	52	13	−6.867	51
Cingulum_Ant_L	−4	52	4	−6.3922	31
Frontal_Med_Orb_L	−11	64	−5	−5.729	42
Frontal_Med_Orb_R	4	62	−2	−4.8874	37
Frontal_Sup_L	−12	58	22	−4.68	42
Frontal_Sup_R	18	33	39	−5.6635	38
Postcentral_L	−60	−15	24	−7.4313	46
Parietal_Inf_R	40	−58	48	−5.633	41
Angular_R	36	−60	48	−5.8183	32
Temporal_Mid_R	63	0	21	−6.4505	35
Occipital_Inf_L	−41	−80	−6	−4.9689	46
Lingual_L	−9	−78	−6	−5.4952	46
Fusiform_L	−21	−80	−8	−5.2044	32

Functional couplings based on the modulatory interaction analysis

The regions representing significant BOLD activation and deactivation were involved in the modulatory interaction analysis. Then, group-level modulation-dependent functional couplings with statistical significance between ROIs were illustrated in Fig. 4 (P < 0.05, corrected for FWE). Generally, negative modulatory interactions dominated the functional couplings between regions (Fig. 4 (A)). Under EEG modulation, the frontal regions and ACC showed decreased functional coupling to the SMA, paracentral lobule regions, and the cerebellum showed decreased functional coupling to both the frontal (frontal_sup and frontal_medial) and sensorimotor-related (SMA and paracentral lobule) areas. Moreover, the thalamus showed decreased connectivity to the frontal regions. In addition, a few positive modulation interactions were also found (Fig. 4 B), which represented increased functional coupling within the frontal regions (frontal_sup, frontal_medial, and ACC), as well as the increased connectivity from the frontal regions to the cerebellum.

Fig. 4

Modulatory interaction for EEG-network variation and BOLD functional coupling. (A) Negative modulation effects exist from thalamus to frontal areas, from cerebellum to frontal and sensorimotor areas, as well as from frontal to sensorimotor areas. (B) Positive modulation effects exist within frontal areas, and from frontal to cerebellum. (C) Summarization of EEG-dependent functional coupling resulting from the modulatory interaction among four important epileptic central districts. Abbreviations of the regions according to AAL atlas: Front_S: Frontal_Sup; Front_M: Frontal_Med_Orb; ACC: anterior cingulate cortex; Front_Sm: Frontal_Sup_Medial; SMA: Supplementary motor area; Parac_l: Paracentral_lobule.

Taken together, four important epileptic central districts, including the thalamus, frontal regions, sensorimotor-related regions (i.e., SMA, and paracentral lobule), and cerebellum were identified, and the modulation-dependent functional couplings among them were shown in Fig. 4 (C). The results demonstrated that the main circuit in JME was associated with the epileptic dynamic network organization.

Relationship between modulatory interactions and clinical features

The EEG modulation effects were correlated with age of epilepsy onset, the duration of epilepsy, and the seizure frequency (times per month) (P < 0.05, Fig. 5). Significant negative relationships were found between the negative output effects from the thalamus, ACC, with the age of epilepsy onset, as well as the negative input effects into precuneus with age of epilepsy onset. Moreover, more prominent negative input effects in cerebellum, the medial and superior frontal areas were associated with longer epilepsy duration. In addition, frontal areas also showed larger negative input effects together with the seizure frequency increasing.

Fig. 5

Relationship between the modulation effects and clinical features. Significant negative correlation was found between negative modulation effects with age of epileptic onset, duration and seizure frequency (times per month). Note that the negative effect on y-axis represents reverse EEG modulation resulting to decreased functional coupling between regions, and the modulation becomes more prominent along with the negative effect enhancing.

Discussion

In this study, temporal variation of the connectivity between EEG electrodes was used for the first time to predict the BOLD response in JME. BOLD activation and deactivation was associated with different levels of network variation. Then, functional coupling between BOLD responding regions was conducted based on the EEG network-dependent modulatory interaction analysis, and predominantly directional negative modulatory effects among the thalamus, cerebellum, frontal and sensorimotor-related areas were found.

Combined EEG-fMRI has been successfully applied in epilepsy studies to capture the pathophysiological mechanisms associated with epileptiform discharges. Generalized discharge-related BOLD changes commonly showed thalamocortical activation and deactivation in the posterior and anterior cortex (Aghakhani et al., 2004; Gotman et al., 2005; Laufs et al., 2006). The thalamus, midfrontal and parietal areas are known to play an important role in the generation and spread of generalized epileptic activity (Meeren et al., 2002; Hamandi et al., 2006), and the BOLD activations in these regions were supposed to stem from the synchronized neuronal activity represented by the spike-and-wave discharges (Gotman et al., 2005). In this study, consistent results were found applying temporal variation of scalp EEG network to predict the BOLD response. The results also showed that the BOLD response was correlated with the level of network variation. The thalamus, and sensorimotor areas, including the SMA, MCC, paracentral lobule were involved with BOLD activation related to high temporal variation of the scalp network. The high variation of the scalp network may be partly due to the epileptic discharges, as we found about 47% IED onsets were included with the high network variation (Supplemental Material, Fig. S5) in this study. Moreover, comparing to a healthy group of 17 subjects, significant larger network variation (two-sample t-test, p = 0.0021), and increased activations in MCC, cerebellum, SMA in JME group were found (Supplementary material, Tables S1, S2). Therefore, we thought the high scalp network variation may be associated with the neural activity synchronization and desynchronization resulting from epileptic activity as well as the altered functional organization in JME. Compared with the results from the traditional IED-informed fMRI analysis, the complementary network-informed BOLD responses enabled a further understanding of the epileptic brain at the level of dynamic functional organization. Another interesting finding was the fMRI deactivations in bilateral frontal, parietal, and occipital areas both in JME and healthy group, which were correlated with medium EEG-network variation. The dominant deactivation may be attributed to the rhythm oscillation. As previous studies have reported, low frequency is negatively correlated with BOLD changes (Laufs et al., 2006; Scheeringa et al., 2008; Scheeringa et al., 2011), which benefits the normal network organization (Gotman et al., 2005). In the current study, a broad low frequency band was selected, which may contribute to the negative BOLD response. Therefore, the BOLD deactivation may be a basic state of the human brain, which originated from rhythmic oscillations and thus served for network reorganization and switching, while in JME, the basic state was altered (Supplemental Material, Tables S3, S4).

The modulatory functional coupling associated with network organization of scalp EEG provides new insight into understanding the cortical and subcortical circuits in epileptic brain. In this study, the modulatory interaction analysis identified the functional couplings among four important epileptic central districts, including the thalamus, frontal regions, sensorimotor-related regions, and cerebellum. Predominantly negative modulation effects may be related to the network inhibition hypothesis of subcortical structures in human epilepsy; that is, generalized discharges disrupt function in thalamus and cerebellum, then leading to widespread inhibition to the other cortical regions (Charlton et al., 1977; Norden and Blumenfeld, 2002). Interestingly, correlations were found between the modulatory effect of thalamus, cerebellum and clinical features, which points to the crucial role for subcortical structures in epilepsy pathology (Norden and Blumenfeld, 2002). Meanwhile, the deactivated frontal areas may act as the relay in the epilepsy circuit providing the excitatory drive to catalyze the disruption of the motor and premotor areas, resulting in the myoclonic seizures in JME, while the cerebellum output may be involved in controlling the generalized spike-wave discharge occurrences by regulating the GABA (Kros et al., 2015). The directional circuits integrated the mainstream ideas concerned in JME, and further, the key nodes in these circuits with clinical features were linked. Correlations between the negative modulatory effects of the high-order association cortex (the ACC, precuneus) and clinical features may be associated with the impaired consciousness and executive functioning in JME (Blumenfeld, 2012), while the weakened modulation in the case of earlier onset age may be related to neural plasticity during epileptic brain development. Moreover, the medial and superior frontal gyrus demonstrated more prominent modulation effect accompanied with longer epilepsy duration as well as increased seizure frequency in epilepsy patients. Here, epilepsy duration and seizure frequency in patients have significant positive correlation (R = 0.524, P < 0.05), representing two interdependent features of epilepsy. As the key area of the cognitive system, frontal gyrus was consistently involved in JME patients with functional and structural abnormality (Meschaks et al., 2005; Cao et al., 2013), and the interaction between cognitive system and motor system facilitates the seizures (Vollmar et al., 2011). The functional couplings between frontal and other regions in this study may be important biomarker of severity of the disease. Additionally, the positive modulation within the frontal areas represented increased functional coupling within one module, suggesting hyper-connectivity in frontal areas in JME. The positive modulation effect from the frontal areas to the cerebellum, may be associated with the feedback loop. Therefore, we inferred that the epileptic modulation circuit may commence within the thalamus and cerebellum, and then through inhibitory outputs to the frontal and frontoparietal area, consequently leading to the dysfunction of cortex.

Here we investigated the correlation between BOLD activity and EEG functional connectivity representing dynamic brain states. This approach could be a supplement to the traditional EEG-fMRI applications in epilepsy, which mostly depends on the onsets of generalized spike-wave discharges identified by neurologists. Moreover, the modulatory interaction analysis gave new insight to explore the functional circuits in epilepsy. We envisage in the future these methods may be incorporated with source imaging based on high-density EEG to infer on the focus and directed propagation of the epileptic activity. Moreover, this kind of multivariate, time-varying, and frequency-resolved measurement may constitute a refinement of the brain dynamics before and after discharge events, and provide new perspective for different seizure types of the heterogeneous disorder. In addition, network variation-based epileptic activity identification may weaken its sensitivity to MRI artifacts, which disturbs the simultaneously acquired EEG signals, as the gradient and BCG artifacts are generally considered relatively stable during the consecutive occurrences.

However, there are several methodological issues and limitations in this study. The first limitation is the lack of a direct relationship between the scalp network variation and epileptogenicity. The scalp time-varying network is an indirect measurement and the consequence may be influenced by the epileptic discharges and other disturbance, such as head motion. To verify the discharge detection ability of the network variation, we compared the IEDs onsets and the time points of high network variation for each subject, and about 47% IEDs onsets were involved. In addition, although a lot of measures were taken to avoid movement influence, correlations were still found in 4 subjects between head movement signal and network variation. However, correlations between head movement and IEDs onsets were also found in these subjects. Therefore, removing the movement artifact completely from the scalp dynamic network is hard, especially for epilepsy patients with IEDs during MRI scanning. In the future studies, additional correction using a motion-tracking system may provide promising results for removal of the artifacts (Qin et al., 2009; LeVan et al., 2013). The second possible concern is the selection of high and medium network variation to predict BOLD responses. In addition to the standard deviation and mean values chosen as different levels of the network variation, we also constructed the time series by selecting different proportions of the network variation (from 10% to 80%). Similar results were found with the proportion of top 10%, 20%, and 30% network variation as the high network variation, and the BOLD responses tended to be negative with increasing proportions. In addition, we applied PPI analysis to detect the modulation effect of the dynamic scalp network between EEG electrodes on the BOLD connectivity. This is the first attempt in the application of modulation interaction from one modality to another modality. Finally, considering that most epileptic discharge power occurred below 20 Hz and the narrow frequency bands may vary across patients, a broad lower frequency band 1-20 Hz rather than the endogenous EEG frequency bands (i.e. theta, alpha, beta, gamma) was extracted during time-varying EEG network construction. Additional analyses in future work using narrow patient-specific frequency bands would be interesting to explore epileptic activity.

Conclusion

In conclusion, the study provided a novel approach to predict fMRI activity based on the inherent dynamic network characteristics of scalp EEG, and BOLD responses correlated with different level of the network variation were found. Furthermore, the functional couplings among the thalamus, frontal regions, cerebellum, and sensorimotor-related regions were involved associating with dynamic EEG organization, which may be an important representation in epilepsy generation and propagation. This study extended the traditional EEG-fMRI fusion in epilepsy, and may provide new insight into the understanding of the pathophysiological mechanism and intervene target for JME.