Neural network oscillations are a fundamental mechanism for cognition, perception and consciousness. Consequently, perturbations of network activity play an important role in the pathophysiology of brain disorders. When structural information from non-invasive brain imaging is merged with mathematical modelling, then generative brain network models constitute personalized in silico platforms for the exploration of causal mechanisms of brain function and clinical hypothesis testing. We here demonstrate with the example of drug-resistant epilepsy that patient-specific virtual brain models derived from diffusion magnetic resonance imaging have sufficient predictive power to improve diagnosis and surgery outcome. In partial epilepsy, seizures originate in a local network, the so-called epileptogenic zone, before recruiting other close or distant brain regions. We create personalized large-scale brain networks for 15 patients and simulate the individual seizure propagation patterns. Model validation is performed against the presurgical stereotactic electroencephalography data and the standard-of-care clinical evaluation. We demonstrate that the individual brain models account for the patient seizure propagation patterns, explain the variability in postsurgical success, but do not reliably augment with the use of patient-specific connectivity. Our results show that connectome-based brain network models have the capacity to explain changes in the organization of brain activity as observed in some brain disorders, thus opening up avenues towards discovery of novel clinical interventions.

Introduction

Neural network oscillations are a fundamental mechanism for the establishment of precise spatiotemporal relationships between neural responses that are in turn relevant for cognition, memory, perception and consciousness. When neurons discharge, the subsequent oscillatory activity propagates through the network recruiting other brain regions, thereby dynamically binding widely distributed sets of neurons into functionally coherent ensembles, hypothesized to represent neural correlates of a cognitive or behavioural content (Singer, 1999). As the transient synchronization wave evolves, it establishes a spatiotemporal pattern characteristic for cognitive processes, sensory (Engel et al., 1991), motor and sensorimotor tasks (Kelso, 1995), resting state (Allen et al., 2014; Hansen et al., 2014) and stimulation paradigms (Spiegler et al., 2016). Alterations of the spatiotemporal organization of these network oscillations play an important role in the pathophysiology of brain disorders. In particular focal epilepsy shows complex alterations of resting state functional connectivity patterns, responsible for some of the symptomatology (Bettus et al., 2009; Liao et al., 2011; Haneef et al., 2012; Pittau et al., 2012; Ridley et al., 2015). Simple activation paradigms lack the functional complexity to explain the richness of observed spatiotemporal behaviours linked to these brain dynamics (Bressler, 1995), leaving it essentially to network processes to explain the origin of the emergent functional and pathological spatiotemporal patterns. When structural connectivity is merged with mathematical modelling, then generative brain network models can be constructed allowing for the exploration of the causal network mechanisms of brain function and clinical hypothesis testing.

Large-scale brain network models (BNMs) emphasize the network character of the brain and bring dynamical properties to the structural data of individual brains. In BNMs, a network region is a neural mass model of lumped neuronal activity and is connected to other regions via a connectivity matrix representing white matter fibre tracts of the brain. This form of virtual brain modelling (Jirsa et al., 2002, 2010) exploits the explanatory power of measured network connectivity imposed as a constraint upon network dynamics and has provided important insights into the mechanisms underlying the emergence of the resting-state networks dynamics (Ghosh et al., 2008; Deco et al., 2011) of healthy subjects, stroke (Falcon et al., 2016) and schizophrenic patients (Cabral et al., 2013). So far, these studies have exploited generic or averaged connectomes to uncover basic principles of brain network functioning. What yet needs to be demonstrated, however, is the influence of individual structural variations of the connectome on the large-scale brain network dynamics of the models. The impact for personalized medicine would be substantial, allowing exploiting the predictive value with regard to the pathophysiology of brain disorders, and their associated abnormal brain imaging patterns. A personalized BNM derived from non-invasive structural imaging data, potentially fit to non-invasive functional imaging data, would allow testing of clinical hypotheses and exploration of novel therapeutic approaches. To explore this capacity of personalized BNMs to serve as a clinical validation and exploration tool in brain network disorders, we here systematically test the virtual brain approach along the example of epilepsy. So far, neural mass models have proven successful in explaining the biophysical and dynamical nature of seizure onsets and offsets (Robinson et al., 2002; Wendling et al., 2002; Lopes da Silva et al., 2003; Breakspear et al., 2006; Kalitzin et al., 2010; Touboul et al., 2011; Kramer et al., 2012; Jirsa et al., 2014). Only recently, however, has propagation of epileptic seizures started to be studied using BNMs (Terry et al., 2012; Taylor et al., 2013; Hutchings et al., 2015). Partial seizures have been reported to propagate through large-scale networks in humans (Bartolomei et al., 2013) and animal models (Toyoda et al., 2013). Around 30% of the patients with focal epilepsies are drug-resistant. A possible treatment for these patients is the surgical resection of the epileptogenic zone, a localized region or network where seizures arise, before recruiting secondary networks, called the propagation zone (Talairach and Bancaud, 1966; Bartolomei et al., 2001; Spencer, 2002). As a part of the standard presurgical evaluation, stereotactic EEG is used to help correctly delineate the epileptogenic zone (Bartolomei et al., 2002). Alternative imaging techniques such as structural MRI, EEG, MEG and PET help the clinician to outline the epileptogenic zone. Recently, diffusion MRI and the derived streamlines reflecting the connectivity between different brain regions started to be investigated to better understand the pathophysiology of temporal lobe epilepsy, revealing reduced fractional anisotropy (Ahmadi et al., 2009; Bernhardt et al., 2013) or other structural alterations in the connectome of epileptic patients (Bonilha et al., 2012; Besson et al., 2014; DeSalvo et al., 2014). However, the usefulness of diffusion MRI in the delineation of the epileptogenic zone and propagation zone during presurgical evaluation of epilepsy remains elusive.

We sharpen our virtual brain hypothesis with regard to personalized large-scale BNMs as clinical tools for the case of epilepsy and will explore their predictive power with regard to the organization of the epileptogenic zone and the propagation zone prior to epilepsy surgery. In this article, we use connectivity matrices derived from patient-specific diffusion MRI scans to build BNMs for a cohort of 15 epileptic patients, simulate the individual seizure propagation patterns and validate the analytical and numerical predictions of the propagation zone against clinical diagnosis and stereotactic EEG signals. Our results demonstrate that personalized virtual brain models reliably predict the propagation zone for a given epileptogenic zone. A correlation of virtual brain-based simulations and surgical outcomes further underlines the predictive power of this approach.

Materials and methods

Patient selection and data acquisition

We selected 15 drug-resistant patients (seven males, mean age 33.4 years, range 22–56) with different types of partial epilepsy accounting for different epileptogenic zone localizations. Patients were retrospectively enrolled based on the following criteria: (i) stereotactic EEG recordings were performed and the epileptogenic zone was clearly defined from ictal recordings; and (ii) all these patients had available MRI with diffusion MRI sequences. All patients underwent a presurgical evaluation (

). The first phase in the evaluation of each patient is not invasive and comprises the patient clinical record, neurological examinations, PET, and EEG along with video monitoring. T1-weighted anatomical images (MPRAGE sequence, repetition time = 1900 ms, echo time = 2.19 ms, 1.0 × 1.0 × 1.0 mm, 208 slices) and diffusion MRI images (DTI-MR sequence, angular gradient set of 64 directions, repetition time = 10.7 s, echo time = 95 ms, 2.0 × 2.0 × 2.0 mm, 70 slices, b-weighting of 1000 s/mm2) were also acquired on a Siemens Magnetom Verio 3 T MR-scanner. From the gathered data clinicians conclude potential epileptogenic zones. Further elaboration on the epileptogenic zone was done in the second phase, which is invasive and comprises the placement of stereotactic EEG electrodes in or close to the suspected regions. These electrodes have 10 to 15 contacts that are 1.5 mm apart. Each contact is 2 mm of length and 0.8 mm in diameter. The stereotactic EEG was recorded by a 128 channel DeltamedTM system using a 256 Hz sampling rate. The stereotactic EEG recordings were band-pass filtered between 0.16 and 97 Hz by a hardware filter. All the chosen patients showed seizures in the stereotactic EEG starting in one or several localized areas, i.e. the epileptogenic zone before recruiting distant regions: the propagation zone. The position of the electrodes was determined by performing a CT scan or an MRI after implanting the electrodes. Two-thirds of the patients were operated, in agreement with current rate of surgery after stereotactic EEG recordings.

Additionally, five healthy control subjects were studied to evaluate the predictive power of subject-specific connectivity derived from diffusion MRI data and used in personalized brain network models. All controls signed an informed consent form according to the rules of the local ethics committee [Comité de Protection des Personnes (CPP) Marseille 2] and underwent the same MRI protocol.

Data processing

To import structural and diffusion MRI data in The Virtual Brain, the data were processed using SCRIPTS (Proix et al., 2016). This processing pipeline makes use of various tools such as FreeSurfer (Fischl, 2012), FSL (Jenkinson et al., 2012), MRtrix3 (Tournier) and Remesher (Fuhrmann et al., 2010) to reconstruct the individual cortical surface and large-scale connectivity. The surface was reconstructed using 20 000 vertices. Cortical and volumetric parcellations were performed using the Desikan-Killiany atlas with 70 cortical regions and 17 subcortical regions (Desikan et al., 2006). Two additional parcellations were used to quantify the uncertainty in the exact size of the epileptogenic zone due to the sparse sampling issue of stereotactic EEG. To do so, we subdivided each cortical regions of the Desikan-Killiany atlas in two and four to obtain 157 and 297 regions, respectively (Zalesky et al., 2010). The diffusion data were corrected for eddy-currents and head motions using eddy-correct FSL functions. Fibre orientation estimation was performed with constrained spherical deconvolution (Tournier et al., 2007), and improved with anatomically constrained tractography (ACT; Smith et al., 2012). Tractography was performed using 5 × 106 streamlines and were corrected using spherical-deconvolution informed filtering of tractograms (SIFT; Smith et al., 2013) to obtain 2.5 × 106 streamlines. The connectivity matrix was obtained by summing track counts over each region of the parcellation, and normalized so that the maximum value of the connectivity matrix was one (connection density: 45.7%). No threshold was used to prune weaker edges; however, the BNM is sensitive to connection strength, thereby effectively discarding the effect of smaller weights.

The CT or MRI scan performed after electrode placements were aligned with the structural MRI recorded before the surgery using the FLIRT function of FSL, with six degrees of freedom and a mutual information cost function. Each contact surface was reconstructed and assigned to the region of the corresponding parcellations containing the most of the contact volume.

Definition of the epileptogenic zone and the propagation zone

The epileptogenic zone was evaluated by a clinical expert (F.B.) based on the patient comprehensive evaluation data gathered throughout the two-step procedure (non-invasive and invasive). Non-invasive data were used to define the hypothesis about the epileptogenic zone and to guide stereotactic EEG implantation and strategy. Invasive stereotactic EEG signals are used to define the epileptogenic zone by the mean of visual analysis (regions involved at seizure onset, generally characterized by low voltage fast activity) and by quantification of the Epileptogenicity Index (Bartolomei et al., 2008). The propagation zone was as a first approach defined by a clinical expert (F.B.) by visual analysis of the stereotactic EEG recordings (delayed, rhythmic modifications, low values of epileptogenicity index) (referred thereafter as ‘stereotactic EEG clinician estimation’ of the propagation zone). As a secondary alternative approach, the propagation zone was defined through quantification of stereotactic EEG signals (referred thereafter as ‘stereotactic EEG signal quantification’ of the propagation zone). In particular, for each patient, all seizures were isolated in the stereotactic EEG time series. The bipolar stereotactic EEG was considered (between pairs of electrode contacts) and filtered between 1 and 50 Hz using a Butterworth band-pass filter. The energy of the signal is defined as Ei=xi2[n], where xi[n] is the nth value of the ith channel of the stereotactic EEG time series. A contact was considered to be in the propagation zone if its signal energy was responsible for at least 30% of the maximum signal energy over all the contacts excluding the ones in the epileptogenic zone. The corresponding region was then assigned to the propagation zone.

Comparing the estimates of the propagation zone

Two different scores were computed to compare predicted propagation zone with the estimated propagation zone as describe above. The binary score S1 simply counts the accordance of each region found in the simulated propagation zone (ensemble PZS) and the predicted propagation zone (ensemble PZP):  

(1)
S1=i=1NPZS{1if RiPZP0if RiPZP
The distance score S2 quantitatively estimates the L1-norm between the normalized probability of the predicted propagation xiPZS and the strength of the stereotactic EEG signal energy (for clinician prediction, the strength was set to 1):  
(2)
S2=i=1NPZS(1|xiPZSEiPZP|)

Chance level to predict areas in the propagation zone

The chance level was the probability of obtaining k regions in the propagation zone by drawing randomly n regions from the set of all regions in the parcellations, multiplied by the score obtained for k regions, summed over all possibilities for k (Haken, 2000). This can be mathematically written as:  

(3)
k=1mscore(k)·P(X=k)=k=1mkm·(mk)·(Nmnk)(Nm)
with mΖ the number of regions that are in the clinical estimation of the propagation zone, n the number of regions drawn randomly from the parcellations, N the total number of regions in the parcellation.

Modelling

Brain network model

Based on recordings of epileptic seizures in different species, Jirsa et al. (2014) identified the dominant bifurcations involved at seizure onset and offset amid the theoretically 16 possible classes of bifurcation pairs for bursting activity predicted by dynamical system theory (Izhikevich, 2000). This consideration resulted in a phenomenological model, called the Epileptor. The model autonomously switches between interictal and ictal states because of a slow permittivity variable that could be related to tissue oxygenation (Suh et al., 2006), metabolism (Zhao et al., 2011), and extracellular levels of ions (Heinemann et al., 1986). The Epileptor is mathematically a set of two coupled oscillators linked together by the slow permittivity variable z. The two oscillators account for the fast discharges (variables x1 and y1) and spike and wave events (variables x2 and y2) observed in electrographic seizure recordings. The Epileptor model was also shown to reproduce refractory status epilepticus and depolarization block (El Houssaini et al., 2015). Using time-scale separation as well as evidence from stereotactic EEG recordings, Proix et al. (2014) suggested considering seizure recruitment among brain regions on a slow time scale and defined a permittivity-based coupling. Such a coupling function can be expressed as a linear difference coupling term, subsuming first order deviations from the homeostatic equilibrium of the slow permittivity variable, expressing perturbations of ion regulation in remote regions. Possible biophysical mechanisms for such coupling include potassium spatial buffering by glial cells (Amzica et al., 2002), extracellular diffusion of potassium (Durand et al., 2010), and increase of firing rate in remote regions leading to increase of remote extracellular potassium (Avoli et al., 2016). We used the same approach here in the general case of N coupled Epileptors, which reads for each Epileptor i:  

(3)
x˙1,i=y1,if1(x1,i,x2,i)zi+I1y˙1,i=15x1,i2y1,iz˙i=1τ0(4(x1,ix0,i)zij=1NKij·(x1,jx1,i))x˙2,i=y2,i+x2,ix2,i3+I2+0.002g(x1,i)0.3(zi3.5)y˙2,i=1τ2(y2,i+f2(x2,i))
where  
f1(x1,i,x2,i)={x1,i33x1,i2if x1,i<0(x2,i0.6(zi4)2)x1,iif x1,i0f2(x2,i)={0if x2,i<0.256(x2,i+0.25)ifx2,i0.25g(x1,i)=t0teγ(tτ)x1,i(τ)dτ
and τ0=2857;τ2=10;I1=3.1; I2=0.45; γ=0.01. Kij is the connection strength between Epileptor i and Epileptor j as given by the connectivity matrix. The degree of excitability of each Epileptor is represented by the value x0,i that we varied in this study. To simplify the interpretation, we define Δx0,i=x0,ix0C with x0C=2.1 the critical value of excitability. If Δx0,i>0, a brain region is epileptogenic and seizures are triggered autonomously. Otherwise, Δx0,i<0 and regions are in a healthy equilibrium state.

Two-dimensional reduction of the Epileptor

Taking advantage of time scale separation and focusing on the slower time scale, the five-dimensional Epileptor reduces to (Proix et al., 2014):  

(4)
{x˙1,i=x1,i32x1,i2+1zi+I1,iz˙i=1τ0(4(x1,ix0,i)zij=1NKij(x1,jx1,i))
for each Epileptor i with τ0 = 2857 and I1,i = 3.1.

One-dimensional reduction of the Epileptor

Using averaging methods to project the activity on the slow manifold (Haken, 1987), we can further reduce the Epileptor to a one-dimensional system (

):  
(5)
z˙i=1τ0(4x0,i+4F(zi)zij=1NKij(F(zj)F(zi)))
with F(z)=1/4(16/38z629.6/27)

Numerical simulations

The simulations were performed with The Virtual Brain (TVB; Sanz Leon et al., 2013) using a Heun integration scheme (time step: 0.04 ms). A zero mean white Gaussian noise with a variance of 0.0025 was added to the variables x2 and y2 to make time series reproduces more closely the stereotactic EEG recorded seizures. Two hundred and fifty-six time steps correspond to 1 s of physical time for realistic seizure durations. In Fig. 2, simulated signals are filtered with an order 5 bandpass Butterworth filter (0.16 Hz–97 Hz) to reproduce hardware filters used in stereotactic EEG data acquisition.

Model-based prediction of propagation zone

We predict the propagation zone as a function of the spatial location of the epileptogenic zone, the network model and its chosen parameterization via the excitability values Δx0,i, and the connectivity matrix. To do so, we estimated the propagation zone by identifying the dominating subnetworks involved in the transition toward the seizure state via a linear stability analysis (

). The centre manifold theorem in non-linear dynamical system theory predicts that the dominating subnetwork is closest to the onset of instability at the bifurcation point. To make use of this insight and determine the propagation direction as equivalent to the strongest direction of instability, we performed the linear stability analysis on a reduced 2D Epileptor model and computed numerically the corresponding Jacobian matrix. We confirmed the validity of our approach with respect to the simulated Epileptor system ().

Surrogate models

Scores for each surrogate model were obtained by replacing the Epileptor models in the BNM by the following surrogate models, while keeping the connectivity matrices and the spatial location of the epileptogenic zone unchanged.

Fast coupling

Epileptors in the BNM are coupled via a slow permittivity coupling function. We implemented a surrogate model with coupling on the fast time scale. The coupling function in Equation 4 was set to operate on the fast time scale, with opposite sign to keep the coupling function excitatory:  

(6)
{x˙1,i=x1,i32x1,i2+1zi+I1,i+j=1NKij(x1,jx1,i)z˙i=1τ0(4(x1,ix0,i)zi)

Time-scale separation

The Epileptor model assumes a time-scale separation, where the slow permittivity variable drives the transition between the interictal and the ictal state. We implemented a surrogate model with no time scale separation. The time-scale separation was suppressed by setting τ0 = 1 in Equation 4.

Generic saddle-node bifurcation

The seizure onset in the Epileptor model is represented by a saddle-node bifurcation. We also implemented a generic normal form of a saddle-node bifurcation  

x˙i=xi2+x0,i+2.1+j=1NKijx1,j
and validated to which degree it could predict the propagation zone.

Surrogate connectivities

Scores for each surrogate connectivity matrix were obtained by replacing the connectivity matrix of each patient by the surrogate connectivity matrix, while keeping the Epileptor models and the spatial location of the epileptogenic zone unchanged.

Connectivity from control subjects

Connectivity matrices were obtained for five control subjects. For each patient, scores were computed for a surrogate model using the connectivity matrix of the individual control subjects, and then averaged over all the control subjects.

Shuffled connectivity

Connectivity matrices were shuffled while preserving the degree and weight distributions (Rubinov and Sporns, 2011). In short, the procedure iteratively rewires the connectivity matrix links by randomly swapping connected node pairs, and then attributes the weights to the new network while attempting to preserve each node strength, i.e. the sum of the weighted connections to each node. Five shuffled connectivity matrices were obtained.

Changing the weights

Each non-zero weight kij of the connectivity was summed to a draw of a uniform random distribution, whose values were taken between ɛ·kij and +ɛ·kij, with ɛ a chosen percentage.

Log of the connectivity matrix

We simply redefined kij=log(kij+1) as the log of the connectivity matrix.

Statistics

We compared the different group medians for different conditions using first a Kruskal-Wallis test, followed by pairwise Mann Whitney U-tests.

Results

Partial seizures were recorded with stereotactic EEG electrodes in 15 drug-resistant epileptic patients undergoing presurgical evaluation. The clinical characteristics of each patient are given in

. Two examples of epileptic seizures of Patient CJ are shown in Fig. 1A. On the left of Fig. 1A the seizure is symptomatic and propagates from the epileptogenic zone, i.e. a part of the left lateral occipital cortex (channel highlighted in yellow), to the propagation zone (channels highlighted in red). On the right of Fig. 1A the seizure is asymptomatic and stays limited to the epileptogenic zone. Structural and diffusion MRI were preprocessed to construct individualized virtual patients, comprising cortical surface, surface and volumetric parcellation, electrode positions, and structural connectivity matrix (Fig. 1B and C). We used three different parcellation scales with 17 subcortical regions and 70, 140 or 280 cortical regions, accounting for the uncertainty in the exact size of the epileptogenic zone due to the sampling issue of stereotactic EEG. These virtualized patients were then imported into The Virtual Brain (Sanz Leon et al., 2013), a neuroinformatics platform for large-scale brain simulation.
Figure 1

Stereotactic EEG data and reconstruction of virtual Patient CJ. (A) Two examples of partial seizures recorded with stereotactic EEG in this patient. Left, the seizure propagates from the epileptogenic zone (yellow) to the propagation zone (red). Right, the seizure is limited to the epileptogenic zone. The normalized energy of each channel is shown in colour on the side, with below the corresponding colour bar. Spatiotemporal activation patterns are shown at different time points of the seizures. lSPC = left superior parietal cortex; lIPC = left inferior parietal cortex; lLgG = left lingual gyrus; lLOC1 = part 1 of the left lateral occipital cortex; lFuG = left fusiform gyrus; lLOC2 = part 2 of the left lateral occipital cortex; lITG = left inferior temporal gyrus. (B) Coregistration of the T1 MRI (levels of grey), the parcellation with 157 regions (colours) and the intracranial electrodes (red strips). (C) Connectivity matrix obtained from diffusion MRI for this parcellation. Ls = left subcortical; Rs = right subcortical.

Figure 1

Stereotactic EEG data and reconstruction of virtual Patient CJ. (A) Two examples of partial seizures recorded with stereotactic EEG in this patient. Left, the seizure propagates from the epileptogenic zone (yellow) to the propagation zone (red). Right, the seizure is limited to the epileptogenic zone. The normalized energy of each channel is shown in colour on the side, with below the corresponding colour bar. Spatiotemporal activation patterns are shown at different time points of the seizures. lSPC = left superior parietal cortex; lIPC = left inferior parietal cortex; lLgG = left lingual gyrus; lLOC1 = part 1 of the left lateral occipital cortex; lFuG = left fusiform gyrus; lLOC2 = part 2 of the left lateral occipital cortex; lITG = left inferior temporal gyrus. (B) Coregistration of the T1 MRI (levels of grey), the parcellation with 157 regions (colours) and the intracranial electrodes (red strips). (C) Connectivity matrix obtained from diffusion MRI for this parcellation. Ls = left subcortical; Rs = right subcortical.

Modelling seizure propagation

A BNM (Sanz-Leon et al., 2015) was constructed by placing at each node of the parcellation a neural mass model able to reproduce the temporal seizure dynamics and to switch autonomously between interictal and ictal states, the so-called Epileptor model (Jirsa et al., 2014). Epileptors were connected together via a permittivity coupling acting on a slow time-scale (i.e. seconds), which is sufficient to describe the recruitment of other brain regions in the seizure (Proix et al., 2014). For each patient, the epileptogenic zone localization was evaluated by a clinical expert (F.B.) using patient presurgical evaluation data, and the corresponding regions were set with a high excitability value such that the Epileptors trigger seizures autonomously (Fig. 2A, regions in yellow). Figure 2B shows an example of a simulation for Patient CJ, reproducing both symptomatic and asymptomatic seizures without any parameter changes. For this simulation, the brain regions in the propagation zone were set with different excitability values to correctly reproduce this recruitment scenario. However, choosing these parameters is a difficult and computationally costly task. Here we presented a systematic method to estimate the propagation zone directly from the knowledge of the epileptogenic zone and the large-scale connectivity matrix (see ‘Materials and methods’ section).

Figure 2

Simulations of the BNM for the Patient CJ. (A) A network of Epileptor models is build using the connectivity matrix. The nodes in the epileptogenic zone are epileptogenic (Δx0,i > 0, in yellow), the nodes in the propagation zone have different excitability values (0.5 > Δx0,i> −0.5, shades of red), while all the other nodes are not epileptogenic (stable state, Δx0,i< 0, in shades of blue). The blue links represent the anatomical links of the connectivity matrix. (B) Example of time series generated by the simulated BNM with the connectome of Patient CJ. Without changing any parameters, the propagation zone is not always recruited, reproducing the two seizures types of this patient as shown in Fig. 1A.

Figure 2

Simulations of the BNM for the Patient CJ. (A) A network of Epileptor models is build using the connectivity matrix. The nodes in the epileptogenic zone are epileptogenic (Δx0,i > 0, in yellow), the nodes in the propagation zone have different excitability values (0.5 > Δx0,i> −0.5, shades of red), while all the other nodes are not epileptogenic (stable state, Δx0,i< 0, in shades of blue). The blue links represent the anatomical links of the connectivity matrix. (B) Example of time series generated by the simulated BNM with the connectome of Patient CJ. Without changing any parameters, the propagation zone is not always recruited, reproducing the two seizures types of this patient as shown in Fig. 1A.

We systematically applied our method to predict seizure propagation for the 15 patients. We set the regions in the epileptogenic zone close to the epileptogenic state (Δx0,j = −2.2, Figs 3A, 4A and 5A, regions in yellow), and we set all the other regions not epileptogenic with the same excitability (Δx0,i = −2.5, Figs 3A, 4A and 5A, regions in blue). Examples of propagation zones predicted by the analytical analysis for Patients CJ, GC, and ML are shown in Figs 3B, 4B and 5B, respectively. The model prediction of the propagation zone was compared to (i) the clinical estimation of the propagation zone based stereotactic EEG acquired during the patient evaluation (Figs 3C, 4C and 5C); and (ii) the total energy of the stereotactic EEG signal over each channel during a seizure to evaluate the propagation zone (Litt et al., 2001) (Figs 3D, 4D and 5D). For the comparison, two different scores were applied to measure the accuracy of the predicted propagation zone: (i) a normalized binary score; and (ii) a normalized distance between the predicted propagation zone values and the reference propagation zone (i.e. stereotactic EEG clinician estimation or stereotactic EEG signal quantification). In addition, we used three different parcellations to explore the influence of different sizes of the epileptogenic zone around the stereotactic EEG contacts considered to be in the epileptogenic zone. Figure 6 shows individual scores obtained for both reference measures with a normalized binary score and a parcellation of 158 regions for Patients CJ, GC, and ML. The average results across all patients and all parcellation types are given in

. The individual results are given in and the abbreviations taxonomy in . In each case, the chance level was found by computing the score for randomly selected regions (see ‘Materials and methods’ section), and is shown as a dashed line in Fig. 6 and .
Figure 3

Prediction of the propagation zone by linear stability analysis for Patient CJ. (A) A network of Epileptor models is built using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 3

Prediction of the propagation zone by linear stability analysis for Patient CJ. (A) A network of Epileptor models is built using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 4

Prediction of the propagation zone by linear stability analysis for Patient GC. A network of Epileptor models is built using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 4

Prediction of the propagation zone by linear stability analysis for Patient GC. A network of Epileptor models is built using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 5

Prediction of the propagation zone by linear stability analysis for Patient ML. A network of Epileptor models is build using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 5

Prediction of the propagation zone by linear stability analysis for Patient ML. A network of Epileptor models is build using the connectivity matrix. The nodes in the epileptogenic zone (EZ) are epileptogenic (Δx0,i< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,i< −0.5, in blue). Examples of the localization of the epileptogenic zone (yellow) and the propagation zone (PZ) (shade of red, as indicated by the colour bar) in Patient CJ such as found by: (B) linear stability analysis using the patient connectome; (C) stereotactic EEG clinician estimations; (D) stereotactic EEG signal quantifications. The propagation zone values are all the same for the stereotactic EEG clinician estimations. SEEG = stereotactic EEG.

Figure 6

Prediction scores for patient, control subject and shuffled connectivity matrices. Results for Patients CJ, GC and ML, compared to five control subjects and to shuffled connectivity matrix, for propagation zone (PZ) location according to (A) stereotactic EEG clinician estimation, and (B) stereotactic EEG signal quantification. The dashed line indicates the level of chance. In box plots, boxes represent the 25th and 75th percentiles, centre line indicates the median and the whiskers extend to the most extreme data points not considered outliers. SEEG = stereotactic EEG.

Figure 6

Prediction scores for patient, control subject and shuffled connectivity matrices. Results for Patients CJ, GC and ML, compared to five control subjects and to shuffled connectivity matrix, for propagation zone (PZ) location according to (A) stereotactic EEG clinician estimation, and (B) stereotactic EEG signal quantification. The dashed line indicates the level of chance. In box plots, boxes represent the 25th and 75th percentiles, centre line indicates the median and the whiskers extend to the most extreme data points not considered outliers. SEEG = stereotactic EEG.

To compare our results with the surgical outcome, for each patient, we estimated the number of regions, which were found in the propagation zone by our model and not explored by stereotactic EEG. These are the regions that were not taken into account by the clinicians in the presurgical evaluation (regions in green in Fig. 7A). We found that a large extent of this propagation zone was significantly correlated with poor seizure prognosis according to the Engel classification (Engel, 1993), classifying postoperative outcomes for epilepsy surgery (Kruskal-Wallis P = 0.064, Mann-Whitney U-test class I versus class III P < 0.05, class II versus class III P < 0.1, Fig. 4B using clinical prediction, 158 regions and the distance score).

Figure 7

Prediction of the surgical outcome. (A) Example for Patient CJ showing stereotactic EEG signal quantification of the propagation zone (PZ) overlaid in green by the regions found in the propagation zone by the linear stability analysis but not explored by stereotactic EEG, therefore not considered in the propagation zone by the clinical expertise. (B) Comparison for all patients of the size of the unexplored regions predicted as in the propagation zone by the analytical model and the Engel classification (stereotactic EEG clinician estimation, 158 regions, distance score). Significant P-values are indicated. **P < 0.05, *P < 0.1; Mann-Whitney U-test.

Figure 7

Prediction of the surgical outcome. (A) Example for Patient CJ showing stereotactic EEG signal quantification of the propagation zone (PZ) overlaid in green by the regions found in the propagation zone by the linear stability analysis but not explored by stereotactic EEG, therefore not considered in the propagation zone by the clinical expertise. (B) Comparison for all patients of the size of the unexplored regions predicted as in the propagation zone by the analytical model and the Engel classification (stereotactic EEG clinician estimation, 158 regions, distance score). Significant P-values are indicated. **P < 0.05, *P < 0.1; Mann-Whitney U-test.

Surrogate connectivities

To gain more confidence in our results we computed the following surrogate connectivities: structural and diffusion MRI of five different control subjects were preprocessed to construct five different control connectivity matrices. Using each patient-specific epileptogenic zone, we computed iteratively the average score obtained by these five control subject connectivity matrices. Scores obtained for the individualized connectivity matrices were higher than control subject connectivities for nine patients (Fig. 6 and

), but the group level comparison did not show significant differences of the medians (Mann-Whitney U-test, , control boxplot). With the same procedure, we examined the effects of shuffling the weights of the connectivity matrix, i.e. rewiring differently the connectome by changing the topology of the network. The shuffling conserved the weight and degree distribution of the connectivity matrices (Rubinov and Sporns, 2011). The scores were always higher with the individualized connectivity matrices (Fig. 3E and ), and the group level comparison showed significant differences of the medians (Kruskal-Wallis P < 0.01, Mann-Whitney U-test P < 0.01 for patients versus shuffled, and controls versus shuffled, , shuffle boxplot). Other connectivities [random Erdös-Renyi networks, Strogatz-Watts small-world networks (Watts and Strogatz, 1998) performed very poorly (results not shown)] compared to the patient connectivity matrices. We also examined the effect of randomly changing up to 20% and 40% of the values of the weights with respect to the original values, therefore respecting the topology of the network (, 20% and 40% boxplots). As expected the results did not change significantly. The weights show a truncated power-law distribution (). Other studies have used a rescaling of the distribution of weights by compressing the range of edge weights (Hagmann et al., 2007, 2008), but our results were degraded by taking the logarithm of the weight matrices (, log boxplot).

These results suggest that the topology of the connectivity matrix is significantly important to predict the recruited network.

Surrogate models

We also tested alternative models and coupling functions. First we show that recruited network, as predicted by the BNM, is well approximated by a function of the excitabilities and connection weights to the epileptogenic zone. The reduced Epileptor is a slow-fast system that we further reduced to a one-dimensional system by projecting the dynamics on the slow manifold. Since the permittivity coupling leads to small terms at the first order approximation, the fixed point solution of the system does not depend of the coupling function. Furthermore, the linear stability analysis can be simplified by assuming weak couplings (which is the case when normalizing the connectivity matrix to its maximum weight), and is shown to be only a function of the excitabilities and connection weights to the epileptogenic zone (see

). In the case of the linear stability analysis, since excitability of all the nodes outside the epileptogenic zone is equal, the connections weights of the different regions from the epileptogenic zone directly determine the propagation zone. Our analytical result is confirmed for a generic connectivity matrix () as well as for patient connectivity matrices (, topologic boxplot).

We checked the importance of our assumptions, i.e. time-scale separation, permittivity coupling and weak coupling, by using different surrogate models: without time scale separation, with fast coupling and with a normal form of a saddle-node. In each case, the predictions based on generic (

) or patient-specific (, fast, time-scale, and saddle-node boxplots) connectivity matrices were degraded. We, however, found that by normalizing our connectivity matrix to smaller weights leads to results comparable to our model, demonstrating that weak coupling is a crucial assumption for the prediction of recruited networks.

Discussion

The recruitment of distal brain regions in partial seizures is a large-scale phenomenon spanning multiple time scales. We predicted the recruitment of distant areas by seizures originating from a focal epileptogenic network by constructing large-scale BNMs based on individualized connectomes. We demonstrated that simulations and analytical solutions approximating the BNM behaviour predict the propagation zone as determined by stereotactic EEG recordings and clinical expertise. The robustness of the model predictions was examined by testing various surrogate models and connectivity matrices.

To predict the recruitment network we posited a slow permittivity difference coupling function (Proix et al., 2014), which approximates the effect of local and remote fast neuronal discharges as a perturbation of the slow permittivity variable from the local homeostatic equilibrium. Such mechanisms do not exclude additional couplings on faster time scales that could comprise spike-wave events synchronization. Synaptic and electric coupling alone fail, however, to explain the temporal delays up to tens of seconds long that can be observed using stereotactic EEG, electrocorticography, or microelectrode arrays during network recruitments (Bartolomei et al., 2008; Schevon et al., 2012; Martinet et al., 2015). The slow permittivity variable is supported by a variety of biophysical mechanisms that act on slow time scales such as tissue oxygenation (Suh et al., 2006), extracellular level of ions (Heinemann et al., 1986), and metabolism (Zhao et al., 2011), which were found to vary in mouse (Jirsa et al., 2014), cat (Moody et al., 1974) and baboon (Pumain et al., 1985) models of epileptic seizures.

Most computational models of seizure propagation focus on small continuous spatial scales (Ursino and La Cara, 2006; Kim et al., 2009; Hall and Kuhlmann, 2013) or population of neurons (Miles et al., 1988; Golomb and Amitai, 1997; Compte et al., 2003; Fröhlich et al., 2007; Bazhenov et al., 2008). To our best knowledge, modelling seizure propagation at a large-scale (i.e. based on diffusion MRI) was to date only used to investigate absence seizures (Taylor et al., 2013) or temporal lobe epilepsy (Hutchings et al., 2015). Other modelling studies focused on small networks to investigate the role of the topology and localization of the epileptogenic zone (Terry et al., 2012). We proposed the large-scale connectome to be a major determinant of recruitment networks, which can be investigated in large-scale brain models based on patient-specific data. Diffusion MRI has revealed a quantitative decrease of regional connectivity around the epileptogenic zone that is associated with a network reorganization (Ahmadi et al., 2009; Bonilha et al., 2012; Bernhardt et al., 2013; Besson et al., 2014; DeSalvo et al., 2014) and cognitive impairments (Leyden et al., 2015). Histological studies provide evidence of white matter alterations in temporal lobe epilepsy (Thom et al., 2001; Blanc et al., 2011). Functional, volumetric and electrographic data suggest a broad reorganization of the networks in epileptic patients (Lieb et al., 1987, 1991; Cassidy and Gale, 1998; Rosenberg et al., 2006; Bettus et al., 2009). Hemispherectomy (Beier and Rutka, 2013) or regional disconnection procedures (Mohamed et al., 2011), by disrupting tracts have shown the interest of cutting seizure propagation pathways. Altogether, this evidence highlights the large-scale character of partial seizure propagation in the human brain. On an ad hoc basis, clinicians routinely factor knowledge of white matter anatomy into the interpretation of stereotactic EEG data and seizure spread but no framework exists to probe the underlying mechanisms. In this article, we used patient-specific diffusion MRI data to systematically test the relevance of large-scale network modelling in predicting seizure recruitment networks.

Estimating the weights of the connectivity matrix by counting the streamlines among areas is controversial. The number of streamlines may not directly reflect the strength of signal transmission between two areas but simply the probability that a streamline between two regions is found by the tractography algorithm (Jbabdi and Johansen-Berg, 2011; Jones et al., 2013). However, using ACT (Smith et al., 2012) and SIFT (Smith et al., 2013) methods for the tractography have been shown recently to robustly quantify the total cross-sectional area of the white matter fibre bundles between areas and match white matter statistics estimated from post-mortem studies (Smith et al., 2015). Additionally, we have shown that varying the weights of the connectivity matrix has less influence on the recruited network than the topology of the connectivity matrix. As a consequence, tracks that have a large weight in the connectivity matrix are crucial in determining the recruitment network.

The epileptogenic zone and propagation zone estimation by the clinician is based not only on stereotactic EEG but also on prior knowledge gathered throughout the patient’s comprehensive presurgical evaluation (Bartolomei et al., 2002), and the final estimated epileptogenic zone may differ from the direct stereotactic EEG signal energy estimation. Several biomarkers are used in the daily clinical routine such as the Epileptogenicity Index (Bartolomei et al., 2010), which is based on the appearance of fast discharges (beta-gamma bands) and delays of recruitment, or electrical stimulations to explore tissue excitability (Valentín et al., 2005). Analyses of BNMs have shown that such biomarkers depend on the coupling assessing the threshold level separating ictal from interictal state (Jirsa et al., 2014; Proix et al., 2014). For this reason, any interpretation of these biomarkers as epileptogenicity must be performed with caution, keeping the network character of the brain in mind. The extent of the epileptogenic zone has been linked to the surgical prognosis but is difficult to estimate. Stereotactic EEG exploration uses only a limited number of electrodes (10 to 15) and therefore is sparse and can lead to poor surgical outcome if the epileptogenic zone is underestimated or misidentified. This study focused on validating the BNM methods as a tool to predict the spatial propagation of partial seizure in the human brain using personalized connectivity matrices. A limitation of this study for practical clinical use is that this method cannot directly estimate the epileptogenic zone, but rather validates an epileptogenic zone hypothesis by comparison of the predicted propagation zone with the clinician estimation of the propagation zone. Indirect estimates via data fitting of the BNM to brain imaging data have been obtained previously (Jirsa et al., 2016); or alternatively epileptogenic zone estimates may be obtained based on other metrics such as graph theoretical estimates (Burns et al., 2014). Here, by helping to predict the propagation zone based on the epileptogenic zone, our model can aid the clinician throughout the decision-making process regarding the localization of the epileptogenic regions in the brain, for instance by discarding epileptogenic zone hypotheses leading to spatial recruitment patterns in contradiction with stereotactic EEG or EEG recordings, hence improving stereotactic EEG electrode implantation and delineation of surgical resection limits. For example, if the propagation zone predicted by the BNM is larger than the propagation zone assessed by the clinical expert, larger extent of the epileptogenic zone and therefore of the surgical resection are worth considering (Fig. 4B). The contribution to epileptogenic zone of subcortical regions, which are linked to loss of consciousness and surgical prognosis (Arthuis et al., 2009), is notoriously difficult to estimate and can be quantified through the probability of recruitment of unexplored regions with the BNM. For generic connectivity, BNMs perform reasonably well in predicting the recruited network pattern as a first approximation, establishing the possibility to compute a catalogue of all possible propagation patterns, thus accelerating the presurgical evaluation process without resorting to a diffusion MRI scan. Predictions based on the individualized connectivity matrices were not significantly different at the group level from predictions based on control subjects’ connectivity matrices (

); however, we do observe a trend towards better scores for the individualized connectivity matrices. Diffusion MRI images were acquired using a sequence compatible with a clinical exam. However, although the sequence parameters used are still of good standard in the clinical context, recent advances in diffusion MRI acquisitions and processing (e.g. higher b-weighting values, tractography methods for multi-shell schemes, etc.) will help going towards better individualized predictions. We therefore suggest that for a detailed patient evaluation, the individual structural connectivity is essential for predicting seizure spatial propagation, paving the way for broader applications in personalized medicine (see Jirsa et al., 2016). More generally, several previous studies have demonstrated a link between anatomical networks and pathologies, in particular neurodegenerative diseases (Seeley et al., 2009). Therefore, our personalized BNM approach, as it is driven by structural data, establishes a novel way of thinking of how to study structure–function alterations in neurological and psychiatric disorders in general.

Supplementary material

is available at Brain online.

Funding

This work was supported by the Brain Network Recovery Group through the James S. McDonnell Foundation, FHU EPINEXT [A*MIDEX project (ANR-11-IDEX-0001-02) funded by the ‘Investissements d’Avenir’ French Government] and the European Union's Horizon 2020 research and innovation programme under grant agreement No. 720270.

Abbreviation

    Abbreviation
  • BNM

    brain network model

References

Ahmadi
ME
Hagler
DJ
McDonald
CR
Tecoma
ES
Iragui
VJ
Dale
AM
et al
.
Side matters: diffusion tensor imaging tractography in left and right temporal lobe epilepsy
.
AJNR Am J Neuroradiol
 
2009
;
30
:
1740
7
.
Allen
EA
Damaraju
E
Plis
SM
Erhardt
EB
Eichele
T
Calhoun
VD
.
Tracking whole-brain connectivity dynamics in the resting state
.
Cereb Cortex
 
2014
;
24
:
663
76
.
Amzica
F
Massimini
M
Manfridi
A
.
Spatial buffering during slow and paroxysmal sleep oscillations in cortical networks of glial cells in vivo
.
J Neurosci
 
2002
;
22
:
1042
53
.
Arthuis
M
Valton
L
Régis
J
Chauvel
P
Wendling
F
Naccache
L
et al
.
Impaired consciousness during temporal lobe seizures is related to increased long-distance cortical-subcortical synchronization
.
Brain
 
2009
;
132
:
2091
101
.
Avoli
M
De Curtis
M
Gnatkovsky
V
Gotman
J
Köhling
R
Lévesque
M
et al
.
Specific imbalance of excitatory/inhibitory signaling establishes seizure onset pattern in temporal lobe epilepsy
.
J Neurophysiol
 
2016
;
115
:
3229
37
.
jn.01128.2015
.
Bartolomei
F
Chauvel
P
Wendling
F
.
Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG
.
Brain
 
2008
;
131
:
1818
30
.
Bartolomei
F
Cosandier-Rimele
D
McGonigal
A
Aubert
S
Régis
J
Gavaret
M
et al
.
From mesial temporal lobe to temporoperisylvian seizures: a quantified study of temporal lobe seizure networks
.
Epilepsia
 
2010
;
51
:
2147
58
.
Bartolomei
F
Guye
M
Gavaret
M
Régis
J
Wendling
F
Raybaud
C
et al
.
The presurgical evaluation of epilepsies
.
Rev Neurol
 
2002
;
158
:
4S55
64
.
Bartolomei
F
Guye
M
Wendling
F
.
Abnormal binding and disruption in large scale networks involved in human partial seizures
.
EPJ Nonlinear Biomed Phys
 
2013
;
1
:
4
.
Bartolomei
F
Wendling
F
Bellanger
JJ
Régis
J
Chauvel
P
.
Neural networks involving the medial temporal structures in temporal lobe epilepsy
.
Clin Neurophysiol
 
2001
;
112
:
1746
60
.
Bazhenov
M
Timofeev
I
Fröhlich
F
Sejnowski
TJ
.
Cellular and network mechanisms of electrographic seizures
.
Drug Discov Today Dis Models
 
2008
;
5
:
45
57
.
Beier
AD
Rutka
JT
.
Hemispherectomy: historical review and recent technical advances
.
Neurosurg Focus
 
2013
;
34
:
E11
.
Bernhardt
BC
Hong
S
Bernasconi
A
Bernasconi
N
.
Imaging structural and functional brain networks in temporal lobe epilepsy
.
Front Hum Neurosci
 
2013
;
7
:
624
.
Besson
P
Dinkelacker
V
Valabregue
R
Thivard
L
Leclerc
X
Baulac
M
et al
.
Structural connectivity differences in left and right temporal lobe epilepsy
.
Neuroimage
 
2014
;
100
:
135
44
.
Bettus
G
Guedj
E
Joyeux
F
Confort-Gouny
S
Soulier
E
Laguitton
V
et al
.
Decreased basal fMRI functional connectivity in epileptogenic networks and contralateral compensatory mechanisms
.
Hum Brain Mapp
 
2009
;
30
:
1580
91
.
Blanc
F
Martinian
L
Liagkouras
I
Catarino
C
Sisodiya
SM
Thom
M
.
Investigation of widespread neocortical pathology associated with hippocampal sclerosis in epilepsy: a postmortem study
.
Epilepsia
 
2011
;
52
:
10
21
.
Bonilha
L
Nesland
T
Martz
GU
Joseph
JE
Spampinato
MV
Edwards
JC
et al
.
Medial temporal lobe epilepsy is associated with neuronal fibre loss and paradoxical increase in structural connectivity of limbic structures
.
J Neurol Neurosurg Psychiatry
 
2012
;
83
:
903
9
.
Breakspear
M
Roberts
JA
Terry
JR
Rodrigues
S
Mahant
N
Robinson
PA
.
A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis
.
Cereb Cortex
 
2006
;
16
:
1296
313
.
Bressler
SL
.
Large scale cortical networks and cognition
.
Brain Res Rev
 
1995
;
20
:
288
304
.
Burns
SP
Santaniello
S
Yaffe
RB
Jouny
CC
Crone
NE
Bergeyc
GK
et al
.
Network dynamics of the brain and influence of the epileptic seizure onset zone
.
Proc Natl Acad Sci USA
 
2014
;
111
:
E5321
30
.
Cabral
J
Fernandes
HM
Van Hartevelt
TJ
James
AC
Kringelbach
ML
Deco
G
.
Structural connectivity in schizophrenia and its impact on the dynamics of spontaneous functional networks
.
Chaos
 
2013
;
23
:
046111
.
Cassidy
RM
Gale
K
.
Mediodorsal thalamus plays a critical role in the development of limbic motor seizures
.
J Neurosci
 
1998
;
18
:
9002
9
.
Compte
A
Sanchez-Vives
MV
McCormick
DA
Wang
XJ
.
Cellular and network mechanisms of slow oscillatory activity (<1 Hz) and wave propagations in a cortical network model
.
J Neurophysiol
 
2003
;
89
:
2707
25
.
Deco
G
Jirsa
VK
McIntosh
AR
.
Emerging concepts for the dynamical organization of resting-state activity in the brain
.
Nat Rev Neurosci
 
2011
;
12
:
43
56
.
DeSalvo
MN
Douw
L
Tanaka
N
Reinsberger
C
Stufflebeam
SM
.
Altered structural connectome in temporal lobe epilepsy
.
Radiology
 
2014
;
270
:
842
8
.
Desikan
RS
Ségonne
F
Fischl
B
Quinn
BT
Dickerson
BC
Blacker
D
et al
.
An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest
.
Neuroimage
 
2006
;
31
:
968
80
.
Durand
DM
Park
EH
Jensen
AL
.
Potassium diffusive coupling in neural networks
.
Philos Trans R Soc Lond B Biol Sci
 
2010
;
365
:
2347
62
.
Engel
AK
Konig
P
Kreiter
AK
Singer
W
.
Interhemispheric synchronization of oscillatory neuronal responses in cat visual cortex
.
Science
 
1991
;
252
:
1177
9
.
Engel
J
.
Surgical treatment of the epilepsies
 .
Philadelphia, PA:
Lippincott Williams & Wilkins
;
1993
.
Falcon
MI
Riley
JD
Jirsa
V
McIntosh
AR
Elinor Chen
E
Solodkin
A
.
Functional mechanisms of recovery after chronic stroke: modeling with the virtual brain
.
Eneuro
 
2016
;
3
.
ENEURO.0158–15.2016
.
Fischl
B
.
FreeSurfer
.
Neuroimage
 
2012
;
62
:
774
81
.
Fröhlich
F
Bazhenov
M
Timofeev
I
Sejnowski
TJ
.
Maintenance and termination of neocortical oscillations by dynamic modulation of intrinsic and synaptic excitability
.
Thalamus Relat Syst
 
2007
;
3
:
147
.
Fuhrmann
S
Ackermann
J
Kalbe
T
Goesele
M
.
Direct resampling for isotropic surface remeshing
. In:
Vision, modeling and visualization
 .
Siegen, Germany:
The Eurographics Association
;
2010
.
Ghosh
A
Rho
Y
McIntosh
AR
Kötter
R
Jirsa
VK
.
Noise during rest enables the exploration of the brain’s dynamic repertoire
.
PLoS Comput Biol
 
2008
;
4
:
e1000196
.
Golomb
D
Amitai
Y
.
Propagating neuronal discharges in neocortical slices: computational and experimental study
.
J Neurophysiol
 
1997
;
78
:
1199
211
.
Hagmann
P
Cammoun
L
Gigandet
X
Meuli
R
Honey
CJ
Wedeen
VJ
et al
.
Mapping the structural core of human cerebral cortex
.
PLoS Biol
 
2008
;
6
:
e159
.
Hagmann
P
Kurant
M
Gigandet
X
Thiran
P
Wedeen
VJ
Meuli
R
et al
.
Mapping human whole-brain structural networks with diffusion MRI
.
PLoS One
 
2007
;
2
:
e597
.
Haken
H
.
Advanced synergetics: instability hierarchies of self-organizing systems and devices
 .
2nd edn
Berlin
:
Springer
;
1987
.
Haken
H
.
Information and self-organization
 
2nd edn
.
Springer Berlin Heidelberg
;
2000
.
Hall
D
Kuhlmann
L
.
Mechanisms of seizure propagation in 2-dimensional centre-surround recurrent networks
.
PLoS One
 
2013
;
8
:
e71369
.
Haneef
Z
Lenartowicz
A
Yeh
HJ
Engel
J
Stern
JM
.
Effect of lateralized temporal lobe epilepsy on the default mode network
.
Epilepsy Behav
 
2012
;
25
:
350
7
.
Hansen
ECA
Battaglia
D
Spiegler
A
Deco
G
Jirsa
VK
.
Functional connectivity dynamics: modeling the switching behavior of the resting state
.
Neuroimage
 
2014
;
105
:
525
35
.
Heinemann
U
Konnerth
A
Pumain
R
Wadman
WJ
.
Extracellular calcium and potassium concentration changes in chronic epileptic brain tissue
.
Adv Neurol
 
1986
;
44
:
641
61
.
El Houssaini
K
Ivanov
AI
Bernard
C
Jirsa
VK
.
Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities
.
Phys Rev E Stat Nonlin Soft Matter Phys
 
2015
;
91
:
010701
.
Hutchings
F
Han
CE
Keller
SS
Weber
B
Taylor
PN
Kaiser
M
.
Predicting surgery targets in temporal lobe epilepsy through structural connectome based simulations
.
PLoS Comput Biol
 
2015
;
11
:
e1004642
.
Izhikevich
EM
.
Neural excitability, spiking and bursting
.
Int J Bifurcat Chaos
 
2000
;
10
:
1171
266
.
Jbabdi
S
Johansen-Berg
H
.
Tractography: where do we go from here?
.
Brain Connect
 
2011
;
1
:
169
83
.
Jenkinson
M
Beckmann
CF
Behrens
TEJ
Woolrich
MW
Smith
SM
.
FSL
.
Neuroimage
 
2012
;
62
:
782
90
.
Jirsa
VK
Jantzen
KJ
Fuchs
A
Kelso
SJA
.
Spatiotemporal forward solution of the EEG and MEG using network modeling
.
IEEE Trans Med Imaging
 
2002
;
21
:
493
504
.
Jirsa
VK
Proix
T
Perdikis
D
Woodman
MM
Wang
H
Bénar
CG
et al
.
The virtual epileptic patient: individualized whole-brain models of epilepsy spread
.
Neuroimage
 
2016
;
145
:
377
88
.
Jirsa
VK
Sporns
O
Breakspear
M
Deco
G
McIntosh
RA
.
Towards the virtual brain: network modeling of the intact and the damaged brain
.
Arch Ital Biol
 
2010
;
148
:
189
205
.
Jirsa
VK
Stacey
WC
Quilichini
PP
Ivanov
AI
Bernard
C
.
On the nature of seizure dynamics
.
Brain
 
2014
;
137
:
2210
30
.
Jones
DK
Knösche
TR
Turner
R
.
White matter integrity, fiber count, and other fallacies: the do’s and don’ts of diffusion MRI
.
Neuroimage
 
2013
;
73
:
239
54
.
Kalitzin
SN
Velis
DN
Lopes Da Silva
FH
.
Stimulation-based anticipation and control of state transitions in the epileptic brain
.
Epilepsy Behav
 
2010
;
17
:
310
23
.
Kelso
SJA
.
Dynamic patterns: the self-organization of brain and behavior (Complex Adaptive Systems)
 .
Cambridge, MA
:
MIT Press
;
1995
.
Kim
JW
Roberts
JA
Robinson
PA
.
Dynamics of epileptic seizures: evolution, spreading, and suppression
.
J Theor Biol
 
2009
;
257
:
527
32
.
Kramer
MA
Truccolo
W
Eden
U
Lepage
KQK
Hochberg
LR
Eskandar
EN
et al
.
Human seizures self-terminate across spatial scales via a critical transition
.
Proc Natl Acad Sci USA
 
2012
;
109
:
21116
21
.
Leyden
KM
Kucukboyaci
NE
Puckett
OK
Lee
D
Loi
RQ
Mcdonald
CR
.
What does diffusion tensor imaging (DTI) tell us about cognitive networks in temporal lobe epilepsy?
.
Quant Imaging Med Surg
 
2015
;
5
:
247
63
.
Liao
W
Zhang
Z
Pan
Z
Mantini
D
Ding
J
Duan
X
et al
.
Default mode network abnormalities in mesial temporal lobe epilepsy: a study combining fMRI and DTI
.
Hum Brain Mapp
 
2011
;
32
:
883
95
.
Lieb
J
Hoque
K
Skomer
C
Song
X
.
Inter-hemispheric propagation of human mesial temporal lobe seizures: a coherence/phase analysis
.
Electroencephalogr Clin Neurophysiol
 
1987
;
67
:
101
19
.
Lieb
JP
Dasheiff
RM
Engel
J
.
Role of the frontal lobes in the propagation of mesial temporal lobe seizures
.
Epilepsia
 
1991
;
32
:
822
37
.
Litt
B
Esteller
R
Echauz
J
D’Alessandro
M
.
Epileptic seizures may begin hours in advance of clinical onset: a report of five patients
.
Neuron
 
2001
;
30
:
51
64
.
Lopes da Silva
FH
Blanes
W
Kalitzin
SN
Parra
J
Suffczynski
P
Velis
DN
.
Dynamical diseases of brain systems: different routes to epileptic seizures
.
IEEE Trans Biomed Eng
 
2003
;
50
:
540
8
.
Martinet
LE
Ahmed
OJ
Lepage
KQ
Cash
SS
Kramer
MA
.
Slow spatial recruitment of neocortex during secondarily generalized seizures and its relation to surgical outcome
.
J Neurosci
 
2015
;
35
:
9477
90
.
Miles
R
Traub
RD
Wong
RK
.
Spread of synchronous firing in longitudinal slices from the CA3 region of the hippocampus
.
J Neurophysiol
 
1988
;
60
:
1481
96
.
Mohamed
AR
Freeman
JL
Maixner
W
Bailey
CA
Wrennall
JA
Harvey
AS
.
Temporoparietooccipital disconnection in children with intractable epilepsy
.
J Neurosurg Pediatr
 
2011
;
7
:
660
70
.
Moody
WJ
Futamachi
KJ
Prince
DA
.
Extracellular potassium activity during epileptogenesis
.
Exp Neurol
 
1974
;
42
:
248
63
.
Pittau
F
Grova
C
Moeller
F
Dubeau
F
Gotman
J
.
Patterns of altered functional connectivity in mesial temporal lobe epilepsy
.
Epilepsia
 
2012
;
53
:
1013
23
.
Proix
T
Bartolomei
F
Chauvel
P
Bernard
C
Jirsa
VK
.
Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy
.
J Neurosci
 
2014
;
34
:
15009
21
.
Proix
T
Spiegler
A
Schirner
M
Rothmeier
S
Ritter
P
Jirsa
VK
.
How do parcellation size and short-range connectivity affect dynamics im large-scale brain network models?
.
Neuroimage
 
2016
;
142
:
135
49
.
Pumain
R
Menini
C
Heinemann
U
Louvel
J
Silva-Barrat
C
.
Chemical synaptic transmission is not necessary for epileptic seizures to persist in the baboon Papio papio
.
Exp Neurol
 
1985
;
89
:
250
8
.
Ridley
BGY
Rousseau
C
Wirsich
J
Le Troter
A
Soulier
E
Confort-Gouny
S
et al
.
Nodal approach reveals differential impact of lateralized focal epilepsies on hub reorganization
.
Neuroimage
 
2015
;
118
:
39
48
.
Robinson
P
Rennie
C
Rowe
D
.
Dynamics of large-scale brain activity in normal arousal states and epileptic seizures
.
Phys Rev E Stat Nonlin Soft Matter Phys
 
2002
;
65
:
041924
Rosenberg
DS
Mauguière
F
Demarquay
G
Ryvlin
P
Isnard
J
Fischer
C
et al
.
Involvement of medial pulvinar thalamic nucleus in human temporal lobe seizures
.
Epilepsia
 
2006
;
47
:
98
107
.
Rubinov
M
Sporns
O
.
Weight-conserving characterization of complex functional brain networks
.
Neuroimage
 
2011
;
56
:
2068
79
.
Sanz Leon
P
Knock
SA
Woodman
MM
Domide
L
Mersmann
J
McIntosh
AR
et al
.
The virtual brain: a simulator of primate brain network dynamics
.
Front Neuroinform
 
2013
;
7
:
10
.
Sanz-Leon
P
Knock
SA
Spiegler
A
Jirsa
VK
.
Mathematical framework for large-scale brain network modelling in the virtual brain
.
Neuroimage
 
2015
;
111
:
385
430
.
Schevon
CA
Weiss
SA
McKhann
G
Goodman
RR
Yuste
R
Emerson
RG
et al
.
Evidence of an inhibitory restraint of seizure activity in humans
.
Nat Commun
 
2012
;
3
:
1060
.
Seeley
WW
Crawford
RK
Zhou
J
Miller
BL
Greicius
MD
.
Neurodegenerative diseases target large-scale human brain networks
.
Neuron
 
2009
;
62
:
42
52
.
Singer
W
.
Neuronal synchrony: a versatile code for the definition of relations
.
Neuron
 
1999
;
24
:
49
65
.
Smith
RE
Tournier
JD
Calamante
F
Connelly
A
.
Anatomically-constrained tractography: improved diffusion MRI streamlines tractography through effective use of anatomical information
.
Neuroimage
 
2012
;
62
:
1924
38
.
Smith
RE
Tournier
JD
Calamante
F
Connelly
A
.
SIFT: spherical-deconvolution informed filtering of tractograms
.
Neuroimage
 
2013
;
67
:
298
312
.
Smith
RE
Tournier
JD
Calamante
F
Connelly
A
.
The effects of SIFT on the reproducibility and biological accuracy of the structural connectome Article
.
Neuroimage
 
2015
;
104
:
253
65
.
Spencer
SS
.
Neural networks in human epilepsy: evidence of and implications for treatment
.
Epilepsia
 
2002
;
43
:
219
27
.
Spiegler
A
Hansen
ECA
Bernard
C
Mcintosh
AR
Jirsa
VK
.
Selective activation of resting-state networks following focal stimulation in a connectome-based network model of the human brain
.
eNeuro
 
2016
;
3
:
1
17
.
Suh
M
Hongtao
M
Mingrui
Z
Sharif
S
Schwartz
TH
.
Neurovascular coupling and oximetry during epileptic events
.
Mol Neurobiol
 
2006
;
33
:
181
97
.
Talairach
J
Bancaud
J
.
Lesion, ‘irritative’ zone and epileptogenic focus
.
Confin Neurol
 
1966
;
27
:
91
4
.
Taylor
PN
Goodfellow
M
Wang
Y
Baier
G
.
Towards a large-scale model of patient-specific epileptic spike-wave discharges
.
Biol Cybern
 
2013
;
107
:
83
94
.
Terry
JR
Benjamin
O
Richardson
MP
.
Seizure generation: the role of nodes and networks
.
Epilepsia
 
2012
;
53
:
e166
9
.
Thom
M
Sisodiya
S
Harkness
W
Scaravilli
F
.
Microdysgenesis in temporal lobe epilepsy. A quantitative and immunohistochemical study of white matter neurones
.
Brain
 
2001
;
124
:
2299
309
.
Touboul
J
Wendling
F
Chauvel
P
Faugeras
O
.
Neural mass activity, bifurcations, and epilepsy
.
Neural Comput
 
2011
;
23
:
3232
86
.
Tournier
JD
Calamante
F
Connelly
A
.
Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained super-resolved spherical deconvolution
.
Neuroimage
 
2007
;
35
:
1459
72
.
Tournier
JD
.
MRtrix package, Brain Research Institute, Melbourne, Australia [Internet]
.
Toyoda
I
Bower
MR
Leyva
F
Buckmaster
PS
.
Early activation of ventral hippocampus and subiculum during spontaneous seizures in a rat model of temporal lobe epilepsy
.
J Neurosci
 
2013
;
33
:
11100
15
.
Ursino
M
La Cara
GE
.
Travelling waves and EEG patterns during epileptic seizure: analysis with an integrate-and-fire neural network
.
J Theor Biol
 
2006
;
242
:
171
87
.
Valentín
A
Alarcón
G
Honavar
M
García Seoane
JJ
Selway
RP
Polkey
CE
et al
.
Single pulse electrical stimulation for identification of structural abnormalities and prediction of seizure outcome after epilepsy surgery: a prospective study
.
Lancet Neurol
 
2005
;
4
:
718
26
.
Watts
DJ
Strogatz
SH
.
Collective dynamics of ‘small-world’ networks
.
Nature
 
1998
;
393
:
440
2
.
Wendling
F
Bartolomei
F
Bellanger
JJ
Chauvel
P
.
Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition
.
Eur J Neurosci
 
2002
;
15
:
1499
508
.
Zalesky
A
Fornito
A
Harding
IH
Cocchi
L
Yücel
M
Pantelis
C
et al
.
Whole-brain anatomical networks: does the choice of nodes matter?
.
Neuroimage
 
2010
;
50
:
970
83
.
Zhao
M
Nguyen
J
Ma
H
Nishimura
N
Schaffer
CB
Schwartz
TH
.
Preictal and ictal neurovascular and metabolic coupling surrounding a seizure focus
.
J Neurosci
 
2011
;
31
:
13292
300
.

Author notes

See Lytton (doi:10.1093/awx018) for a scientific commentary on this article.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Supplementary data