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.
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 (
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.
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 , 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):
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:
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:
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):
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 (
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 (
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.
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:
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
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.
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 and , with ɛ a chosen percentage.
Log of the connectivity matrix
We simply redefined as the log of the connectivity matrix.
We compared the different group medians for different conditions using first a Kruskal-Wallis test, followed by pairwise Mann Whitney U-tests.
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 inFig. 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.
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).
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 inFig. 6 and
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).
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 andRubinov and Sporns, 2011). The scores were always higher with the individualized connectivity matrices (Fig. 3E and 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 (Hagmann et al., 2007, 2008), but our results were degraded by taking the logarithm of the weight matrices (
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
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 (
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 (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.
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.
brain network model