Disruption, emergence and lateralization of brain network hubs in mesial temporal lobe epilepsy.

Hubs of brain networks are brain regions exhibiting denser connections than others, promoting long-range communication. Studies suggested the reorganization of hubs in epilepsy. The patterns of connector hub abnormalities specific to mesial temporal lobe epilepsy (mTLE) are unclear. We wish to quantify connector hub abnormalities in mTLE and identify epilepsy-related resting state networks involving abnormal connector hubs. A recently developed sparsity-based analysis of reliable k-hubness (SPARK) allowed us to address this question by using resting state functional MRI in 20 mTLE patients and 17 healthy controls. Handling the multicollinearity of functional networks, SPARK measures a new metric of hubness by counting the number (k) of networks involved in each voxel, and identifies which networks are actually associated to each connector hub. This measure provides new information about the network architecture involving connector hubs and a realistic range of k-hubness. We quantified the disruption and emergence of connector hubs in individual epileptic subjects and assessed the lateralization of networks involving connector hubs. In mTLE, we found pathological disruptions of normal connector hubs in the mTL and within the default mode network. Right mTLE had remarkably higher emergence of new connector hubs in the mTL than left mTLE. Different patterns of lateralization of the salience network involving the abnormal hippocampus were found in right versus left mTLE. The temporal, cerebellar, default mode, subcortical and motor networks also contributed to the lateralization of hippocampal networks. We finally observed an asymmetrical connector hub reorganization and overall regularization of epilepsy-related resting state networks in mTLE, characterized by the disruption of distant connections and the emergence of local connections.

Introduction

Mesial temporal lobe epilepsy (mTLE) is characterized by recurrent seizures originating in the mesial temporal structures - amygdala, hippocampus or entorhinal cortex (Engel Jr, 2001). For drug-resistant mTLE patients, surgical resection is often offered when the epileptic focus is well characterized. Long term post-surgical seizure freedom is obtained in ∼65% of patients (Harroud et al., 2012), suggesting the presence of epileptic networks involving other regions. Simultaneous electroencephalography (EEG) and functional MRI (fMRI) studies revealed hemodynamic responses to interictal epileptic discharges beyond the epileptic focus (Kobayashi et al., 2009; Khoo et al., 2017). Functional connectivity (FC) studies brain networks by looking at the temporally correlated fluctuations in the blood oxygen-level dependent (BOLD) fMRI signals in remote brain regions. Using FC, resting state networks (RSNs) are consistently found within the healthy population, such as the default mode network (DMN) (Fox and Raichle, 2007). In mTLE, abnormally decreased FC is also found within the mTL and DMN (Pittau et al., 2012).

Hubs of FC are brain regions presenting with denser connections than others (Bullmore and Sporns, 2009). Normal brain networks have a small-world topology characterized by many regular short-range connections and fewer long-range connections (Bassett and Bullmore, 2006). Connector hubs participating in inter-network connectivity include long-range connections (Power et al., 2013). Being more vulnerable than other hubs, they are primary targets of several neurological disorders including epilepsy (Crossley et al., 2014). In TLE, hubs are disrupted, i.e. hub regions in healthy subjects become non-hubs in patients, changes that are associated with damage to long connections, resulting in the evolution from a small-world network into a regularized network (Stam, 2014). Such evolution was observed in EEG studies during seizures (ictal state) (Schindler et al., 2008) and interictal epileptic discharges, between seizures (Bartolomei et al., 2013). Hub disruptions during resting state (RS) were also seen using fMRI (Vaughan et al., 2016; Ridley et al., 2015). Reorganization of hubs in mTLE correlates with epilepsy duration (van Dellen et al., 2009) and seizure frequency (Douw et al., 2015), and impacts post-surgical outcome (Wilke et al., 2011; van Dellen et al., 2014; Liao et al., 2016). However, the patterns of connector hub reorganization specific to mTLE and the characterization of the epilepsy-related networks involving these abnormal connector hubs remain largely unknown.

We recently developed a “SParse general linear model (GLM)-based Analysis of Reliable k-hubness (SPARK)” to study connector hub organization using RS fMRI (Lee et al., 2016). SPARK handles the multicollinearity of RSNs based on a biologically plausible sparsity assumption that a voxel can be involved in more than one but not all RSNs (Karahanoğlu and Van De Ville, 2015). The multicollinearity of RSNs is difficult to detect using conventional methods based on Pearson's correlation. For example, if regions A and B are temporally correlated through a functional process while regions B and C are correlated through another functional process, Pearson's correlation may connect A and C. The hubness metrics in graph theory, such as degree centrality, count the number of all pairwise significant correlations between voxels, up to several thousands (Bullmore and Sporns, 2009). SPARK measures a discrete hubness by counting the number (k) of networks involved in each voxel that are often <10, and identifies which networks are actually associated to each connector hub. This allows analyzing more specifically the network architecture in individual connector hubs. The robustness and test-retest reliability of SPARK were already validated using realistic simulations and real data (Lee et al., 2016).

We hypothesized that RSN reorganization occurs in mTLE with: (i) hub disruption where brain regions corresponding to connector hubs in healthy subjects, become non-hubs in patients, and (ii) hub emergence where brain regions corresponding to non-hub in healthy subjects, become connector hubs in patients. We focused our analysis on the mTL and DMN, as the involvement of these regions was suggested by other studies (Pittau et al., 2012; Vaughan et al., 2016). We estimated the connector hub organizations from RS fMRI acquired from 14 right (R)-lateralized mTLE and 6 left (L)-lateralized mTLE patients, and compared with the connector hub organization found in 17 healthy controls. To quantify the local reorganization of connector hubs within each volume of interest (VOI), we formulated a linear regression model to estimate the hub disruption index (HDI) and hub emergence index (HEI). We identified for each patient several epilepsy-related RSNs involving the pathological and contralateral hippocampi and assessed the lateralization of such epilepsy-related RSNs.

Material and methods

Subjects

Patients

From the database of drug-resistant patients with epilepsy who participated in EEG-fMRI studies at the Montreal Neurological Institute (MNI) between 2007 and 2016, we identified consecutive patients with mTLE based on the following inclusion criteria: (i) history of seizures with autonomic (i.e., epigastric sensation) or psychic features (e.g., déja vu or dreaming state) as the first ictal symptom, followed by arrest of activity, staring, and automatisms (mainly oroalimentary); (ii) interictal EEG epileptiform abnormalities over temporal or frontotemporal regions and ictal onset over the same regions; (iii) unilateral epileptogenic focus, which was lateralized based on clinical history, interictal and ictal video-EEG recordings, and neuroimaging; (iv) right-handed or left hemispheric dominance for language determined by neuropsychological evaluations and etomidate speech and memory test when indicated (Jones-Gotman et al., 2005).

We excluded some patients from whom we could not find one or more 6-min RS fMRI run showing (i) no indication of sleep, (ii) no ictal and less than three interictal epileptic discharges, (iii) no excessive motion during the 6-min scan. Specifically, for each patient, an expert epileptologist reviewed the simultaneous EEG recordings to select a 6-min run of RS fMRI which did not fulfill the exclusion criterion (i) and (ii). For the remaining runs, after applying motion correction of fMRI time series, we removed the time frames with excessive motion (frame displacement >0.5 mm) from each run by scrubbing (Power et al., 2012). Then, we excluded some runs if >35% of the time-points were removed by scrubbing because we wanted to have a sufficient length of time-course for our analysis. We also excluded some runs that had translation higher than 1 mm or rotation higher than 3°. Finally, we analyzed 6 patients with left (L) mTLE (age: 24 ± 7, 1F/5M) and 14 with right (R) mTLE (age: 32.9 ± 14.5, 11F/4M). 5/6 L mTLE patients and 13/14 R mTLE had MRI abnormalities over unilateral mesial temporal structures (see Table 1 for patient demographics).

Table 1

Patients demography.

ID.	Age	Sex	Syndrome	Handedness	Aetiology (MRI)	Epilepsy duration (yr)	Age of seizure onset (yr)
R1	28	F	R mTLE	R	Enlarged RH	10	18
R2	39	F	R mTLE	R	MTS	16	23
R3	17	F	R mTLE	R	MTS	7	10
R4	18	M	R mTLE	R	Gliosis w/o atrophy	1	17
R5	50	M	R mTLE	R	MTS	50	0.5
R6	25	F	R mTLE	R	MTS	12	13
R7	56	M	R mTLE	R	MTS	11	45
R8	20	F	R mTLE	R	MTS	10	10
R9	17	F	R mTLE	R	MTS	12	5
R10	22	M	R mTLE	R	MTS	13	10
R11	39	F	R mTLE	R	MTS	38	1
R12	53	F	R mTLE	R	Non-lesional	37	16
R13	49	F	R mTLE	R	HS	39	10
R14	27	F	R mTLE	R	MTS	21	6
L1	21	F	L mTLE	R	MTS	4	17
L2	28	M	L mTLE	R	MTS	10	18
L3	19	M	L mTLE	R	MTS	5	14
L4	31	M	L mTLE	R	MTS	24	7
L5	30	M	L mTLE	R	Enlarged LH	2	28
L6	46	M	L mTLE	R	Non-lesional	15	31

F: female, M: male. L/R mTLE: left/right mesial temporal lobe epilepsy. LH/RH: left/right hippocampus, HS: hippocampal sclerosis, MTS: mesial temporal sclerosis.

Controls

We recruited 21 age-matched (two-sample t-test, p < .05) right-handed healthy controls (age: 33 ± 12.5, 13F/8M) for comparison. All control subjects had no history of any neurological or psychiatric disorders and no evidence of structural lesions on T1 weighted MRI. After fMRI preprocessing, four healthy subjects were excluded due to excessive motion evaluated by the same quality control as in patients (see the exclusion criteria (iii) for patients), resulting in 17 healthy controls for our connector hub analysis.

The study was approved by the MNI Research Ethics Board. Each subject gave written informed consent in accordance with the MNI and Hospital Research.

Simultaneous EEG/fMRI acquisition

EEG/fMRI acquisition, processing, and analysis were identical to that performed in previous studies (Khoo et al., 2017; Pittau et al., 2012). EEG was continuously recorded inside a 3 T MRI scanner (Siemens, Trio, Germany). EEG was acquired with 25 MR-compatible scalp electrodes (Ag/AgCl) using the 10–20 (reference FCz) and the 10–10 (F9, T9, P9 and F10, T10, and P10) systems. Two electrodes were placed on the shoulder to record the electrocardiogram. Data were transmitted from a BrainAmp amplifier (Brain Products, Munich, Germany; 5 kHz sampling rate) to the EEG monitor outside the scanner room via optic fiber. A T1-weighted anatomic acquisition was obtained (1 mm slice thickness, 256 × 256 matrix; echo time [TE] = 7.4 ms and repetition time [TR] = 23 ms; flip angle 30°) and used to superimpose functional images and for intersubject coregistration. The functional data were acquired in runs of 6 min using a T2⁎-weighted echo-planar imaging (EPI) sequence, with two acquisition protocols (protocol 1: 64 × 64 matrix; 25 slices, 5 × 5 × 5 mm, TE = 30 ms, TR = 1.7 s, or protocol 2: 64 × 64 matrix; 33 slices, 3.7 × 3.7 × 3.7 mm, TE = 25 ms, TR = 1.9 s; flip angle 90°; from protocol 1 to protocol 2 we increased the number of coil channels to improve signal quality). Patients were instructed not to move and to stay with eyes closed, resting. 10/14 R mTLE patients, 5/6 L mTLE patients, and all healthy controls were scanned using the protocol 2.

EEG preprocessing

The gradient artifacts were offline corrected using Brain Vision Analyzer software (Brain Products), and the remaining artifacts were removed by a 50 Hz low-pass filter. The ballistocardiographic artifact was removed by averaged artifact subtraction-based methods (Allen et al., 1998) or independent component analysis (Bénar et al., 2003). For patients, an expert epileptologist reviewed the EEG recording to mark the interictal epileptic discharges and excluded the runs showing patterns of sleep and more than three interictal epileptic discharges. As a result, one to four 6-min resting-state EEG-fMRI recordings were identified for each subject.

fMRI preprocessing

We used NIAK 0.13.0 for fMRI preprocessing (Bellec et al., 2010)1 with the Minc Tool Kit2. The EPI volumes were first corrected for inter-slice timing differences. Rigid-body motion (3 translations and 3 rotations) was estimated within and between sessions using the median volume of the first session as a target, by maximization of correlation coefficient using Minctracc (Collins et al., 1994). T1 anatomical MRI scans were preprocessed using CIVET pipeline (Ad-Dab'bagh et al., 2006), according to the preprocessing pipeline suggested in Zijdenbos et al. (2002). Linear coregistration of a T1 anatomical scan in native space to the MNI stereotaxic space was performed with a nine parameter transformation (3 translations, 3 rotations and 3 scales) by maximization of correlation coefficient using Minctracc (Collins et al., 1994). The target template brain model was the mni-models-icbm152-nl-2009-1.0 template (Fonov et al., 2011). To correct for local deformation in the volume and shape of brain areas, non-linear coregistration from T1 native to the MNI ICBM 152 template (Collins et al., 1995; Collins and Evans, 1997) was then performed by maximization of correlation coefficient embedded within a multiresolution iterative strategy. Non-uniformity correction (Sled et al., 1998) and brain extraction (Park and Lee, 2009) were applied twice first in native space and then in the MNI template after nonlinear coregistration. Individual coregistration from an EPI volume (the median volume of the first session as defined during the motion estimation) to T1-weighted volume in native space was performed for each subject using Minctracc. To do so, both anatomical and functional volumes were corrected to have a zero mean inside their brain masks, and values outside the masks were set to zero. Those volumes were then linear coregistered using a rigid-body transformation (3 translations and 3 rotations) by maximization of the mutual information between volumes (Maes et al., 1997; Collignon et al., 1995), and this procedure was refined in a multiresolution iterative strategy. Finally, the EPI volumes were resampled in the MNI space by combining a linear transformation from EPI to T1 and non-linear transformation from T1 to the MNI template at a voxel resolution of 4 × 4 × 4 mm3.

Then, we removed the time frames with excessive motion (frame displacement >.5 mm) by scrubbing (Power et al., 2012). Slow time drifts were filtered out with the basis of discrete cosines using a .01 Hz cutoff. Several additional confounds were also regressed out from the fMRI data to reduce the influence of spatially structured noise: motion parameters, the average signals in the white matter and the lateral ventricles. Specifically, a twelve parameters model for motion was built by including the six rigid-body motion parameters (3 translations and 3 rotations) as well as their square (Lund et al., 2006). Then, a principal component analysis was used to retain 95% of the variance of these twelve parameters. Finally, each fMRI run in the MNI stereotaxic space was spatially smoothed using an isotropic Gaussian blurring kernel with 8 mm FWHM, which might allow for an improvement of the signal-to-noise ratio, mitigating the potential impact of misregistration between functional areas across subjects. Several runs were excluded when >35% of the total time frames were removed by scrubbing, translation >1 mm or rotation >3°. Finally, only the fMRI run associated with a minimal motion was selected for each subject for the connector hub analysis.

For these runs, gray matter voxels were selected by coregistering and resampling a modified automated anatomical labeling (AAL) template to the average of the functional images in the MNI space (Tzourio-Mazoyer et al., 2002). By default, NIAK package provides a modified version of AAL template, because NIAK uses the different stereotaxic space for coregistration than the one used by Tzourio-Mazoyer et al. (2002). The modified AAL template was generated by NIAK using non-linear coregistration to the MNI ICBM 152 template (Fonov et al., 2011), using the same coregistration strategy applied for our preprocessing. Therefore, we resampled the modified AAL template at a voxel resolution of 4 × 4 × 4 mm3, the resolution used for resampling of our EPI images, and selected gray matter voxels belonging to the AAL mask. The segmentation of cortical regions defined by the modified AAL template was further used to identify several volume of interests in our analyses.

Individual functional connector hub analysis using SPARK

This section will overview the SPARK procedure to detect connector hubs from individual RS fMRI based on the sparse general linear model (GLM). A summary diagram of the SPARK procedure is illustrated in Fig. 1. After preprocessing of an individual RS fMRI (Y), the SPARK procedure was performed including the following steps. First of all, individual parameters, such as the total number N of expected networks, for the sparse GLM were estimated for the fMRI data Y. In AppendixA, we presented our new development for the automatic estimation of these individual parameters. Using a sparse dictionary learning algorithm for the sparse GLM, the whole brain fMRI data was decomposed into N resting state networks. By doing this, N temporal and spatial patterns of resting state networks were estimated from data. Using the sparse GLM, the BOLD signals in each voxel was modelled as a voxel-specific linear combination of k atoms, selected among N atoms from the data-driven dictionary. This means that “connector hubness” of each voxel was modelled by k resting state networks that are spatio-temporally overlapping and specifically contributing to the signal in this voxel. Details of the sparse GLM are described in Section 2.5.1. An important contribution of SPARK is that statistical reproducibility of the data-driven sparse GLM was estimated and only the reproducible components in our estimation will remain after thresholding. The assessment of statistical reproducibility in SPARK was based on a bootstrap resampling based strategy, which will be described in Section 2.5.2. Taking advantage of this assessment, we could finally analyze only individually reliable connector hubs in each subject.

Fig. 1

The procedure to detect reliable functional connector hub organization from individual resting state fMRI using SPARK. After preprocessing of an individual resting state fMRI (Y), individual parameters for the sparse general linear model (GLM) are estimated. A bootstrap resampling strategy preserving temporal structures in fMRI is applied to generate B resampled datasets (Y(, b = 1, …, B). Each resampled fMRI data is modelled using the sparse GLM in parallel: a dictionary of data-driven temporal features (() and the corresponding sparse coefficient matrix (X() are estimated for each Y(. The outputs of all resampled datasets are collected and a spatial clustering algorithm is applied on the spatial maps. N is the number of atoms in each dictionary and equivalent to the number of average spatial maps estimated by averaging the spatial maps in each cluster. After evaluating and removing statistically unreliable elements in the average spatial maps and noise-related atoms, only the signal atoms representing meaningful networks remain, i.e. N atoms. Finally, the SPARK procedure provides: (i) spatial maps of resting state networks, (ii) k-hubness in each voxel measured by counting the number of atoms (resting state networks) involved in each voxel, and (iii) the identification of which (k) networks are involved in each connector hub for an individual resting state fMRI, which are associated with k-hubness.

The sparse general linear model for SPARK

SPARK provided a multivariate FC hub analysis to identify connector hubs of overlapping brain networks, reliably from individual RS fMRI using a sparse GLM data analysis framework (Lee et al., 2011). The sparse GLM models the whole brain fMRI data Y ∈ ℝ (T time samples × V voxels) using a data-driven global dictionary Ω ∈ ℝ and a sparse coefficient matrix X ∈ ℝ with the corresponding noise E ∈ ℝ:

where y is a time-course measured in a voxel i and ∥x ∥0 indicates the L0 norm of a vector x, i.e. the number of non-zero elements. Each column of is an atom , i.e. a network-characteristic temporal feature that should be shared by a subset of voxels. Each column of X (x) is a sparse code for the voxel i, where the number of non-zeros in x defines the sparsity level k given the total number of atoms (N). The coefficients indicate the signal amplitudes of the atoms. This sparsity level k provides also information on how many networks the time course of a particular voxel i is involved with, therefore providing a discrete measurement of hubness: number of networks associated to this voxel. Therefore, each row of X provides the spatial map, which could be interpreted as a network of voxels sharing the similar time course in each atom. The BOLD signal y measured in the voxel i can be represented by a weighted linear combination of k atoms, indicating the k networks involved in the voxel where the whole brain activities can be represented by N (>k) networks that recruit different sets of voxels. We estimated and X that best described the data using a sparse dictionary learning algorithm, a variant of K-SVD algorithm (vK-SVD) (Aharon et al., 2006; Lee et al., 2016).

The vK-SVD algorithm considered for sparse GLM decomposition actually requires two parameters to be specified: an initial value of k (k) and the total number of atoms (N), i.e. the scale of networks characterizing the whole brain activities (see Fig. 1). In order to find the voxel-specific sparsity level k for each voxel i, the vK-SVD algorithm requires the initialization of sparsity level in all voxels. In this paper, we proposed three-steps procedure for the automatic estimation of the subject-specific scale of whole brain networks (N) and an initial value of k (k) for the sparse GLM using the minimum description length (MDL) criteria (Saito, 1994) (see AppendixA, Fig. A6). The MDL criteria is one of the model order selection methods suggested for sparse models of functional data.

The SPARK procedure to estimate reproducible hub structures

Using SPARK, the individual-level reproducibility of the estimated connector hub structures could be assessed by a bootstrap resampling-based strategy (Lee et al., 2016). First, we applied the circular block bootstrap resampling to the preprocessed RS fMRI data matrix (Y) and generated B = 200 resampled datasets with equal dimension (Y( ∈ ℝ, b = 1, …,B). The length of temporal blocks for each bootstrap sample was randomly chosen between 10 and 30 time samples to preserve the temporal structure in the BOLD signals (Bellec et al., 2010). Then, vK-SVD was applied on each resampled data (Y() in order to estimate a dictionary (() and the corresponding sparse representation (X(). The parallel process generated 200 dictionaries ( (of N time-courses) and 200 sparse coefficient matrices X( (of N spatial maps). The main objective of this parallel process was to find a spatially reproducible set of atoms across the 200 resampled datasets. C nearest neighbor spatial clustering was applied to the collection of all 200 × N spatial maps to identify C = N clusters and the maps were spatially averaged in each cluster.

Then, we removed inconsistent elements in the resultant N maps, i.e. low signal amplitudes after the averaging, assuming a Gaussian distribution of the background noise. Elements that were less reproducible, exhibiting an amplitude lower than the threshold, were considered as background noise. We chose a higher threshold (p < .01) for removing the background noise than the one used in our previous work using real data obtained from healthy controls (p < .05) (Lee et al., 2016). Considering that fMRI data obtained from patients might be less stable than healthy controls due to the presence of pathological brain activity, we used the higher threshold p < .01 to avoid a potential bias in our comparison between the patients and healthy controls. This resulted in some voxels with k = 0, indicating that the involvement of this voxel in those spatially reproducible RSNs was not significantly consistent, while we focused our RSN estimation on only highly reproducible structures. It however does not mean that this voxel was never involved in any functional brain network. We also classified several atom maps that included physiological noise such as remaining significant cardiac, respiratory artifacts or subject motion, and some redundant atoms that had <30 active voxels based on the manual classification criteria given by visual inspection (Griffanti et al., 2017). After all, it resulted in N(< N) reliable RSNs for each subject. Finally, k-hubness was quantified by counting the number of associated networks for each voxel, and connector hubs were identified as the voxels that had k-hubness >1.

A linear regression model for hub reorganization using k-hubness

After finishing the whole SPARK procedure for individual RS fMRI, we ended up with a k-map for each subject, in which the gray matter voxels have the estimated k-hubness value. We were interested in detecting local reorganization of connector hubs within a particular volume of interest (VOI) associated with the pathology in mTLE. We defined three VOIs: the L and R mTL regions - amygdala, hippocampus and parahippocampus, and the DMN - the anterior/posterior cingulate, medial prefrontal, precuneus, angular and mid temporal gyri, and temporal poles. Our definition of the DMN did not include parts of the mTL, which were included in the description of DMN in other studies (Andrews-Hanna et al., 2010; Thomas Yeo et al., 2011), because we wanted to study connector hub reorganization in the mTL specifically. These definitions of the mTL and DMN are used throughout this paper.

To quantify the connector hub reorganization in a patient in comparison to normal connector hubs estimated from the group of controls, we formulated a linear regression model for k-hubness that is similar to the one suggested for graph theory in Achard et al. (2012). Let k denote k-hubness at voxel i measured from a patient, μ and σ denote the mean and standard deviation of the k-hubness measured at voxel i over the group of controls. i is the index of a voxel belonging to a VOI. p or c denotes an individual patient or healthy control subject, and P or C denotes a group of patients or controls. We first estimated the reorganization of k-hubness in each voxel by subtracting the mean k over the healthy controls (μ) from the subject-specific k measured from a patient (k) and then dividing the difference by the standard deviation of k in controls (σ), as follows:

Local trend of the reorganization of k-hubness was estimated from the voxels belonging to each VOI, by plotting this normalized difference of subject-specific k (d) with respect to the normalized group average k in controls (μ/σ). We formulated this model as follows:

where a is a coefficient or slope and b is a y-intercept, assuming Gaussian distribution of k values. In Achard et al. (2012), they suggested drawing a scatter plot where x-axis denotes μ and y-axis denotes k − μ, whereas k indicated circ centrality in their study. The pattern of such scatter plots was shown to exhibit a linear relationship. Then, a slope of the linear fit should give us a measure of circ to which the connector hubs in healthy controls were disrupted in an individual patient. They suggested to define the slope as hub disruption index (HDI), i.e. HDI = a. Note that we proposed in Eqs. 2 and 3 to normalize the k values by the standard deviation of k in controls (σ). This normalization was not considered in previous studies when estimating HDI.

In this model, we actually could define a new metric of connector hub reorganization, the hub emergence index (HEI) using the intercept b, i.e. HEI = b. A negative coefficient a indicates that some connector hubs in controls become non-hub in the patient and a positive a indicates that k-hubness of the normal connector hubs become even higher in patients. Then, a positive intercept b would indicate the emergence of new connector hubs: some non-hubs in controls became connector hubs in the patient. The intercept b cannot have a negative value, because a voxel cannot be involved in a negative number of networks. a = 0 and b = 0 together would indicate that there was no difference in k-hubness between the patient and the group of controls. However, the parameter b has never been utilized in previous studies. For comparison purposes, HDI and HEI estimated for each of the individual healthy subjects (k) deviating from their group average were also computed.

Reorganization and lateralization of patient-specific networks involving the hippocampus

We further investigated the patient-specific network reorganization associated to the hippocampus by measuring the proportion of hippocampal volume involved in each network from the outputs of SPARK. SPARK can provide, for each voxel, the hubness level as well as the characterization of which networks are contributing to the signal of that voxel. We used such a property to estimate the volume proportion by counting the number of hippocampal voxels involved in each RSN j (n) with respect to the total number of voxels in the hippocampus (n).

The identification of the total volume of hippocampus in each preprocessed individual fMRI was performed in the MNI stereotaxic space using the AAL template (Tzourio-Mazoyer et al., 2002). This AAL template itself was also previously coregistered and normalized to the linear and nonlinear MNI space (the same MNI152 template) using the same method for fMRI preprocessing. It is noteworthy to mention that the linear and non-linear coregistration and normalization strategies considered in our study were also used successfully for the segmentation and volumetric analyses of the hippocampus, amygdala and other medial temporal lobe structures, for instance, to study pathological volume changes in patients with Alzheimer's disease (Zandifar et al., 2017) and developmental volume changes from childhood to adolescence (Hu et al., 2013). The AAL template was resampled to (4 mm)3 isotropic voxels. At this point, the EPI image and the AAL template were in the same space with the same resolution. In all patients, the hippocampus was delineated from the normalized template space using the AAL template. Doing so, our approach did not consider any possible volume loss of the pathological hippocampi (in case of hippocampal sclerosis), that was indeed partially compensated by the spatial normalization and AAL template measure. Assessing hub reorganization using SPARK as a function of patient specific volume loss in the hippocampus is a topic of great interest that will be considered in our future investigations.

Several major RSNs involved within the left or right hippocampus (LH/RH) were identified for each individual by selecting the networks for which >10% of the hippocampal volume was involved, i.e. %vol ≥ 10. These RSNs are the following: the temporal network - anterior, lateral, mesial, and posterior parts of the temporal lobes, the salience network - anterior cingulate and insula cortices, the subcortical network - basal ganglia, limbic system, thalamus, and hypothalamus, the visual network - primary and secondary visual area in the occipital lobes, the motor network - primary and secondary motor area in the pre/post-central gyri and supplementary motor area, and the cerebellar network involving the cerebellum. Note that our definition of the DMN in Section 2.6 includes the middle temporal gyri and temporal poles, which are also included in the definition of the temporal network. However, the definition of DMN had its main focus in the midline core components including the medial prefrontal, precuneus, posterior cingulate and angular gyri.

Second, we measured the lateralization index (LI) for each network j by computing the difference of the volume proportion between the ipsi- and contra-lateral hippocampi, as follows:

We also estimated the absolute lateralization index, i.e. ∣LI∣ =  ∣  % vol −  % vol ∣, because we wanted to study if there was a significant change in the overall amplitudes of lateralization indices measured in all types of hippocampal networks in patients, when compared to controls. Finally, to study more specifically which hippocampal networks were actually contributing to the lateralization, we drew a diagram for each patient showing the k-hubness map and all major hippocampal networks involved within the LH and RH.

Results

Using SPARK, the individual fMRI was decomposed into several temporal features, each being associated with a specific consistent RSN. During this decomposition process, the BOLD signal in each voxel was modelled as a weighted linear combination of a small number of features (sparse number k) mainly contributing to the signal, among all features characterizing the whole brain (Lee et al., 2011). A voxel-specific level of sparsity (k) indicated how many RSNs were involved in each voxel, therefore providing the k-hubness map. The total number of RSNs characterizing the whole brain activities and the voxel-specific sparsity levels were automatically estimated for each subject. The estimated total numbers of whole-brain RSNs were 19.7 ± 1.5 (mean ± standard deviation; SD) in L mTLE patients, 21.7 ± 3.6 in R mTLE patients, and 21.8 ± 2.5 in controls (Fig. A7). The distributions of the estimated numbers were not significantly different between groups (two-sample t-test). The reliability of SPARK decomposition at the single subject level was ensured using a bootstrap resampling strategy to keep only the most reproducible features (Bellec et al., 2010; Lee et al., 2016).

Disruption and emergence of connector hubs

To quantify the local connector hub reorganization within the mTL and DMN from each patient (p), we formulated a linear regression model for k-hubness, as suggested by Achard et al. (2012). We defined three VOIs - the L ad R mTL regions and the DMN. Using the linear regression model, we estimated the reorganization of k-hubness from the voxels belonging to each VOI in each patient, when compared to the k-hubness values from the same set of voxels estimated in the group of healthy controls. Two patterns of connector hub reorganization, disruption and emergence, were quantified by the slope (hub disruption index, HDI) and intercept (hub emergence index, HEI) of a straight line fitted to the scatter plot of our data, as described in Section 2.6.

Fig. 2 illustrates the estimation of local HDI and HEI in the L and R mTL regions from the groups of L and R mTLE patients when compared to the group of healthy controls. For each group of patients, we drew a scatter plot for k-hubness measured in the voxels belonging to each VOI. For group comparisons, we estimated d = (μ − μ)/σ, which is the normalized difference between the group average k-hubness in patients (μ) and controls (μ). The local HDI and HEI were estimated, at the group-level, using the slope and intercept of the estimated linear fitting line for each scatter plot for each group of patients. L mTLE patients exhibited abnormal hub reorganization in both ipsilateral (HDI: − 0.58, HEI: 0.28) and contralateral (HDI: − 0.53, HEI: .24) mTL. R mTLE patients showed an even greater abnormal hub reorganization in both ipsilateral (HDI: − .89, HEI: 1.12) and contralateral (HDI: − .74, HEI: .93) mTL, consisting of both hub disruption and hub emergence patterns.

Fig. 2

Estimation of local connector hub disruption and connector hub emergence indices using the linear regression model from the groups of L and R mTLE patients in comparison to healthy controls (CONT) in the left (a-b) and right (c-d) mesial temporal lobe regions. For each group of patients, we drew a scatter plot for k-hubness measured in the voxels belonging to each VOI. Each dot denotes a voxel in the gray matter belonging to each VOI. In each plot, x-axis denotes the normalized group average k-hubness in healthy controls (μ/σ) and y-axis denotes d = (μ − μ)/σ, which is the difference between the group average k-hubness in patients (μ) and controls (μ) normalized by the standard deviation (SD) of k in controls (σ). The local HDI and HEI were estimated using the slope (a) and intercept (b) of the estimated linear fitting line (in red) for each scatter plot.

We were interested in estimating the local HDI and HEI in patients at the individual level, taking advantage of SPARK for individual reliability. For comparison purposes, we estimated HDI and HEI on every healthy subject to assess their distributions within the control group. When estimating the HDI and HEI in each patient using d, we also observed a similar linear pattern of scatter plot in every patient, while the plot exhibited a discrete distribution of scattered dots due to the discrete nature of k-hubness (k). As a result, L mTLE patients exhibited significant abnormal hub disruptions in the ipsilateral mTL only (HDI: median = −.6; interquartile range, IQR = .6; two-tailed Wilcoxon rank sum test, p < .005), while for R mTLE patients, significant hub disruptions were found in ipsi- (HDI: median = −.8, IQR = 1; p < .001) and contralateral (HDI: median = −1, IQR = 1.2; p < .01) mTL (Fig. 3). The interquartile range (IQR) refers to the difference between the 75th and 25th percentiles of the values. L mTLE patients did not show significant hub emergence in the mTL when compared to controls, while, again, R mTLE patients showed highly significant hub emergence in both ipsi- (HEI: median = .7, IQR = 1.4; p < .001) and contralateral (HEI: median = .9, IQR = 1.6; p < .005) mTL. Within the DMN, both L and R mTLE patients exhibited abnormal hub disruptions. The median and IQR of HDI were − .3 and .1 in L mTLE (p < .005) and − .2 and 0.3 in R mTLE (p < .005). There was no significant hub emergence in L mTLE, while R mTLE showed a slightly significant emergence of new connector hubs (HEI: median = .3, IQR = .4; p < .05).

Fig. 3

Hub disruption index (HDI) and hub emergence index (HEI) were estimated by measuring the changes of the k-hubness values in each L and R mTLE patient (d) deviating from the group average k-hubness across all the healthy control subjects (CONT) in the voxels belonging to each VOI. We defined three VOIs: the L and R mTLs (left top corner of (a)), and the DMN (left top corner of (e)). Each point indicates the hub index estimated from a single subject. Bold/dotted lines in red: mean/median of the points. The points were laid over a 1.96 standard error of the mean (SEM) (95% confidence interval) in red and a 1 SD in blue. *p < .05, **p < .01, ***p < .005, ****p < .001, two-sided Wilcoxon rank sum test for equal medians.

To evaluate a potential effect of combining data acquired using two acquisition protocols on our hub analysis, we performed the two-tailed Wilcoxon rank sum tests using only the HDI/HEI values estimated from those 10 R mTLE patients, 5 L mTLE patients and 17 healthy controls who were scanned using Protocol 2. The tests were same as those using the entire datasets. As a result, both L and R mTLE patients showed significant connector hub disruption in the ipsilateral mTL region (L: p < .01, R: p < .05). We also found significant connector hub emergence in both the ipsi- (p < .005) and contra-lateral (p < .05) mTL regions in patients with R mTLE, where we did not find any mTL hub emergence in patients with L mTLE. Within the DMN, we found a significant connector hub disruption in both L mTLE (p < .005) and R mTLE (p < .05), while we did not find any significant connector hub emergence in both L and R mTLE. These results being similar than the ones obtained considering the complete database, we assumed that the change of acquisition protocol had overall little impact on the proposed connector hub analysis using SPARK.

Reorganization of the networks involving the hippocampus

For our patients and healthy controls, we found a set of functional RSNs similar to those typically found in the FC literature (Smith et al., 2009). Among those RSNs, we studied more specifically the RSN organizations associated to the L or R hippocampi (LH/RH) (see Section 2.7 for the regions involved in each RSN). We measured the proportion of hippocampal volume (%vol) involved in each RSN, by counting the number of voxels in the hippocampus associated with a specific network to the total number of voxels in each hippocampus. The total volume of each hippocampus in each preprocessed individual fMRI was defined using the AAL template in the same MNI space with the same voxel resolution, i.e. (4 mm)3. The total number of hippocampal voxels identified in the MNI template space was 115 for the LH and 120 for the RH. We identified some important “hippocampal networks”, RSNs involving >10% of the voxels of the hippocampus (%vol > 10%). The occurrence rate for each hippocampal network was computed by counting the number of occurrences of a specific hippocampal network in each population to the total number of all hippocampal networks occurred in this population. For the healthy controls, the LH or RH involved a median of 2 networks (IQR = 2): Temporal (occurrence rate: 47% of the population), salience (31%), DMN (10%), visual (6%), subcortical (3%), and cerebellar (3%) (Fig. 4(a)). For the L mTLE patients, each hippocampus involved a median of 2 networks (IQR = 1): Temporal (occurrence rate: 46%), salience (27%), motor (18%) and subcortical (9%). For the R mTLE patients, each hippocampus involved a median of 2 networks (IQR = 2): Temporal (occurrence rate: 46%), cerebellar (18%), subcortical (14%), salience (11%), default mode (7%), and visual (4%). The consistency of the subject-specific networks (atoms) across subjects was measured by the percentages shown in Fig. 4(a), while the circ of functional involvement of the hippocampus was different in each individual. We observed in some patients that parts of the temporal networks were actually mixed into the salience, subcortical and cerebellar networks, resulting in new patient specific RSNs identified by SPARK. Consequently, we classified such mixed networks as occurrence of the salience, subcortical and cerebellar networks respectively, when drawing the pie-plots presented in Fig. 4(a).

Fig. 4

Reorganization of the networks involving the hippocampus in mTLE. (a) Types of resting state networks involving the L or R hippocampus (LH/RH) from 17 healthy controls (CONT) and 6 L and 14 R mTLE patients. (b) The comparison of absolute lateralization indices estimated from all hippocampal networks presented in (a). Each dot indicates the laterality index estimated for each network. (c–d) The comparisons of lateralization indices specific for the temporal and salience networks involving the hippocampi. For patients, lateralization indices were estimated by the % volume difference between the ipsi- and contralateral hippocampi (%vol − %vol). For controls, (%vol − %vol). Bold/dotted lines in red: mean/median of the points. The points were laid over a 1.96 SEM (95% confidence interval) in red and a 1 SD in blue. **p < .01, ****p < .001, two-tailed Wilcoxon rank sum test.

Lateralization of the networks involving the hippocampus

Among the hippocampal networks shown in Fig. 4(a), 21/39(54%) networks found in patients involved both the LH and RH, whereas 22/32(69%) networks in controls involved both hippocampi. The remaining networks involved either the LH or RH. Moreover, in patients, the proportion of LH versus RH volumes involved in specific RSNs (%vol) was remarkably different even when both were involved in some networks. These results suggest decreased FC between the two hippocampi within a network, i.e. lateralization. We estimated the lateralization index (LI) for each network by the % volume difference between the ipsi- and contralateral hippocampi in patients (LI = %vol − %vol), while for controls, the LI = %vol − %vol. The absolute lateralization indices estimated for all hippocampal networks were significantly larger in patients (∣LI∣: median = 13.1, IQR = 16.9) than in controls (∣LI∣: median = 4.9, IQR = 6.8; p < .001, two-tailed Wilcoxon rank sum test) (Fig. 4 (b)). We also compared the LI estimated for each type of hippocampal networks. For the temporal networks, the SD of LI was larger and more variable in patients (Fig. 4(c)). The direction of laterality (ipsi > contra or ipsi < contra) was however inconsistent. When considering subjects for whom the salience network was functionally connected to the LH or RH, the LI estimated for the salience networks in 2 patients with L mTLE and 3 patients with R mTLE were significantly different from the LI estimated for 10 controls (p < .01 for both L/R mTLE, two-tailed Wilcoxon rank sum test). In one patient with L mTLE, we found two sub-networks of the salience network, resulting in two data points from this patient in Fig. 4(d). The direction of laterality was ipsi > contra in R mTLE and ipsi < contra in L mTLE (Fig. 4(d)). Other types of networks were also analyzed but are not reported because of a too small number of occurrences.

Patient-specific hippocampal network reorganization

We analyzed which networks specifically contributed to such lateralization observed with hippocampal FC in individual patients. Six examples from 4 patients with R mTLE (R1–4) and 2 with L mTLE (L1−2) are illustrated in Fig. 5. For R1, the pathological RH was more involved in the temporal (%vol, RH: 40% > LH: 16.5%) and salience (RH: 30.8% > LH: 7.8%) networks and less in the DMN (RH: 1.7% < LH: 32.2%) than the LH (Fig. 5(a)). For R2, the RH was less involved in the cerebellar network and the DMN (Fig. 5(b)). For R3, the RH was more involved in the salience network and functionally less connected to the subcortical and visual networks (Fig. 5(c)). For R4, the RH and LH were involved in the salience, subcortical, and temporal networks, however exhibiting less lateralization compared to other patients (Fig. 5(d)). For L1, the LH was functionally disconnected from the two sub-networks of salience function whereas the RH remained functionally connected (Fig. 5(e)). Finally, for L2, the LH was functionally disconnected from the salience network and functionally less connected to the temporal network than the RH (Fig. 5(f)).

Fig. 5

Hippocampal networks reorganization in patients with R or L mTLE. The k-hubness map estimated from each patient is shown at the bottom, and the circles indicate the ipsi- (red) and contralateral (black) hippocampi. Each line connects the hippocampus and the network involving this hippocampus, and the hippocampal volume proportion involved in each network is presented beside each line. Beside each k-hubness map, the mean and SD of k-hubness within the ipsi- and contra-lateral hippocampus (〈k〉ipsi, 〈k〉contra) are shown.

Clinical relevance

In our patients, duration of disease was not significantly different between the groups of L and R mTLE patients. We did not find a correlation between disease duration and the HDI and HEI. However, R mTLE patients exhibited a significant correlation between the age of seizure onset and connector hub disruption (the amplitude of HDI) (r = .6, p = .03) in the ipsilateral mTL. We did not find a significant correlation between the age of seizure onset and connector hub emergence (the amplitude of HEI) (r = .5, p = .08). We did not find any correlation in the contralateral mTL. L mTLE patients exhibited a significant correlation between the age of seizure onset and connector hub disruption in the contralateral mTL (r = .9, p = .02). We did not find a correlation between the age of seizure onset and connector hub emergence in L mTLE.

Discussion

In this study, we analyzed the patterns of connector hub reorganization in a group of drug-resistant mTLE patients. We hypothesized that RSN reorganization in mTLE includes normal connector hubs becoming non-hubs when involved in epilepsy-related RSNs and newly emerging connector hubs. We adopted a novel approach called SPARK (Lee et al., 2016) to detect and analyze connector hubs from RS fMRI. We estimated a hub disruption index (HDI) and defined a new metric of connector hub emergence (HEI). We also investigated, at the individual level, the lateralization of epileptic resting state networks involving the pathological hippocampus, taking advantage of the ability of SPARK to identify which networks are actually involved in each connector hub. We found significant pathological connector hub disruptions within the mTL and DMN. Also, our results suggest an asymmetrical connector hub reorganization in mTLE: Patients with R mTLE showed remarkably higher emergence of new connector hubs in the mTLs compared to those with L mTLE. Different patterns of lateralization of the salience network involving the pathological hippocampus were found in R vs L mTLE in those patients who exhibited functional involvement of the hippocampus in the salience network (i.e. 2 patients with L mTLE and 3 patients with R mTLE). The temporal, cerebellar, subcortical, default mode and motor networks also contributed to the lateralization of hippocampal networks. Finally, the patterns of connector hub reorganization indicate a regularization of the epilepsy-related RSNs, characterized by the loss of distant connections and emergence of local connections. Another important consequence of our findings is that we do not recommend flipping the brain images in the L-R direction for patients with mTLE in order to increase the size of the population, since one cannot assume similar patterns of FC in both hemispheres.

Specifically, connector hubs in the mTL regions were disrupted bilaterally in mTLE, although ipsilateral disruptions were slightly stronger than contralateral disruptions in L mTLE. Our results agree with a RS fMRI study using graph theory, which showed that R-sided epilepsy exhibited more connector hub disruption within the whole brain than L-sided epilepsy (Ridley et al., 2015). We however analyzed a more homogeneous group of patients with unilateral mTLE, and found a local connector hub disruption index (HDI) larger than the whole brain HDI in the group of patients they studied. Several explanations to the disruptions of mTL connector hubs can be suggested: (i) intra-hemispheric disconnection between the mTL and their neighboring temporal cortices (Bettus et al., 2009; Maccotta et al., 2013), (ii) inter-hemispheric disconnection between the ipsi- and contralateral temporal structures (Vaughan et al., 2016; Pittau et al., 2012), or (iii) disconnection between the mTL and distal regions, such as those involved in the DMN (Pittau et al., 2012).

An important finding of our study is that the emergence of new connector hubs in bilateral mTL was prominent only in R mTLE, while there was no emergence of mTL connector hubs in L mTLE when compared to healthy controls, suggesting it as a candidate biomarker for mTLE lateralization. An increased FC between the mTL structures or between the mTL and distributed networks can be at the origin of such hub emergence in R mTLE. In our subsequent analyses, the variability of RSNs associated to the LH or RH in individuals was lower in L mTLE than in R mTLE, when compared to healthy subjects (Fig. 4). A RS fMRI study found decreased ipsilateral FC within the mTL structures with increased FC in the non-epileptic side in mTLE (Bettus et al., 2009). A high-density EEG study reported more abnormal directed connections from ipsilateral to contralateral temporal lobes in R TLE than in L TLE, while the abnormal FC changes were associated with cognitive impairments only in R TLE (Coito et al., 2015). Such findings including our study raise the possibility of epileptic or compensatory reorganization of FC associated with the emergence of new connector hubs in the bilateral mTL specific to R mTLE.

Connector hubs within the DMN were significantly disrupted in both L and R mTLE, but R mTLE exhibited significantly more hub emergence in the DMN regions compared to L mTLE. Disruptions of connector hubs in the DMN are likely due to the loss of long-range connections, as connector hubs are often found in DMN regions of healthy subjects (Power et al., 2013). Decreases in hubness in mTLE was already described in parts of the DMN (Vaughan et al., 2016), and other studies suggested decreased FC within the DMN and between the mTL structures and the DMN (Pittau et al., 2012; Maccotta et al., 2013; Liao et al., 2011). Also, decreased FC between the anterior and posterior parts of DMN was associated with increased seizure frequency (Douw et al., 2015).

We observed clear reorganization of RSNs associated with the hippocampus in mTLE patients when compared to healthy controls (Fig. 4). The types of RSNs involving the hippocampus were less variable in L mTLE than in controls, and for patients with R mTLE, the types of hippocampal networks were similar to those found in controls. However, the occurrence rate of the salience, subcortical and cerebellar networks were clearly different. We also demonstrated significantly higher laterality index in hippocampal networks from patients than those from controls. These results suggest decreases in FC between the ipsi and contralateral hippocampi when involved in epilepsy-related RSNs. Our specific illustrations, taking advantage of the ability of SPARK to identify which networks are actually involved in each hub, allowed us to further characterize the lateralization of hippocampal RSNs (Fig. 5). In a sense, hub disruptions and emergence found in the mTL and DMN could be further understood by the trends suggested by these detailed analyses, showing FC of hippocampal sub-regions involved in different epilepsy-related RSNs.

Even though we could study the laterality of the salience network reorganization only in a small number of patients (3 R and 2 L mTLE patients, see Fig. 5), the direction of this reorganization was clearly distinguishable between the L and R mTLE. In healthy subjects, the salience network includes both insula and anterior cingulate cortex. It functions to segregate the most relevant among internal and external stimuli, switching between other networks to facilitate access to attention (Menon and Uddin, 2010). We observed that the ipsilateral hippocampus is more functionally connected to the salience network in R mTLE than the contralateral one, with an opposite pattern in L mTLE (Fig. 4(d)). These preliminary results suggest different mechanisms involving the salience network in R versus L mTLE. Maccotta and colleagues (Maccotta et al., 2013) showed that the affected temporal lobe has increased FC with the ipsilateral insula and neighboring subcortical regions, suggesting a facilitated pathway of seizure propagation involving the insula (Isnard et al., 2000). The involvement of the insula in the epilepsy-related RSN has been also suggested to participate to the concept of temporal plus epilepsy, which is characterized by the involvement of neighboring regions in addition to a temporal epileptic focus (Barba et al., 2017). Considering more bilateral hub abnormalities in R TLE than in L TLE, a possible scenario may be that the insula takes a role in the spread of FC abnormality to broader areas in R TLE. Such scenario may be associated with the more frequent impairment in attention in R TLE than in L TLE (Helmstaedter, 2013).

We observed that the hippocampal FC involved in the temporal networks was more lateralized in mTLE patients than in controls (Figs. 4(c) and 5). The direction of this laterality was however inconsistent in patients, and we were not able to propose a pattern to characterize L versus R mTLE. In some patients, parts of the temporal networks were mixed with the subcortical, cerebellar, and salience networks, suggesting enhanced FC between temporal and extra-temporal structures. We confirmed the involvement of DMN in the abnormal lateralization of hippocampal FC in mTLE. In two patients with R mTLE, the RH was less functionally connected with the DMN than the LH (Fig. 5(a,b)). These results are in agreement with previous work showing decreased FC between the hippocampi and between the ipsilateral hippocampus and DMN, as well as increased FC between the contralateral hippocampus and the DMN (Pittau et al., 2012; Voets et al., 2012; McCormick et al., 2013). The lateralization of hippocampal FC in the cerebellar network was also interesting, since it was mainly found in R mTLE patients (4/14 patients, Fig. 4(b), results not shown). Such lateralization of cerebellar networks from RS fMRI in mTLE might be explained by similar processes involved in crossed cerebellar diaschisis seen with interictal glucose hypometabolism and ictal blood flow (Nelissen et al., 2006).

Finally, we estimated if there was any correlation between patients' clinical information and the measures of connector hub disruption and emergence. We did not find a significant correlation between disease duration and connector hub disruption and emergence, whereas we found a significant correlation between the age of seizure onset and connector hub disruption in both L and R mTLE patients. Although we had a small sample size, it is possible that connector hubs may be particularly sensitive to disruption at an early age, when they are in development. This hypothesis should be further validated by studying a larger population of mTLE patients.

Our results using the HDI and HEI suggest that SPARK can be a promising tool to study hub reorganization in epilepsy and also in different neurological disorders. Indeed, our connector hub analyses in the hippocampus clearly showed that some voxels that were hubs in controls became non-hubs in patients and several other voxels that were non-hubs in controls became hubs in patients, suggesting different patterns of connector hub reorganization in sub-regions of the hippocampus. While many functional connectivity studies considered only the average signal within the whole hippocampus as a seed region or single node in graph theory, using either atlas-based methods or patient specific segmentations (Liao et al., 2010; Pereira et al., 2010), our analysis of the HDI and HEI using SPARK considered that functional reorganization associated with neuropathology might exhibit different patterns in sub-regions of the conventionally defined nodes in graph theory. Whereas in our proposed analysis, the same preprocessing was applied to all data (non-linear co-registration, spatial sampling), further characterization of sub-field of the hippocampus on such connector hub reorganization, obtained from patient specific segmentation, would be of great interest and will be considered in future investigations. For example, it would be interesting to study if the anterior versus posterior hippocampus in individual patients show different patterns of connector hub reorganization using SPARK with a higher spatial resolution (i.e. <2 mm) in individual patient's native space, as suggested from recent studies on the functional organization of the hippocampus (Strange et al., 2014; Abdallah et al., 2017). Whereas we assume that in our proposed analysis, any possible hippocampal volume loss was at least partially compensated during the nonlinear registration, it would be of great interest to correlate reorganization of connector hubs estimated using SPARK with structural abnormalities specific to mTLE (Voets et al., 2012), but this was out of the scope of present study.

It should be noted that we could not evaluate the effect of drugs on connector hubs due to the heterogeneity of antiepileptic drug therapy. It should be also discussed that the measure of k-hubness provides discrete k values that are often bounded to a narrow interval such as [0, 5]. Therefore, k-hubness values are not normally distributed, making standard voxel-wise statistical tests for k-hubness difficult. We already validated the test-retest reliability and spatial consistency of k-hubness using real data from 25 healthy controls (Lee et al., 2016). Future work with larger cohorts may enhance the statistical power to perform a voxel-wise group comparison and enable us to detect the spatial consistency of the changes in k-hubness in patients compared to healthy controls. From a clinical point of view, it would be interesting to correlate in individual patients connector hub disruption with a possible cognitive disruption involving multiple left hemisphere functions, and connector hub emergence with preserved cognitive function in the right hemisphere. Whereas we acknowledge that studying a more extended population should be considered in future studies, our proposed method SPARK together with HDI and HEI was able to find unique features characterizing L versus R mTLE, suggesting a possible candidate biomarker of mTLE at the single subject level.

Conclusion

In this study, we analyzed resting state functional MRI in 20 mTLE patients and 17 healthy controls using SPARK. SPARK provides a realistic measure of hubness, by counting the number (k) of networks involved in each voxel, with the unique ability to identify the networks associated to each connector hub. In mTLE, we found pathological connector hub disruptions within the mTL and default mode network. The emergence of new connector hubs in the mTL was prominent in right but not in left mTLE. Hippocampal connector hubs were differently reorganized between left and right mTLE, when involved in the salience network. This study suggests an asymmetrical connector hub reorganization and network regularization in mTLE, including the disruption of distant connections and the emergence of aberrant local connections.