Individual brain structure and modelling predict seizure propagation.

See Lytton (doi:10.1093/awx018) for a scientific commentary on this article.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 ), motor and sensorimotor tasks (Kelso, 1995), resting state (Allen ; Hansen ) and stimulation paradigms (Spiegler ). 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 ; Liao ; Haneef ; Pittau ; Ridley ). 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 , 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 ; Deco ) of healthy subjects, stroke (Falcon ) and schizophrenic patients (Cabral ). 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 ; Wendling ; Lopes da Silva ; Breakspear ; Kalitzin ; Touboul ; Kramer ; Jirsa ). Only recently, however, has propagation of epileptic seizures started to be studied using BNMs (Terry ; Taylor ; Hutchings ). Partial seizures have been reported to propagate through large-scale networks in humans (Bartolomei ) and animal models (Toyoda ). 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 ; Spencer, 2002). As a part of the standard presurgical evaluation, stereotactic EEG is used to help correctly delineate the epileptogenic zone (Bartolomei ). 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 ; Bernhardt ) or other structural alterations in the connectome of epileptic patients (Bonilha ; Besson ; DeSalvo ). 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 (Supplementary Table 1). 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 ). This processing pipeline makes use of various tools such as FreeSurfer (Fischl, 2012), FSL (Jenkinson ), MRtrix3 (Tournier) and Remesher (Fuhrmann ) 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 ). 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 ). 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 ), and improved with anatomically constrained tractography (ACT; Smith ). Tractography was performed using 5 × 106 streamlines and were corrected using spherical-deconvolution informed filtering of tractograms (SIFT; Smith ) 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 ). 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 , where x[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 S simply counts the accordance of each region found in the simulated propagation zone (ensemble PZS) and the predicted propagation zone (ensemble PZP): The distance score S quantitatively estimates the L1-norm between the normalized probability of the predicted propagation and the strength of the stereotactic EEG signal energy (for clinician prediction, the strength was set to 1):

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: with 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 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 ), metabolism (Zhao ), and extracellular levels of ions (Heinemann ). 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 x and y) and spike and wave events (variables x and y) observed in electrographic seizure recordings. The Epileptor model was also shown to reproduce refractory status epilepticus and depolarization block (El Houssaini ). Using time-scale separation as well as evidence from stereotactic EEG recordings, Proix 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 ), extracellular diffusion of potassium (Durand ), and increase of firing rate in remote regions leading to increase of remote extracellular potassium (Avoli ). We used the same approach here in the general case of N coupled Epileptors, which reads for each Epileptor i: where and . K 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, that we varied in this study. To simplify the interpretation, we define with the critical value of excitability. If , a brain region is epileptogenic and seizures are triggered autonomously. Otherwise, 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 ): for each Epileptor i with τ0 = 2857 and I1, = 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 (Supplementary material): with

Numerical simulations

The simulations were performed with The Virtual Brain (TVB; Sanz Leon ) 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.

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, > 0, in yellow), the nodes in the propagation zone have different excitability values (0.5 > Δx0,> −0.5, shades of red), while all the other nodes are not epileptogenic (stable state, Δx0,< 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.

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,, 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 (Supplementary material). 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 (Supplementary Fig. 1A–C).

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 5 was set to operate on the fast time scale, with opposite sign to keep the coupling function excitatory:

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 5.

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 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 k of the connectivity was summed to a draw of a uniform random distribution, whose values were taken between and , with ɛ a chosen percentage.

Log of the connectivity matrix

We simply redefined 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 Supplementary Table 1. 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 ), 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.

Modelling seizure propagation

A BNM (Sanz-Leon ) 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 ). 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 ). 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).

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, = −2.2, Figs 3A, 4A and 5A, regions in yellow), and we set all the other regions not epileptogenic with the same excitability (Δx0, = −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 ) (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 Supplementary Fig. 2. The individual results are given in Supplementary Table 2 and the abbreviations taxonomy in Supplementary Table 3. 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 Supplementary Fig. 2.

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,< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,< −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,< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,< −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,< 0, in yellow), while all the other nodes are equally far from the epileptogenicity threshold (Δx0,< −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.

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.

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 Supplementary Fig. 3A), but the group level comparison did not show significant differences of the medians (Mann-Whitney U-test, Supplementary Fig. 2, 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 Supplementary Fig. 3B), 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, Supplementary Fig. 2, 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 (Supplementary Fig. 6, 20% and 40% boxplots). As expected the results did not change significantly. The weights show a truncated power-law distribution (Supplementary Fig. 7). Other studies have used a rescaling of the distribution of weights by compressing the range of edge weights (Hagmann , 2008), but our results were degraded by taking the logarithm of the weight matrices (Supplementary Fig. 6, log boxplot).

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

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 Supplementary material). 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 (Supplementary Fig. 8B) as well as for patient connectivity matrices (Supplementary Fig. 6, 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 (Supplementary Fig. 8C–E) or patient-specific (Supplementary Fig. 6, 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 ), 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 ; Schevon ; Martinet ). The slow permittivity variable is supported by a variety of biophysical mechanisms that act on slow time scales such as tissue oxygenation (Suh ), extracellular level of ions (Heinemann ), and metabolism (Zhao ), which were found to vary in mouse (Jirsa ), cat (Moody ) and baboon (Pumain ) models of epileptic seizures.

Most computational models of seizure propagation focus on small continuous spatial scales (Ursino and La Cara, 2006; Kim ; Hall and Kuhlmann, 2013) or population of neurons (Miles ; Golomb and Amitai, 1997; Compte ; Fröhlich ; Bazhenov ). 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 ) or temporal lobe epilepsy (Hutchings ). Other modelling studies focused on small networks to investigate the role of the topology and localization of the epileptogenic zone (Terry ). 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 ; Bonilha ; Bernhardt ; Besson ; DeSalvo ) and cognitive impairments (Leyden ). Histological studies provide evidence of white matter alterations in temporal lobe epilepsy (Thom ; Blanc ). Functional, volumetric and electrographic data suggest a broad reorganization of the networks in epileptic patients (Lieb , 1991; Cassidy and Gale, 1998; Rosenberg ; Bettus ). Hemispherectomy (Beier and Rutka, 2013) or regional disconnection procedures (Mohamed ), 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 ). However, using ACT (Smith ) and SIFT (Smith ) 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 ). 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 ), 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 ), 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 ). Analyses of BNMs have shown that such biomarkers depend on the coupling assessing the threshold level separating ictal from interictal state (Jirsa ; Proix ). 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 ); or alternatively epileptogenic zone estimates may be obtained based on other metrics such as graph theoretical estimates (Burns ). 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 ), 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 (Supplementary Fig. 2); 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 ). More generally, several previous studies have demonstrated a link between anatomical networks and pathologies, in particular neurodegenerative diseases (Seeley ). 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.

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.

Supplementary material

Supplementary material is available at Brain online.

Click here for additional data file.