Skip to main content
Advertisement
  • Loading metrics

An optimal strategy for epilepsy surgery: Disruption of the rich-club?

  • Marinho A. Lopes ,

    Roles Conceptualization, Formal analysis, Methodology, Writing – original draft, Writing – review & editing

    m.lopes@exeter.ac.uk

    Affiliations Living Systems Institute, University of Exeter, Exeter, United Kingdom, Wellcome Trust Centre for Biomedical Modelling and Analysis, University of Exeter, Exeter, United Kingdom, EPSRC Centre for Predictive Modelling in Healthcare, University of Exeter, Exeter, United Kingdom

  • Mark P. Richardson,

    Roles Funding acquisition, Project administration, Writing – review & editing

    Affiliations EPSRC Centre for Predictive Modelling in Healthcare, University of Exeter, Exeter, United Kingdom, Institute of Psychiatry, Psychology and Neuroscience, Kings College London, London, United Kingdom

  • Eugenio Abela,

    Roles Data curation, Investigation, Resources, Writing – review & editing

    Affiliations Institute of Psychiatry, Psychology and Neuroscience, Kings College London, London, United Kingdom, Support Center for Advanced Neuroimaging (SCAN), University of Bern, Bern, Switzerland

  • Christian Rummel,

    Roles Data curation, Investigation, Resources, Writing – review & editing

    Affiliation Support Center for Advanced Neuroimaging (SCAN), University of Bern, Bern, Switzerland

  • Kaspar Schindler,

    Roles Data curation, Funding acquisition, Investigation, Resources, Writing – review & editing

    Affiliation Department of Neurology, Inselspital, Bern, Switzerland

  • Marc Goodfellow ,

    Contributed equally to this work with: Marc Goodfellow, John R. Terry

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

    Affiliations Living Systems Institute, University of Exeter, Exeter, United Kingdom, Wellcome Trust Centre for Biomedical Modelling and Analysis, University of Exeter, Exeter, United Kingdom, EPSRC Centre for Predictive Modelling in Healthcare, University of Exeter, Exeter, United Kingdom

  • John R. Terry

    Contributed equally to this work with: Marc Goodfellow, John R. Terry

    Roles Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review & editing

    Affiliations Living Systems Institute, University of Exeter, Exeter, United Kingdom, Wellcome Trust Centre for Biomedical Modelling and Analysis, University of Exeter, Exeter, United Kingdom, EPSRC Centre for Predictive Modelling in Healthcare, University of Exeter, Exeter, United Kingdom

Abstract

Surgery is a therapeutic option for people with epilepsy whose seizures are not controlled by anti-epilepsy drugs. In pre-surgical planning, an array of data modalities, often including intra-cranial EEG, is used in an attempt to map regions of the brain thought to be crucial for the generation of seizures. These regions are then resected with the hope that the individual is rendered seizure free as a consequence. However, post-operative seizure freedom is currently sub-optimal, suggesting that the pre-surgical assessment may be improved by taking advantage of a mechanistic understanding of seizure generation in large brain networks. Herein we use mathematical models to uncover the relative contribution of regions of the brain to seizure generation and consequently which brain regions should be considered for resection. A critical advantage of this modeling approach is that the effect of different surgical strategies can be predicted and quantitatively compared in advance of surgery. Herein we seek to understand seizure generation in networks with different topologies and study how the removal of different nodes in these networks reduces the occurrence of seizures. Since this a computationally demanding problem, a first step for this aim is to facilitate tractability of this approach for large networks. To do this, we demonstrate that predictions arising from a neural mass model are preserved in a lower dimensional, canonical model that is quicker to simulate. We then use this simpler model to study the emergence of seizures in artificial networks with different topologies, and calculate which nodes should be removed to render the network seizure free. We find that for scale-free and rich-club networks there exist specific nodes that are critical for seizure generation and should therefore be removed, whereas for small-world networks the strategy should instead focus on removing sufficient brain tissue. We demonstrate the validity of our approach by analysing intra-cranial EEG recordings from a database comprising 16 patients who have undergone epilepsy surgery, revealing rich-club structures within the obtained functional networks. We show that the postsurgical outcome for these patients was better when a greater proportion of the rich club was removed, in agreement with our theoretical predictions.

Author summary

Epilepsy is a chronic neurological disorder that affects about 1% of people worldwide. The administration of antiepileptic drugs is the preferable treatment, but in around one third of cases, drugs do not stop seizures, and these patients are potential candidates for surgery. Epilepsy surgery however is too often unsuccessful, with around one half of patients continuing to experience seizures. In this work we use mathematical models to study epilepsy surgery so to inform surgeons concerning the brain tissue that should be considered for surgery resection. We show that functional networks derived from data of epileptic patients considered for surgery present rich-club organization. For this kind of network structure, we propose an optimal surgery strategy that consists of disrupting the rich-club.

Introduction

Epilepsy is a chronic neurological disorder that affects about 1% of people worldwide [1]. Antiepileptic drugs are the preferred treatment, but in around one third of cases, drugs do not stop seizures, and patients for whom this is the case are potential candidates for surgery [2]. Surgeons use an array of data modalities, including intra-cranial electroencephalogram (iEEG), in an attempt to map regions of the brain thought to be crucial for the generation of seizures [3]. If these regions of the brain are amenable to surgery (e.g. they do not overlie eloquent cortex), then they are removed with the hope that the individual is rendered seizure free as a consequence. However, long-term success rates from surgery may be as low as 15%, presumably in part due to failures of the assumptions used in the decision making process [4,5]. It is therefore crucial to advance our understanding of the mechanisms that generate seizures and the reasons why removing regions of brain tissue may or may not lead to seizure freedom.

In this regard, seizures are increasingly recognised as arising in large-scale brain networks [69]. Emerging from such networks, both healthy and pathological dynamics are observed, for example through EEG, MEG or fMRI. These dynamics emerge due to the interplay between intrinsic properties of brain areas, structural connectivity, and modulating influences across multiple temporal and spatial scales [1012]. This networks paradigm has led to imaging or electrographic data being used to inform network representations of the brain (for example structural or functional brain networks), and graph theoretical measures are used to characterise the topology of these networks [1319]. Studies analysing graph theoretical properties of networks have reported differences between functional and structural networks derived from healthy individuals versus people with epilepsy [2026].

The emerging field of dynamics on networks is complementary to these traditional, “static” network analyses [27,28], and moves beyond the study of the topology of networks. In this approach, mathematical models are used to link networks and the intrinsic properties of individual nodes to dynamic data [29], which provides an avenue to understand the relationship between structure and function [30]. In particular, mathematical models that can recreate elements of pathological dynamics, for example the occurrence of seizures, have been used to understand the network mechanisms of disorders such as epilepsy [9,3137]. Such approaches are also being used in translational applications, for example providing additional information to complement clinical interpretation, namely within the diagnosis of epilepsy [35,38]. Crucially, a dynamics on networks approach can be extended to study perturbations to networks. On one hand, lesions and traumatic brain injury can lead to the emergence of pathological brain activity, on the other hand, perturbations such as pharmacological treatment, single pulse electrical stimulation (and other electrical stimulations), transcranial magnetic stimulation, thermocoagulation, among others, can transform brain dynamics from pathological to healthy states [36,3941], therefore revealing potential avenues for therapy. In the case of epilepsy surgery, we have demonstrated that a network model derived from iEEG data could provide relevant predictions for the outcome of epilepsy surgery [42]. Our findings have been recently replicated in an independent cohort of 16 people with pharmacoresistant epilepsy [43] offering further support to a dynamics on network approach. However, the ways in which networks with different topologies respond to perturbations is at present unknown. For example, in analogy to epilepsy surgery, it is unclear whether particular networks are amenable to a reduction in pathological dynamics upon removing nodes and if so which nodes would be best to target.

Here, we use a dynamics on networks approach to study the generation of pathological activity in networks and how the removal of nodes can restore healthy dynamics. Our starting point is a neural mass model that has previously been shown to generate epileptiform rhythms in focal seizures [32,37,44], and that we have successfully used to quantify and predict the outcome of epilepsy surgery [42]. It has been shown that the model, when placed close to a saddle-node on invariant circle (SNIC) bifurcation, can generate spontaneous, recurrent transitions to epileptiform dynamics (both inter-ictal spikes as well as seizures) when driven by noise [32,37,42]. In our framework, the neural mass model describes the dynamics of a single node within a wider network. The systematic exploration of node removal in brain networks is computationally demanding, and hence we seek a computationally efficient version of this model that preserves the quantification of the effect of removing nodes. We show that a modification of the theta-neuron model [45] is appropriate for this purpose since it is the canonical form of the bifurcation under consideration. This model is capable of generating spiking dynamics, which here represents seizure-like activity.

The computational benefits of the theta-neuron model allow us to study the emergence of spiking dynamics in different types of networks and also to systematically quantify the effect of removing different nodes. Here, we study small-world, random, rich-club and scale-free and find that rich-club and scale-free networks more readily generate spiking dynamics, since they require a lower strength of coupling between connected nodes to do so. In terms of the contribution of nodes, we find that rich-club, random and scale-free networks possess a small number of nodes that drive spiking dynamics, whereas the propensity of generate spiking dynamics is more evenly distributed across nodes in small-world networks. Collectively, this suggests that patients whose brain networks display rich-club properties should be particularly amenable to current surgery paradigms. In order to test the relevance of these findings, we analyse data from patients who underwent surgery and for whom postoperative outcome is known. We demonstrate that functional networks inferred from iEEG during seizures display a rich-club connectivity structure and that the proportion of rich-club nodes removed correlates with the success of surgery.

Methods

Ethics statement

This study was approved by the Internal Review Board of the Inselspital (approval No. 159399, dated 26th of November, 2013). All patients gave written informed consent that imaging and EEG data may be used for research purposes.

Wendling model

In order to model epilepsy surgery, we consider large-scale brain networks, where each network node is capable of generating epileptiform activity but will do so depending on the connectivity structure of the network. In this framework, a node putatively represents a portion of brain tissue potentially responsible for the emergence of seizure activity across the network. We assume that the dynamics of each node can be described by a neural mass model, such as the Wendling model [37,42]. The model depicts the dynamics of a macroscopic circuit in which a population of excitatory pyramidal neurons interacts with three populations of interneurons (representing one excitatory and two inhibitory populations). The two inhibitory populations are classed as slow and fast, representing dendritic-projecting GABAA and somatic-projecting GABAA interneurons, respectively. The dynamics is described by the following 10 first-order ordinary differential equations (ODEs): where z1-z5 are the output potentials in mV of the neuronal populations, namely z1, z2, z3, and z4 are the outputs of the pyramidal cells, excitatory population, slow inhibitory population, and fast inhibitory population, respectively. z5 is the output of the slow inhibitory population that interacts with the fast inhibitory population. z6-z10 are auxiliary variables, S is a sigmoid function, and A, a, B, b, G, g, C1-C7, p, e0, r, and ν0 are parameters (see Table 1 for their biophysical interpretation and values).

thumbnail
Table 1. Model parameter values and biophysical interpretations.

https://doi.org/10.1371/journal.pcbi.1005637.t001

The output of the model z2(t)–z3(t)–z4(t) corresponds to the aggregated membrane potential of the excitatory cell population and its bifurcations have been extensively characterized [47]. In particular, a SNIC bifurcation has been identified as one mechanism for the generation of epileptiform rhythms observed in typical focal epilepsies [32]. This model and bifurcation were also previously employed to estimate brain network ictogenicity to predict the outcome of epilepsy surgery [42]. Therefore the parameters A, B and C were chosen so that the neural mass is in a steady state close to the SNIC bifurcation that gives rise to spiking dynamics which we consider a proxy for the patho-phenotype of the epileptic brain (see the third figure, left panel, in [32]). p is an extrinsic input parameter that represents stimuli from other areas of the cortex.

Although the neural mass model described above represents the dynamics of four interacting neuronal populations, at the scale we are interested in, it describes the dynamics of a single node in a wider network consisting of other interacting neural masses. Following previous studies [33,48], we account for the coupling between neural masses (nodes) using the extrinsic input parameter p. We make the input of the j-th node both time and node dependent as follows,

Here the index j denotes node j (j = 1,2,…,N, where N is the number of nodes). is used to control the distance to the SNIC bifurcation; ξ(j)(t) represents noisy inputs from other areas of the cortex outside of the network under consideration; λij is the coupling strength from node i to node j; and aij is the i,jth entry of the adjacency matrix (the node receives the outputs of all his in-neighbours) [33,48]. We consider Gaussian noise with mean and where is the variance. A node is in a resting state if pj(t) < pc, where pc is the critical point at which the SNIC bifurcation takes place.

Canonical model

Since the Wendling Model (WM) becomes computationally expensive for studying large networks, we look for a parsimonious representation for spiking dynamics in brain networks. Taking into account that nodes of WM are operating in the vicinity of a SNIC bifurcation, we substitute networks of neural masses with networks in which each node is represented by the normal form of the SNIC, i.e. the theta-neuron model [45]. It is important to stress that although this model is traditionally used to describe the dynamics of a neuron, here we use it (as effectively the canonical form of the SNIC bifurcation) to represent the dynamics of a neural mass in an epileptic spiking regime. The canonical model (CM) is an alternative formulation of a quadratic integrate and fire neuron. It comprises the following ODE: where θj is the phase of node j, and Ij(t) is its input current. The SNIC bifurcation occurs at Ic = 0. At Ij < Ic, the phase oscillator is resting, whereas at Ij > Ic it is oscillating.

We define the coupling between the “canonical neural masses” analogous to the coupling defined within the WM, where Ij is the input current of node j, represents noisy inputs coming from other areas, wij is the coupling strength from node i to node j, and aij is the i,jth entry of the adjacency matrix. As in the WM, we consider Gaussian noise (mean , and variance ). We define the output of the in-neighbour i as , where is its steady state, so that if the node is resting its output is zero, and if it reaches , its output is maximum. This uncoupled steady state is obtained from setting ,

We take the real part so that at . At , there are two fixed points: is a stable fixed point, and is an unstable fixed point. A similar coupling in networks of theta-neurons was recently studied in [49]. Other authors have considered delta-like interactions [50], or rapid rises in the synaptic gating variable [51], which are a reasonable approximation for neurons, but inappropriate for neural masses. Note that the output of a neural mass is an average over the activity of a population of neurons, and so it displays properties of a low-pass filter [52].

Parameter comparison

For simplicity, we consider homogeneous nodes in both models, i.e., all nodes in a network are at the same distance to the SNIC bifurcation ( and ), and have the same coupling strength (λij = λ and wij = w). This is a strong assumption that enables us to focus explicitly on the contribution of the network structure to the network ictogenicity. Thus, there are three free parameters in each model: (p0,σp,λ) in WM, and (I0,σI,w) in the CM. Since our aim is to consider whether the two network models display similar changes in dynamics upon the removal of nodes, it is important that these parameters are comparable between models. Taking into account that we require that the node dynamics switch between the resting state and the spiking dynamics, the three parameters are interdependent. For example, as parameter values of the nodes move closer to the SNIC the required noise variance to elicit spikes becomes smaller. Note, however, that the variance of the noise should not be too large as we wish to ensure that network interactions play a role in the emergent dynamics. Thus, we define and to scale the effect of noise by the distance to the SNIC bifurcation so that the effect of the noise on the dynamics of both models is comparable. In order to establish a relation between the coupling strength and the noise, we also define λ* = 2e0/(p) and w* = 2cw/(I), where c is the mean degree of the network. These relations compare the noise to the average maximum input that a node can receive, 2e0/N and 2cw/N for WM and the CM, respectively. It provides a scale that compares noise perturbations to inputs received from in-neighbours.

Note that, with respect to the input parameter, the dynamics of a node j change from resting to spiking in WM if pj(t) > pc, and likewise, in the CM a node j transitions to spiking if Ij(t) > Ic. Thus, in both models we have the following condition for a node j to be in the parameter region corresponding to a spiking regime at time t, where is the noise, C the homogeneous coupling strength, Yi(t) the output of node i, and T the bifurcation point. If we assume that the network is in the resting state with an average node output of 〈Y〉, then we can estimate the critical coupling Cc at which on average a certain node starts to spike, where is the in-degree of node j (〈ξ(j)(t)〉 = 0 in the case of Gaussian noise). Therefore, for a given network of size N, the larger the in-degree, the smaller is Cc, meaning that nodes with higher in-degree are more likely to transition to spiking. This is valid in both models. Similarly, one can find the critical distance to the SNIC bifurcation, which is smaller than T due to the inputs from the network (〈Y〉 > 0).

Artificial networks and measurements of nodes

The adjacency matrix encodes the network structure on top of which the nodes interact. We consider random, scale-free, small-world and rich-club networks, both directed and undirected (we discarded networks with disconnected components) [53,54]. In order to quantify the “importance” of each node, we analyze the following traditional measures: degree, average neighbour degree, eigenvector centrality, betweenness centrality, closeness centrality, clustering coefficient, and local efficiency [55,56]. Additionally, we also consider eigencentrality based on Jaccard dissimilarity [57] and dynamical importance [58]. In the case of directed networks, we also consider in-degree, out-degree, as well as the sum and product of these measures.

Definition of spiking dynamics and treatment perturbations

We focus our analysis upon two measurements that are relevant for our purposes of studying epileptic dynamics and surgery in silico, namely Brain Network Ictogenicity (BNI) [23,42,59], and Node Ictogenicity (NI) [42]. BNI is a practical approach for quantifying the tendency of a network to generate spiking dynamics. It measures the average fraction of time spent in spiking dynamics by each node [23,42,59]:

Specifically, in the WM, first we extract the spikes generated by a node by applying a threshold to the average absolute amplitude of the model output over a sliding window of 0.05 s. Then, contiguous epochs of spiking dynamics are identified by evaluating the overlap of 1 s time windows centred in each spike. Finally, the time spent in spiking dynamics corresponds to the total time of these spiking epochs [42]. In the CM we use the same method, with similar time scales (we use as conversion time scale the ratio of the full widths at half maximum of the spikes in each model).

NI quantifies the contribution of each node to the ictogenicity of the network by measuring the relative difference in BNI upon removing node i from the network: where BNIpre corresponds to the BNI over the network prior to node resection and is the BNI after the removal of node i. Note that NIi = 1 means that the removal of node i renders the network free of spiking dynamics, whereas NIi = 0 means that the resection of node i made no difference to the BNI. In practice, this quantity measures the success of a given surgery resection in silico, and it may have the potential to guide the search for an optimal surgical strategy. In general, this quantity may also be useful to quantify the result of temporary ablation, assuming that the ablation takes place in a much slower time scale than the network dynamics. In this paper we set BNIpre = 0.5 (we have confirmed that the results are qualitatively the same for other reference values of BNIpre).

To evaluate if the CM can be used as a proxy of the WM in this framework, we compare the NI ordering of the two models for a number of networks. Note that NI is essentially a vector with N entries quantifying the result of removing each node individually, being of particular interest the relative impact of each node removal compared to the others, rather than the absolute value of each one (which is parameter dependent). We use a weighted Kendall's rank correlation measure [60,61], which is defined as follows. Given two rankings (NI) of the same items (nodes of the network), we calculate where P is the number of items in the same order in the two rankings, and Q counts the number of items in reverse order. When τ = 1 the two rankings predict the same ordering, whereas τ = −1 means a reverse order of all items. Here we consider a weighted measure to take into account the relative values of NI: each NI comparison between two nodes i and j is weighted by the product of the distances in NI predicted by the two models, , (where and are the NIs of node i calculated using WM and the CM, respectively). We assume that there are no ties.

Patient data and functional networks

We focus on patients with pharmacoresistant epilepsy, since such patients are candidates for surgery. Data were collected from 16 patients (11 female, mean age 31, and median post-surgical follow up 3 years) who underwent pre-surgical monitoring at Inselspital Bern [42,62]. Following epilepsy surgery, six patients fell into Engel class I (free of disabling seizures), five into Engel class II (rare disabling seizures) and five into Engel class IV (no worthwhile improvement). All patients gave written informed consent that imaging and iEEG data may be used for research purposes. Other details about the data can be found elsewhere [42,62]. Before analysis, the signals were down-sampled to a sampling rate of 512 Hz and re-referenced against the median of all the channels free of permanent artefacts as judged by visual inspection by an experienced epileptologist (K.S.). For each patient, two peri-ictal epochs were considered, which included three minutes before seizure onset, the seizure itself and three minutes after seizure termination (seizure onset and offset were identified by visual inspection (K.S.)). Following the methods described in [62,63], first we applied a band-pass filter between 0.5 and 120 Hz and a notch filter (48 to 52 Hz) using a Butterworth filter. Each epoch was divided in a set of 8 seconds segments (the segments were chosen 1 second apart from each other). For each segment we obtained 10 univariate iterated amplitude adjusted Fourier transform (IAAFT) surrogates independently. Next, the segments were divided in 10 subsegments of 1024 sampling points (2 seconds) distributed with minimal overlap. Thus, we generated an ensemble of 10 subsegments for the original time series, and 100 subsegments for the surrogates (10 for each surrogate). To estimate the correlations between the time series of each iEEG channel, we used the Pearson's equal-time (zero-lag) cross-correlation coefficient ρ, and a non-parametric Mann-Whitney-Wilcoxon U-test was performed to assess the significance of different medians of ρ between the original time series (ρo) and the surrogates (ρsurr). We further applied Bonferroni-Holm corrections to account for multiple comparisons. Finally, we obtained a surrogate-corrected correlation matrix using the heuristic formula [63,64] where s = 1 if the null hypothesis of the statistical test is rejected, or s = 0 otherwise. Using this method, we derived 102 ± 18 functional networks based on cross-correlation for each patient, depending on the duration of each seizure epoch.

Rich-club organization

The organization of functionally derived networks into rich-clubs [6567] was studied using a weighted rich-club parameter ϕ(k) [66]. The richness parameter is the degree k, and the procedure consists in finding groups (clubs) of nodes whose richness is larger than k. For a given degree k, we counted the number of connections E>k of the club, and summed their weights W>k. We then calculated the fraction of weights shared by the club out of the maximum edge weights that the club could have if they were linked by the strongest connections of the network, i.e., where are the ranked weights of the network. This fraction is not enough to verify the existence of a rich-club, since even random networks can have an increasing function ϕw(k) as a result of chance alone (nodes with higher degree are more likely to be connected). Therefore, ϕw(k) is normalized relative to ϕrand(k) obtained from a set of comparable random networks,

Thus, a network exhibits rich-club organization if there is a range of degree k for which ϕ(k) > 1 [65,66]. We generated 100 random networks by applying a reshuffle procedure to the weights while keeping the topology of the original network intact, followed by a link and weight reshuffle procedure that preserves the original degree distribution [56,67]. ϕrand(k) was calculated as the average rich-club coefficient for each level of k. Finally, we evaluate the statistical significance of rich-club organization using a permutation test [67], by testing whether ϕ(k) was statistically significantly larger than ϕrand(k) (a one-sided p value was calculated as the percentage of the distribution of ϕrand(k) that exceeds ϕ(k)). We measured the rich-clubs of the average functional networks of the pre-seizure, seizure, post-seizure, and whole peri-ictal epochs for each patient separately.

Results

Predictions for the effect of node removals are preserved in a canonical model

We compared the dynamics of the CM to the WM in terms of the effect that model parameters have on BNI and the profile of NI for a suite of networks. Fig 1 demonstrates typical dynamics of each model applied to the same network (a directed random network with N = 10, and mean degree c = 1.6). Both models display spiking dynamics, with a heterogeneous distribution of activity across nodes. For each model, nodes 2, 5, 7, 9 and 10 show a greater extent of spiking dynamics than other nodes; thus the distribution of activity across the network is preserved in the canonical model. On the other hand, it is clear that the resting state is noisier in the CM. A predominant feature accounting for this is that the ratio of amplitude of the spiking trajectory to noise is larger in the WM. Moreover, one should also realize that whereas in the WM only positive inputs can move the system towards the SNIC bifurcation, in the CM both positive and negative inputs displace the phase, θ, from the resting state. Furthermore, the output of the resting state in the CM is zero, but non-zero in the WM.

thumbnail
Fig 1.

Node activities in the (a) WM and (b) CM. Although the spikes have different forms between models and the CM shows higher noisy activity, both models identify the same set of nodes as having a larger propensity to spike. The network is the same for both models, a directed random network with N = 10 nodes and mean degree c = 1.6. Parameters: and .

https://doi.org/10.1371/journal.pcbi.1005637.g001

As described in the Methods, although we have identified three free parameters, the network dynamics are in fact affected by only two competing factors: the distance to the SNIC bifurcation and the coupling strength. The strength of noise required to elicit spikes is correlated with the distance to the SNIC bifurcation (that is smaller noise variance is required to elicit spikes the closer the system is to the bifurcation). Thus, we can fix the noise variance and consider BNI as a function of the distance to the SNIC bifurcation and coupling strength. Fig 2 provides an evaluation of this function for an ensemble of 10 random networks with 15 nodes and demonstrates that the smaller the distance to the bifurcation, the easier it is to generate spiking dynamics, and consequently BNI is larger. In addition, BNI grows with increases in coupling strength. Fig 2 demonstrates that the shape of the BNI surface is similar for the two models, which provides evidence that the normal form of the SNIC is appropriate for the study of the propensity of a network to generate spiking dynamics. Similar results were obtained for both smaller and larger networks.

thumbnail
Fig 2. BNI as function of coupling strength and distance to the SNIC bifurcation.

The left and right columns correspond to the WM and CM, respectively. The resemblance between the two surfaces suggests that the CM is a good parsimonious model of WM for the of study BNI. The red lines on the surfaces correspond to BNIpre = 0.5. The surfaces in panels (a) and (b) represent an average over an ensemble of 10 directed random networks, and the error bars in panels (c)-(f) account for the variability of BNI between different network realizations. Parameters: , and c = 1.6. Additionally, in panels (c) λ* = 1.09, (d) w* = 1.33, (e) pcp0 = 1.45, and (f) IcI0 = 1.62.

https://doi.org/10.1371/journal.pcbi.1005637.g002

Our results thus far indicate that despite some expected quantitative differences, network dynamics, and in particular the way that BNI changes with respect to system parameters, are qualitatively similar across the WM and the CM. However, our primary focus is to determine whether the CM would provide the same prediction for the effect on BNI of node removals (i.e. NI). With the application of surgical resections in mind, we are predominantly interested in how comparable the ordering of NI is between the two models. In order to investigate this, we calculated the distribution of NI for a suite of random networks of size N = 15,30 and 50 and calculated the similarity in ordering of NI using Kendall's τ (see Methods). Fig 3 shows that within models, the NI distribution is robust across different parameter sets for which BNIpre = 0.5, which is our starting point for the calculation of NI (see Methods) and defines a line in the surface of Fig 2. Across different choices of parameters within the WM, we find τ > 0.97 for all networks considered, indicating a strong preservation of the ordering of NI when different parameters are used. Within the CM, τ > 0.89 and thus there is slightly more variation across NI orderings for this model. Fig 3C and 3D show τ for comparisons of NI orderings between the two models. Fig 3C demonstrates that when model parameters are chosen randomly, the ordering of NI is preserved between the two models for small networks, but differences in predictions between the two models arise in larger networks (for example with 50 nodes). However, Fig 4D demonstrates that a parameter set for each model can be found such that NI distributions are preserved across models in the larger networks studied (50 nodes, τ = 0.85 ± 0.09).

thumbnail
Fig 3. Consistency of NI within and between models and computational gain.

The consistency of ordering nodes according to their NI value, quantified by Kendall correlation rank τ within (a) the WM, (b) the CM, and (c),(d) comparison between the two models. The comparison is performed for networks of size N = 15, 30 and 50. In each scenario 10 realizations of directed random networks with mean degree c = 1.6, 3.2, and 5.3 are considered. For each case, the NI ranking is contrasted across 5 different parameter realizations. In panel (c) all combinations of parameters are compared, whereas in panel (d) the parameter set that maximizes τ is used. Although the distribution of NI tends to differ between the models as N increases, panel (d) shows that it is possible to find parameters such that the NI ordering is consistent between models. The red bars in panels (c) and (d) correspond to undirected scale-free networks (10 network realizations, N = 50, γ = 2, and c = 6). Panel (e) demonstrates the computational gain of using the CM relative to the WM. The gain corresponds to the computational time the WM takes to compute the BNI divided by the time the CM takes to perform the same task.

https://doi.org/10.1371/journal.pcbi.1005637.g003

We note that as N increases, nodes become topologically similar in a random network, and therefore one can expect a homogeneous distribution of NI. However, we are primarily interested in networks for which nodes exist that should be resected to reduce the presence of spiking dynamics. We therefore study networks for which we might expect the distribution of NI to be heterogeneous (as we will show below). A natural choice is a scale-free network characterized by a power law degree distribution P(k)∼kγ with a small exponent (γ < 3) [54]. Fig 3 demonstrates that for scale-free networks arbitrary choices of parameters yield a strong similarity in ordering of NI (τ = 0.87 ± 0.16) and that model parameters exist for which the ordering is essentially identical (τ = 0.996 ± 0.003).

Fig 3E demonstrates the computational advantage gained by using the CM over the WM. We find that the ratio of computational time of the WM to the CM when estimating BNI is 4.6 for networks of size N = 15, 4.9 for N = 30, and 6.2 for N = 50. Note that this gain does not correspond to the ratio of floating point operations needed by each model to simulate a time step because the time scales are different between these two models. Crucially, such gain will be very useful when applying this framework in the clinical setting, as it represents a speed-up in the computational time from days to hours.

Having demonstrated similarity in the ordering of NI across the WM and CM, we proceed in the following sections to use the CM to study how NI varies across different types of network. We fix the number of nodes that we consider to be 64, in line with a typical number of iEEG and depth electrodes used in pre-surgical planning applications [42,62,68].

The role of topology in network ictogenicity

Fig 4 shows how BNI varies as a function of the coupling strength under different choices of network topology. Scale-free networks are the most prone to transit to spiking dynamics since BNI becomes non-zero for smaller coupling strengths relative to the other topologies. The effect is more noticeable in networks with smaller exponent γ, which have a greater degree variance (i.e. are more heterogeneous). However, the maximal value of BNI is less than one for scale-free networks, in particular in directed networks. Fig 4 demonstrates that rich-club networks exhibit a similar profile of increases in BNI with increases coupling strength features to scale-free networks, which is presumably a consequence of similarities in the degree distributions of these networks. Small-world and random networks have similar profiles, implying that the high clustering coefficient of small-world networks has little impact on a network’s ictogenicity (BNI). Similar results were obtained for smaller and larger networks (up to N = 128), as well as for sparser (c = 4) and denser networks (c = 10), and in the case of small-world networks for smaller rewiring probabilities (p = 0.1). Fig 4B shows that whilst the profile of random and small-world networks is similar for directed and undirected networks, the profile changes for rich-club and scale-free networks, most significantly for scale-free networks, whose BNI in directed networks has a very gradual increase with increasing coupling strength. This is likely due to the disparity between the in- and out-degree of nodes (note that nodes can have a high out-degree but a low in-degree, meaning they can influence the network activity, but not be influenced by it, and vice versa).

thumbnail
Fig 4.

BNI(w*) for different network structures in the CM, namely in panel (a) undirected and in panel (b) directed networks. The red curves correspond to random networks, the blue and green to scale-free (SF) networks with γ = 3 or γ = 2, respectively, the orange to small-world (SW) networks, and the brown to rich-club (RC) networks. Note that scale-free and rich-club networks show a non-zero BNI value at smaller coupling strengths, and take higher strengths to get to BNI = 1 in comparison to random and small-world networks. Each curve is an average over 10 network realizations of the same topology. Parameters: N = 64, c = 6, probability of rewiring for SW networks p = 0.5, and the rich-club has 10 nodes (with p1 = 0.7, and p2 = 0.2 the connection probabilities within and from the rich-club to the rest of the network, respectively).

https://doi.org/10.1371/journal.pcbi.1005637.g004

The role of topology in node ictogenicity

Fig 5 demonstrates the way that NI is distributed amongst nodes in networks with different topologies, and furthermore how NI correlates with graph theoretical properties of nodes. The first column of Fig 5 shows that random, scale-free and rich-club networks each have a skewed NI distribution, with a small subset of nodes having large NI, whereas small-world networks have a flat distribution, with small values of NI across nodes. Nodes in rich-clubs were found to have high NI. Similarly, we found several measures of node importance to correlate with NI, but in a topology-specific way. For example, whilst node degree correlates to a great extent with NI in scale-free networks, it does not for random networks, particularly when the mean degree is large. Eigenvector centrality and dynamical importance were found to be good predictors of NI (R2 > 0.80, see Fig 6) in all networks except those with small-world topology and random networks with large mean degree (c = 10). We found that small-world networks required all considered measures to achieve an adequate prediction of NI (R2 > 0.85 for multiple regressions with all the considered node properties, for c = 4 and c = 6). We note that directed networks typically did not contain nodes with NI > 0.5 and furthermore that NI did not correlate with graph theoretical measures in directed networks (R2 < 0.3, see S1 Fig). We tested the robustness of these results for other reference values of BNIpre (BNIpre = 0.3 and BNIpre = 0.7) in all the networks (c = 6). We obtain similar results, although the networks become less sensible to perturbations for BNIpre = 0.7.

thumbnail
Fig 5. Comparison of NI across different undirected network structures and its relation to node properties in the CM.

Each row corresponds to a different network topology: (a)-(d) random network; (e)-(h) scale-free network (SF) (γ = 3); (i)-(l) scale-free network (γ = 2); (m)-(p) small-world network (SW); (q)-(t) rich-club (RC) network. In the first column the nodes are sorted by their NI, so that it is monotonically increasing; in the second column NI is sorted by the degree of the nodes; in the third column it is sorted by the betweenness centrality (BC); and in the fourth column it is sorted by the eigenvector centrality (EC). A linear correlation is particularly clear between NI and eigenvector centrality. The shaded areas in the first column and the dots in the other panels correspond to 10 different network realizations of the same topology (the line represent the mean). Parameter choices are as in Fig 4.

https://doi.org/10.1371/journal.pcbi.1005637.g005

thumbnail
Fig 6. Coefficient of determination (R2) of linear regressions of NI for a set of node measurements applied in different network topologies for the CM.

Each row corresponds to a node measurement (degree, average neighbour degree, eigenvector centrality, eigencentrality based on Jaccard dissimilarity, dynamical importance, betweeness centrality, closeness centrality, clustering coefficient, and local efficiency), whereas each column corresponds to a different network topology (within topologies, we present results for three different mean degrees). The last row is the R2 of a multiple linear regression combining all measurements. Eigenvector centrality, eigencentrality based on Jaccard dissimilarity, and dynamical importance are the node measures whose correlation with NI is higher. Colour code: green corresponds to R2 > 0.95; blue to 0.90 < R2 ≤ 0.95; orange to 0.80 < R2 ≤ 0.90; and red to R2 ≤ 0.80. All networks are undirected and parameter choices are as in Fig 4.

https://doi.org/10.1371/journal.pcbi.1005637.g006

Targeted and random node removal

Depending upon the particular choice of network representation, resections from brain networks could include more than one node. We therefore sought to gain insight into how many nodes have to be removed from a network in order to render it incapable of generating spiking dynamics. In order to resect a minimum number of nodes while reducing BNI as much as possible, it seems sensible to target the highest NI nodes first. Fig 7 demonstrates how this strategy compares to random node removal (we use the eigenvector centrality as a proxy of NI and therefore we target nodes with highest eigenvector centrality, but similar results were obtained targeting nodes with highest NI). Note that in the case of targeted node removal, when one node is removed, the whole network changes, and therefore the distribution of NI may change as well. Therefore, we recalculated the eigenvector centrality of each node of the new network after each node removal. The figure shows that a targeted node removal is much more effective than a random strategy in all topologies except small-world networks, in which the two strategies give similar outcomes. This is to be expected taking into account the highly homogenous distribution of NI in small-world networks (see Fig 5). Accordingly, the difference between the two strategies is particularly noticeable in scale-free networks, because of its heterogeneous NI distribution.

thumbnail
Fig 7. BNI as function of the number of nodes removed from a network in the CM.

Each panel corresponds to a different undirected network topology: (a) random networks, (b) scale-free networks (γ = 3), (c) small-world networks, and (d) rich-club networks. The two curves represent two strategies of node removal: the blue curves correspond to targeted node removal according to highest eigenvector centrality, whereas the red curves correspond to random node removal. Target node removal is much more effective to reduce BNI in all network topologies except for small-world networks in which both strategies provide a similar outcome. The probability of rewiring for SW networks was p = 0.1, and other parameters choices were as in Fig 4.

https://doi.org/10.1371/journal.pcbi.1005637.g007

Rich-club organization in iEEG-derived functional networks

Our results thus far demonstrate the distribution of NI throughout a network is dependent upon network structure. In particular, rich-club, or networks with highly connected hubs were found to contain nodes with high NI, even though all nodes were equivalently parameterized and therefore those nodes were not apparently pathological. In a practical setting those nodes would be natural targets for epilepsy surgery. We therefore sought to understand whether typically used clinical data would yield network representations of the brain with these properties. We thus quantified the presence of rich-club organization in functional networks derived from iEEG recordings from patients that were considered for epilepsy surgery. We considered peri-ictal recordings, and we found evidence of rich-club structure in the functional networks of pre-seizure, seizure and post-seizure epochs in each patient. Fig 8 shows rich-club functions for 3 representative functional networks from seizure epochs of 3 different patients. A rich-club coefficient larger than one over a range of degree k indicates the presence of rich-club organization [65,66] and we found this to be the case in all patients. Note that scale-free networks also display rich-club organization [65], and so both types of network are identified by this type of analysis.

thumbnail
Fig 8. Rich-club functions ϕ(k) of functional networks of pharmacoresistant epileptic patients.

Each curve ϕ(k) was obtained from an average functional network of a different patient. ϕ(k) > 1 indicates the presence of rich-club organization.

https://doi.org/10.1371/journal.pcbi.1005637.g008

Disruption of rich-club predicts postsurgical outcome

Next we extended this analysis taking into account the location of resections relative to the placement of iEEG electrodes and the postoperative outcome [42,62]. Our model analysis led us to hypothesise that if the rich-club was partially or totally resected, the outcome for patients would likely be favourable, since we would expect nodes in the rich-club to have high NI. To test this, we estimated which nodes were members of the rich-club in each functional network (over the pre-seizure, seizure, and post-seizure epochs and all combined) as the collection of nodes belonging to the ‘richest’ club, i.e., the nodes with degree larger than kr, where kr corresponds to the maximum of ϕ(k) (see Methods). Fig 9 demonstrates that for networks derived from pre-seizure, post-seizure, or full peri-ictal epochs, the fraction of rich-club resected did not correlate with the outcome of surgery. However, in functional networks derived from the seizure epoch, there was a significant difference in fraction of rich-club resected between patients with good (Engel classes I and II) and poor (Engel class IV) post-operative outcome (p = 0.038, Kruskal-Wallis test). Patients with good postoperative outcome (seizure free or almost seizure free) had a significantly larger disruption to the rich-club than those with no postoperative improvement.

thumbnail
Fig 9. Fraction of rich-club (RC) resected versus seizure outcome for 16 patients.

Panels (a), (b), (c), and (d) correspond to a quantification based on functional networks derived from pre-seizure, seizure, post-seizure, and whole epochs, respectively. Braces indicate the p-value of Kruskal-Wallis one-way analysis of variance to quantify the statistical difference between good and bad responders.

https://doi.org/10.1371/journal.pcbi.1005637.g009

Discussion

In this study we used a canonical form of the Wendling model to systematically explore the influence of network topology on the generation of spiking dynamics and the effect that removing nodes from a network has on its ability to generate such dynamics (i.e. its ictogenicity). We demonstrated that networks with scale-free and rich-club topology are more ictogenic in the sense that they require a smaller coupling strength between connected nodes to lead to the onset of spiking dynamics. Furthermore, we showed that on the whole, the ictogenicity of nodes within an undirected network correlate with graph theoretical measures most notably degree, eigenvector centrality and dynamical importance [58]. This led us to hypothesise that disruption of rich-clubs in networks should lead to diminished ictogenicity. We tested this hypothesis by first demonstrating the presence of rich-clubs in functional networks derived from iEEG data of people with pharmacoresistant epilepsy, and further showing a significantly greater extent of disruption to rich-club structures in patients who had good postoperative outcome, compared to those with poor postoperative outcome.

It has recently been suggested that there exists a local pathological hub near the epileptic focus responsible for spreading epileptic activity, which should be resected by surgery [69]. This hypothesis is supported by evidence demonstrating that betweenness centrality correlates with resected cortical regions in patients who had a favourable surgery outcome [70,71]. Our results are in agreement with this hypothesis. Indeed, the rich-club comprises a group of nodes with high betweenness centrality that can work as pathological hub spreading epileptic activity.

Functional networks can be thought of as representing communication pathways in the brain that are active, or open, at a given moment. In the case of brain disorders such as epilepsy, we assume that there exist pathological pathways that support the emergence of seizures (i.e. ictogenic networks). Since it is natural to assume that these networks are also expressed in electrographic data and can be characterized in terms of functional connectivity, we use such networks in our models. This is in contrast to modeling based on the structural connectivity of the brain [73]. Structural networks place a constraint on the rhythms of the brain that can emerge. However, network models of brain rhythms built on structural networks typically neglect biochemical processes taking place over multiple different time scales preventing us from knowing, at any given time, which connections are actually active. That is why we have chosen functional connectivity based modeling in this study. A potential limitation of this approach, however, is that the network studied depends on the choice of epoch and the method used to characterize functional connectivity. Interestingly, in our analysis a significant difference in fraction of the rich-club resected was only found for functional networks derived from epochs containing seizures. This is in line with our previous analysis in which the effect of resections could be predicted based on functional networks derived from seizure epochs [42]. Furthermore, our previous study utilised functional networks derived from a nonlinear channel association measure (the surrogate corrected mutual information [62,63]), and we have therefore demonstrated that information relevant to the ictogenic network can be extracted from spatiotemporal seizure dynamics using both linear and nonlinear measures to infer the connectivity structure between nodes. Previous analyses of scalp EEG or MEG have demonstrated that information relevant to epilepsy is also present in background epochs (i.e. those not containing seizures) [22,23,35,38,59,72]. In particular, the predictive power of models for resective surgery has also been studied using inter-ictal (e.g. away from seizure) epochs [43]. Therefore future work should seek to ascertain whether information capable of guiding surgical strategies can be extracted from background data recorded using iEEG. Additionally, it is necessary to examine what neuroimaging modalities contain the most significant information to infer the ictogenic network.

For our theoretical analysis of the impact of network structure on ictogenicity we quantitatively compared both the canonical model and the full neural mass (Wendling model) [37,46]. Our analysis demonstrated that the canonical model, which is essentially a normal form representation of the SNIC bifurcation present in the neural mass model, is a useful parsimonious model, particularly for the purpose of finding the distribution of node ictogenicity (NI). In fact, we found that the distribution of NI across nodes of a network is almost independent from the particular choice of parameters for sufficiently small networks (N < 30). For larger networks (N = 50) whose topology yields a non-uniform distribution of NI, the two models also return similar predictions of NI without the need for parameter calibration. This implies that NI depends predominantly on the presence of a bifurcation to spiking dynamics and network structure. We have further shown that the computational gain increases nonlinearly with increasing network size and it is therefore significant for networks such as those inferred from iEEG. The reduction of complexity of models to study fundamental mechanisms of epilepsy [9,74,75] or healthy brain dynamics [76,77] is becoming a well-accepted approach, and will be particularly important in computationally intensive applications, such as the study of perturbations to high-dimensional networks.

Although here we have focused on noise-driven models close to a SNIC bifurcation to generate relevant dynamics, previous studies have suggested the use of other bifurcations, such as saddle-node and homoclinic bifurcations to model seizure onset and offset, respectively [75]. In particular, Jirsa et al. [75] used a data-driven approach to identify these bifurcations, under the assumption of a slowly changing control variable moving the model through parameter space. Further work is required to understand to what extent can data reveal which bifurcation underlies the observed dynamical transitions and how different bifurcation mechanisms can influence on the quantification of NI. Patient-specific assessment of the type of bifurcation that best describes the data may lead to further improvements of this modeling framework.

We demonstrated that networks with high degree variance are more likely to seize for relatively smaller coupling strengths, whereas more homogeneous networks reach BNI = 1 within a more confined range of coupling strengths. This provides a potential explanation for observations of increased degree variance in functional networks derived from epilepsy patients, compared to healthy controls [23]. We also uncovered differences in the interplay between global and local spiking generating mechanisms in networks with different topologies: random and small-world networks display a switch-like mechanism for the emergence of spiking dynamics with respect to changes in global coupling, but a more gradual response to removal of nodes. Nodes in random or small-world networks have smaller degree variance, whereas nodes in scale-free or rich-club networks are more heterogeneous in their degree. Therefore, ictogenicity in the latter networks is likely to be concentrated within a few nodes, and thus larger connectivity strengths are required for spiking dynamics to be present in the whole network, which is required here for BNI to be large. These results are in agreement with findings that the critical coupling for the Kuramoto model decreases as the exponent of scale-free networks decreases [78,79]. It is important to note that we assumed all nodes equivalently excitable to focus on the contribution of the network structure to the emergence of spiking dynamics. Future work should consider the potential existence of pathological nodes with higher excitability that may drive the ictogenicity of the network and whose resection may be preferable.

Our analysis suggests that if a brain network under consideration does not have rich-club organization, or if the rich-club were to overlap with eloquent cortex, a resection of a much greater number of nodes would be required. Note that in this context the brain network mapped from iEEG does not correspond to the whole brain, instead it corresponds to a clinically predetermined brain region under investigation as potential surgical targets. Interestingly, our results suggest that a considerable BNI reduction could be attained in most 64 node networks upon removal of 10 nodes at random in all topologies, which is comparable to the average epilepsy resection [62]. This is in agreement with findings that most patients undergoing surgery experience some reduction in seizures, even if they do not achieve seizure freedom [80]. In some cases, the rich-club may comprise non-adjacent nodes making it difficult to resect it through surgery. To tackle these cases, other techniques might be considered such as radiofrequency thermocoagulation [81]. Our approach may be improved by quantitatively assessing predictions for changes in seizures frequency, based on the baseline seizures rates of the individuals.

Interestingly, our findings regarding undirected networks did not extend readily to directed networks. In particular, graph theoretical properties of nodes in directed networks did not correlate with NI, and the effect of node removals was found to be smaller than it would have in “similar” undirected networks. This has important implications for the choice of network representation of the brain used in studies of perturbations. Depending on the data modality under consideration, different approaches should be considered: for modalities that give rise to undirected networks our framework suggests to target nodes according to their eigenvector centrality, whereas if directed networks are derived it is necessary to assess NI using a model (such as the CM). Ultimately, a representation of the brain will be deemed (clinically) useful in the context of our study if it is able to predict the outcome of perturbations. In the current study and our previous work [42] we have demonstrated that useful information is present in an undirected network representation of the brain. However, future work will ascertain whether approaches yielding a directed network may ultimately prove most beneficial.

Supporting information

S1 Fig. NI of different directed network topologies in the CM.

Each row corresponds to a different network topology: (a)-(c) random network; (d)-(f) scale-free network (SF) (γ = 3); (g)-(i) small-world network (SW); (j)-(l) rich-club (RC) network. In the first column the nodes are sorted by their NI, so that NI is monotonically increasing; in the second column NI is sorted by the product of in- and out-degree; and in the third column NI is sorted by the dynamical importance (DI) of the nodes. Note that the correlation between these node measures and NI is not as good as the one found in Fig 5 for undirected networks. The shaded areas in the first column and the dots in the other panels correspond to 10 different network realizations of the same topology (the line represent the mean). Parameter choices are as in Fig 4.

https://doi.org/10.1371/journal.pcbi.1005637.s001

(TIF)

References

  1. 1. Banerjee PN, Filippi D, Hauser WA. The descriptive epidemiology of epilepsy—a review. Epilepsy Res. 2009;85: 31–45. pmid:19369037
  2. 2. Kwan P, Brodie MJ. Early identification of refractory epilepsy. N Engl J Med. 2000;342: 314–9. pmid:10660394
  3. 3. Rosenow F, Lüders H. Presurgical evaluation of epilepsy. Brain. 2001;124: 1683–1700. pmid:11522572
  4. 4. De Tisi J, Bell GS, Peacock JL, McEvoy AW, Harkness WF, Sander JW, et al. The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: A cohort study. Lancet. Elsevier Ltd; 2011;378: 1388–1395. pmid:22000136
  5. 5. Najm I, Jehi L, Palmini A, Gonzalez-Martinez J, Paglioli E, Bingaman W. Temporal patterns and mechanisms of epilepsy surgery failure. Epilepsia. 2013;54: 772–782. pmid:23586531
  6. 6. Richardson MP. New observations may inform seizure models: Very fast and very slow oscillations. Progress in Biophysics and Molecular Biology. 2011. pp. 5–13. pmid:20875832
  7. 7. Richardson MP. Large scale brain models of epilepsy: dynamics meets connectomics. J Neurol Neurosurg Psychiatry. 2012;83: 1238–48. pmid:22917671
  8. 8. Spencer SS. Neural Networks in Human Epilepsy: Evidence of and Implications for Treatment. Epilepsia. 2002;43: 219–227. pmid:11906505
  9. 9. Terry JR, Benjamin O, Richardson MP. Seizure generation: The role of nodes and networks. Epilepsia. 2012;53: 166–169. pmid:22709380
  10. 10. Shah D, Blockx I, Keliris GA, Kara F, Jonckers E, Verhoye M, et al. Cholinergic and serotonergic modulations differentially affect large-scale functional networks in the mouse brain. Brain Struct Funct. 2015; pmid:26195064
  11. 11. Shen K, Hutchison RM, Bezgin G, Everling S, McIntosh AR. Network Structure Shapes Spontaneous Functional Connectivity Dynamics. J Neurosci. 2015;35: 5579–5588. pmid:25855174
  12. 12. Varela C. Thalamic neuromodulation and its implications for executive networks. Front Neural Circuits. 2014;8: 69. pmid:25009467
  13. 13. Bressler SL, Menon V. Large-scale brain networks in cognition: emerging methods and principles. Trends in Cognitive Sciences. 2010. pp. 277–290. pmid:20493761
  14. 14. Bullmore E and OS . Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci. 2009;10: 186–198. pmid:19190637
  15. 15. Bullmore E, Sporns O. The economy of brain network organization. Nat Rev Neurosci. 2012;13: 336–349. pmid:22498897
  16. 16. Fornito A, Zalesky A, Breakspear M. The connectomics of brain disorders. Nat Rev Neurosci. Nature Publishing Group; 2015;16: 159–172. pmid:25697159
  17. 17. Park HJ, Friston K. Structural and functional brain networks: from connections to cognition. Science. 2013;342: 1238411. pmid:24179229
  18. 18. Sporns O, Chialvo DR, Kaiser M, Hilgetag CC. Organization, development and function of complex brain networks. Trends Cogn Sci. 2004;8: 418–425. pmid:15350243
  19. 19. Sporns O. Networks of the Brain. MIT press; 2011.
  20. 20. Bernhardt BC, Chen Z, He Y, Evans AC, Bernasconi N. Graph-theoretical analysis reveals disrupted small-world organization of cortical thickness correlation networks in temporal lobe epilepsy. Cereb Cortex. 2011;21: 2147–2157. pmid:21330467
  21. 21. Bonilha L, Nesland T, Martz GU, Joseph JE, Spampinato M V., 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–909. pmid:22764263
  22. 22. Chavez M, Valencia M, Latora V, Martinerie J. Complex Networks: New Trends for the Analysis of Brain Connectivity. Int J Bifurc Chaos. 2010;20: 1677–1686.
  23. 23. Chowdhury FA, Woldman W, FitzGerald THB, Elwes RDC, Nashef L, Terry JR, et al. Revealing a brain network endophenotype in families with idiopathic generalised epilepsy. PLoS One. 2014;9. pmid:25302690
  24. 24. Horstmann MT, Bialonski S, Noennig N, Mai H, Prusseit J, Wellmer J, et al. State dependent properties of epileptic brain networks: Comparative graph-theoretical analyses of simultaneously recorded EEG and MEG. Clin Neurophysiol. 2010;121: 172–185. pmid:20045375
  25. 25. Quraan MA, McCormick C, Cohn M, Valiante TA, McAndrews MP. Altered Resting State Brain Dynamics in Temporal Lobe Epilepsy Can Be Observed in Spectral Power, Functional Connectivity and Graph Theory Metrics. PLoS One. 2013;8. pmid:23922658
  26. 26. Vaessen MJ, Jansen JFA, Vlooswijk MCG, Hofman PAM, Majoie HJM, Aldenkamp AP, et al. White matter network abnormalities are associated with cognitive decline in chronic epilepsy. Cereb Cortex. 2012;22: 2139–2147. pmid:22038907
  27. 27. Boccaletti S, Latora V, Moreno Y, Chavez M, Hwang DU. Complex networks: Structure and dynamics. Phys Rep. 2006;424: 175–308.
  28. 28. Strogatz SH. Exploring complex networks. Nature. 2001;410: 268–276. pmid:11258382
  29. 29. Kopell NJ, Gritton HJ, Whittington MA, Kramer MA. Beyond the connectome: The dynome. Neuron. 2014. pp. 1319–1328. pmid:25233314
  30. 30. Honey CJ, Kötter R, Breakspear M, Sporns O. Network structure of cerebral cortex shapes functional connectivity on multiple time scales. Proc Natl Acad Sci U S A. 2007;104: 10240–10245. pmid:17548818
  31. 31. Benjamin O, Fitzgerald TH, Ashwin P, Tsaneva-Atanasova K, Chowdhury F, Richardson MP, et al. A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy. J Math Neurosci. 2012;2: 1. pmid:22657571
  32. 32. Blenkinsop A, Valentin A, Richardson MP, Terry JR. The dynamic evolution of focal-onset epilepsies—combining theoretical and clinical observations. Eur J Neurosci. 2012;36: 2188–2200. pmid:22805064
  33. 33. Goodfellow M, Schindler K, Baier G. Intermittent spike-wave dynamics in a heterogeneous, spatially extended neural mass model. Neuroimage. Elsevier Inc.; 2011;55: 920–932. pmid:21195779
  34. 34. Marten F, Rodrigues S, Suffczynski P, Richardson MP, Terry JR. Derivation and analysis of an ordinary differential equation mean-field model for studying clinically recorded epilepsy dynamics. Phys Rev E—Stat Nonlinear, Soft Matter Phys. 2009;79. pmid:19391782
  35. 35. Schmidt H, Petkov G, Richardson MP, Terry JR. Dynamics on Networks: The Role of Local Dynamics and Global Networks on the Emergence of Hypersynchronous Neural Activity. PLoS Comput Biol. 2014;10. pmid:25393751
  36. 36. Stam CJ. Modern network science of neurological disorders. Nat Rev Neurosci. 2014;15: 683–695. pmid:25186238
  37. 37. 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–1508. pmid:12028360
  38. 38. Schmidt H, Woldman W, Goodfellow M, Chowdhury FA, Koutroumanidis M, Jewell S, et al. A computational biomarker of idiopathic generalized epilepsy from resting state EEG. Epilepsia. 2016. pmid:27501083
  39. 39. Alstott J, Breakspear M, Hagmann P, Cammoun L, Sporns O. Modeling the impact of lesions in the human brain. PLoS Comput Biol. 2009;5. pmid:19521503
  40. 40. Honey CJ, Sporns O. Dynamical consequences of lesions in cortical networks. Hum Brain Mapp. 2008;29: 802–809. pmid:18438885
  41. 41. Loh M, Rolls ET, Deco G. A dynamical systems hypothesis of schizophrenia. PLoS Comput Biol. 2007;3: 2255–2265. pmid:17997599
  42. 42. Goodfellow M, Rummel C, Abela E, Richardson MP, Schindler K, Terry JR. Estimation of brain network ictogenicity predicts outcome from epilepsy surgery. Sci Rep. 2016;6: 29215. pmid:27384316
  43. 43. Sinha N, Dauwels J, Kaiser M, Cash SS, Westover MB, Wang Y, et al. Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. Brain. pmid:28011454
  44. 44. Wendling F, Bellanger JJ, Bartolomei F, Chauvel P. Relevance of nonlinear lumped-parameter models in the analysis of depth-EEG epileptic signals. Biol Cybern. 2000;83: 367–378. pmid:11039701
  45. 45. Ermentrout GB, Kopell N. Parabolic Bursting in an Excitable System Coupled with a Slow Oscillation. SIAM J Appl Math. 1986;46: 233–253.
  46. 46. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995;73: 357–366. pmid:7578475
  47. 47. Grimbert F, Faugeras O. Bifurcation analysis of Jansen’s neural mass model. Neural Comput. 2006;18: 3052–3068. pmid:17052158
  48. 48. Goodfellow M, Schindler K, Baier G. Self-organised transients in a neural mass model of epileptogenic tissue dynamics. Neuroimage. 2012;59: 2644–2660. pmid:21945465
  49. 49. Laing CR. Derivation of a neural field model from a network of theta neurons. Phys Rev E—Stat Nonlinear, Soft Matter Phys. 2014;90: 1–5. pmid:25122239
  50. 50. McKennoch S, Voegtlin T, Bushnell L. Spike-timing error backpropagation in theta neuron networks. Neural Comput. 2009;21: 9–45. pmid:19431278
  51. 51. Börgers C, Kopell N. Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity. Neural Comput. 2003;15: 509–538. pmid:12620157
  52. 52. Bédard C, Kröger H, Destexhe A. Model of low-pass filtering of local field potentials in brain tissue. Phys Rev E—Stat Nonlinear, Soft Matter Phys. 2006;73: 1–15. pmid:16802971
  53. 53. Mišić B, Sporns O, McIntosh AR. Communication Efficiency and Congestion of Signal Traffic in Large-Scale Brain Networks. PLoS Comput Biol. 2014;10. pmid:24415931
  54. 54. Newman MEJ. The Structure and Function of Complex Networks. Soc Ind Appl Math Rev. 2003;45: 167–256.
  55. 55. Newman MEJ. The mathematics of networks. New Palgrave Encycl Econ. 2007;2: 1–12.
  56. 56. Rubinov M, Sporns O. Complex network measures of brain connectivity: Uses and interpretations. Neuroimage. Elsevier Inc.; 2010;52: 1059–1069. pmid:19819337
  57. 57. Alvarez-Socorro AJ, Herrera-Almarza GC, González-Díaz LA. Eigencentrality based on dissimilarity measures reveals central nodes in complex networks. Sci Rep. Nature Publishing Group; 2015;5: 17095. pmid:26603652
  58. 58. Restrepo JG, Ott E, Hunt BR. Characterizing the dynamical importance of network nodes and links. Phys Rev Lett. 2006;97: 1–4. pmid:17026366
  59. 59. Petkov G, Goodfellow M, Richardson MP, Terry JR. A critical role for network structure in seizure onset: A computational modeling approach. Front Neurol. 2014;5: 1–7.
  60. 60. Carterette B. On rank correlation and the distance between rankings. Proc 32nd Int Conf Res Dev Inf Retr. 2009; 436.
  61. 61. Kumar R, Vassilvitskii S. Generalized distances between rankings. Proc 19th Int Conf World wide web—WWW ‘10. 2010; 571.
  62. 62. Rummel C, Abela E, Andrzejak RG, Hauf M, Pollo C, Müller M, et al. Resected brain tissue, seizure onset zone and quantitative EEG measures: Towards prediction of post-surgical seizure control. PLoS One. 2015;10. pmid:26513359
  63. 63. Rummel C, Abela E, Müller M, Hauf M, Scheidegger O, Wiest R, et al. Uniform approach to linear and nonlinear interrelation patterns in multivariate time series. Phys Rev E—Stat Nonlinear, Soft Matter Phys. 2011;83. pmid:21797469
  64. 64. Rummel C, Müller M, Baier G, Amor F, Schindler K. Analyzing spatio-temporal patterns of genuine cross-correlations. J Neurosci Methods. 2010;191: 94–100. pmid:20566351
  65. 65. Colizza V, Flammini A, Serrano MA, Vespignani A. Detecting rich-club ordering in complex networks. Nature. 2006;2: 110–115.
  66. 66. Opsahl T, Colizza V, Panzarasa P, Ramasco JJ. Prominence and control: The weighted rich-club effect. Phys Rev Lett. 2008;101: 1–4. pmid:18999722
  67. 67. van den Heuvel MP, Sporns O. Rich-Club Organization of the Human Connectome. J Neurosci. 2011;31: 15775–15786. pmid:22049421
  68. 68. Sinha N, Dauwels J, Wang Y, Cash SS, Taylor PN. An in silico approach for pre-surgical evaluation of an epileptic cortex. 2014 36th Annu Int Conf IEEE Eng Med Biol Soc EMBC 2014. 2014; 4884–4887. pmid:25571086
  69. 69. Stam CJ. Epilepsy: What can we Learn from Modern Network Theories. Epileptologie. 2016;33: 38–43.
  70. 70. Varotto G, Tassi L, Franceschetti S, Spreafico R, Panzica F. Epileptogenic networks of type II focal cortical dysplasia: A stereo-EEG study. Neuroimage. 2012;61: 591–598. pmid:22510255
  71. 71. Wilke C, Worrell G, He B. Graph analysis of epileptogenic networks in human partial epilepsy. Epilepsia. 2011;52: 84–93. pmid:21126244
  72. 72. Rosso OA, Mendes A, Rostas JA, Hunter M, Moscato P. Distinguishing childhood absence epilepsy patients from controls by the analysis of their background brain electrical activity. J Neurosci Methods. 2009;177: 461–468. pmid:19013193
  73. 73. Ritter P, Schirner M, McIntosh AR, Jirsa VK. The Virtual Brain Integrates Computational Modeling and Multimodal Neuroimaging. Brain Connect. 2013;3: 121–145. pmid:23442172
  74. 74. Goodfellow M, Glendinning P. Mechanisms of intermittent state transitions in a coupled heterogeneous oscillator model of epilepsy. J Math Neurosci. 2013;3: 17. pmid:23945016
  75. 75. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014;137: 2210–2230. pmid:24919973
  76. 76. Deco G, McIntosh AR, Shen K, Hutchison RM, Menon RS, Everling S, et al. Identification of Optimal Structural Connectivity Using Functional Connectivity and Neural Modeling. J Neurosci. 2014;34: 7910–7916. pmid:24899713
  77. 77. Freyer F, Roberts JA, Ritter P, Breakspear M. A Canonical Model of Multistability and Scale-Invariance in Biological Systems. PLoS Comput Biol. 2012;8. pmid:22912567
  78. 78. Lopes MA, Lopes EM, Yoon S, Mendes JFF, Goltsev A V. Synchronization in the random field Kuramoto model on complex networks. 2016; 1–7. Available: http://arxiv.org/abs/1605.04733
  79. 79. Yoon S, Sorbaro Sindaci M, Goltsev A V, Mendes JF. Critical behavior of the relaxation rate, the susceptibility, and a pair correlation function in the Kuramoto model on scale-free networks. Phys Rev E Stat Nonlin Soft Matter Phys. 2015;91: 32814. pmid:25871164
  80. 80. Wieser HG, Blume WT, Fish D, Goldensohn E, Hufnagel A, King D, et al. Proposal for a new classification of outcome with respect to epileptic seizures following epilepsy surgery. Epilepsia. 2001;42: 282–286. pmid:11240604
  81. 81. Cossu M, Fuschillo D, Casaceli G, Pelliccia V, Castana L, Mai R, et al. Stereoelectroencephalography-guided radiofrequency thermocoagulation in the epileptogenic zone: a retrospective study on 89 cases. J Neurosurg. 2015;123: 1358–1367. pmid:26090841