Treatment effects in epilepsy: a mathematical framework for understanding response over time

Epilepsy is a neurological disorder characterized by recurrent seizures, affecting over 65 million people worldwide. Treatment typically commences with the use of anti-seizure medications, including both mono- and poly-therapy. Should these fail, more invasive therapies such as surgery, electrical stimulation and focal drug delivery are often considered in an attempt to render the person seizure free. Although a significant portion ultimately benefit from these treatment options, treatment responses often fluctuate over time. The physiological mechanisms underlying these temporal variations are poorly understood, making prognosis a significant challenge when treating epilepsy. Here we use a dynamic network model of seizure transition to understand how seizure propensity may vary over time as a consequence of changes in excitability. Through computer simulations, we explore the relationship between the impact of treatment on dynamic network properties and their vulnerability over time that permit a return to states of high seizure propensity. For small networks we show vulnerability can be fully characterised by the size of the first transitive component (FTC). For larger networks, we find measures of network efficiency, incoherence and heterogeneity (degree variance) correlate with robustness of networks to increasing excitability. These results provide a set of potential prognostic markers for therapeutic interventions in epilepsy. Such markers could be used to support the development of personalized treatment strategies, ultimately contributing to understanding of long-term seizure freedom.

1 Introduction

The response to treatment in epilepsy–such as anti-seizure medication (ASM), neurostimulation, or surgery - often fluctuates over time. Most clinical studies in this context have examined this with respect to the overall long-term probability of seizure freedom for people with epilepsy (e.g., probability of no seizures for a period of at least 1 year) (Brodie et al., 2012; Chen et al., 2018). On shorter time-scales, perhaps the most prominent example of a transient, declining change in treatment outcomes is the so-called “honeymoon effect”. This is broadly characterized by a period of significant reduction in seizure frequency following the intervention: the “honeymoon” period, which can last from a couple of weeks to several months. By definition, the honeymoon period is followed by an increase in seizure frequency, sometimes to levels at least as high as those prior to the intervention (Boggs et al., 2000). Understanding mechanisms that contribute to this phenomenon are therefore critical for improving outcomes for people with epilepsy.

For anti-seizure medications (ASMs), the honeymoon effect (also referred to as acquired tolerance) is most commonly observed in benzodiazepines. As a result they are generally considered unsuitable for use as long-term treatments (Riss et al., 2008). However, people with epilepsy are known to develop tolerance to a wider range of medications (Löscher and Schmidt, 2006; Abou-Khalil et al., 2008). One study published in 2000 saw that as many of 22 out of 80 patients experienced a return to at least a baseline seizure frequency after an initial positive response to medication (Boggs et al., 2000). Animal studies have suggested seizure type may also affect the likelihood of tolerance to certain ASMs, however little supporting evidence exists in humans (Löscher et al., 1996; Löscher and Schmidt, 2006). Further, there is the phenomenon of “cross-tolerance”. Here, tolerance to one medication may lead to tolerance to another (Abou-Khalil et al., 2008). This observation could support the notion that tolerance emerges due to adaptation of the underlying mechanisms of seizure generation and propagation.

In the context of epilepsy surgery, patients undergoing resective surgery often experience a relapse in seizures following an apparently successful period of remission (De Tisi et al., 2011; Petrik et al., 2021). Rates of seizure recurrence and surgical success rates are highly heterogeneous, depending both on the type of seizures experienced, the presence (or otherwise) of visible brain lesions and the nature of the surgery performed (Téllez-Zenteno et al., 2005). Contributing factors can include incomplete resection of the epileptogenic zone, as well as the emergence of new or previously-undetected epileptogenic networks. In this latter scenario, seizures of a different nature to those experienced prior to surgery are commonplace (Petrik et al., 2021). Onset of epilepsy during the first year of life has been associated with late seizure relapse following apparently successful surgery (Petrik et al., 2021), supporting the idea that some brains may have a more ingrained tendency to ictogenicity. Other studies have also shown that seizure recurrence within 6 months of surgery is associated with people whose seizures began earlier in life (Goellner et al., 2013).

Mathematical modelling approaches constitute a powerful tool to investigate how putative physiological factors could drive changes in seizure propensity (Touboul et al., 2011; Goodfellow et al., 2012; Petkov et al., 2014; Jirsa et al., 2017; Proix et al., 2017; Sinha et al., 2017; Olmi et al., 2019; Junges et al., 2020; Gerster et al., 2021). Typically these models endow regions within a network structure with a specific mechanism describing a process of interest (e.g., a transition into a seizure state). Further, the network structures are often inferred or estimated directly from imaging or neurophysiological data (Wang et al., 2014; Chiarion et al., 2023).

Given the relationship between hyperexcitability and seizures, these models typically consider the role of excitability for creating the conditions (or likelihood) for synchronised activity (i.e., a proxy for seizure-like activity). At the macroscale, it is important to disambiguate “excitability” from the activity of excitatory neurons. In this context, neural mass models or mean-field models are typically used to consider the impact on emergent excitability of the balance between excitatory and inhibitory populations. This overall balance has been explored in general terms with respect to seizure generation and termination, and more specifically in the context of seizure-likelihood during treatment with ASM (Kramer et al., 2012; Jiruska et al., 2013; Jirsa et al., 2014; Meisel et al., 2015).

We have previously used dynamic network models of seizure transition to simulate and predict the effects of treatment strategies (Woldman et al., 2019; Junges et al., 2020). In these model frameworks, seizure propensity is fundamentally impacted by a combination of local brain excitability and network configuration. Some key aspects of dynamic brain networks (as represented in these models) include seizure propensity, as well as robustness to change, i.e., how robust a state of low seizure propensity is to potential future changes in network reconfiguration and/or local excitability.

Three scenarios for response to treatment are outlined in Figure 1. The first (Scenario A) describes an effective treatment regime in which seizure frequency is reduced to a sustained extent following the onset of treatment. While a temporary reduction in seizure count following treatment may be attributed to natural variation in seizure rates (Scenario B), observations of initial success of treatment with a later return of seizures suggest some mechanism by which the brain adapts towards increased seizure risk (Scenario C).

Figure 1

Figure 1. Some potential responses to epilepsy treatment over an arbitrary time scale and an arbitrary scale of seizure risk. Seizure risk is represented on the y-axis and time on the x-axis. Time of treatment onset is marked with a vertical dashed purple line, and a threshold for seizure freedom is marked with a horizontal dashed grey line. In scenario (A), treatment reduces seizure risk below the seizure freedom threshold and risk remains low. In scenario (B), seizure risk is naturally fluctuating over time. Treatment does not affect the seizure risk, but fluctuations align with treatment onset by chance. In scenario (C), treatment initially reduces seizure risk to the same extent as in scenario (A), but seizure risk gradually returns to the initial pre-treatment level. This corresponds to the honeymoon effect and is the scenario we propose a framework to describe in this study. Note that in practice scenarios (B,C) may be difficult to disambiguate: the seizure freedom threshold is crossed at the same times in both cases.

The trends in seizure propensity observed over extended periods of time post treatment suggest a diversity of neurophysiological mechanisms may contribute to the return of seizures in people with epilepsy. Whilst a comprehensive exploration of such mechanisms using a mathematical modelling framework is unfeasible, we can use this framework to estimate the effects of general neurophysiological alterations in seizure propensity. For example, exploring how a brain associated with certain model parameters responds to a homogeneous increase in excitability can help us to better understand, and potentially predict, seizure return. Such a hypothesis might be more representative of some cases than others (e.g., generalized versus focal epilepsies), but can nonetheless help identifying features strongly associated to changes in seizure propensity, as well as serve as a basic framework that can be extended to more specific scenarios.

In this study, we systematically explore the relationship between properties of brain networks and the response of these networks to a drive towards increased seizure propensity (considered as an increase in overall excitability). These network properties range in complexity from simple, such as average clustering coefficient or efficiency, to more complex features such as the first transitive component or trophic incoherence. They capture a variety of features of networks which may be relevant for prognosis, for example, the combined effect of network directionality and small-worldness. A drive towards increased seizure propensity may represent the combined effect of neurophysiological mechanisms that underpin transient responses to treatment. Within our proposed framework, we assess robustness to these mechanisms through quantification of the speed and extent that seizure propensity increases when a perturbation is applied, and how properties of the network can impact these changes.

2 Materials and methods

We consider a phenomenological model of seizure transition that permits the existence of two states, reflecting seizure-like behaviour and background activity respectively (Kalitzin et al., 2010; Benjamin et al., 2012). The model is a modified version of the normal form of the subcritical Hopf bifurcation, which incorporates a slowly varying time-dependent variable λ reflecting the level of “excitability” within a brain region. This type of model has been utilised for epilepsy in a variety of contexts (Kalitzin et al., 2010; Benjamin et al., 2012; Hebbink et al., 2017; Junges et al., 2020). Within a single node, the dynamic evolution is defined by:

$\frac{\text{d}z}{\text{d}t}=z\left(\lambda -1+\text{i}\omega +2{\left|z\right|}^2-{\left|z\right|}^4\right)+\alpha \text{dW},(1)$

$\tau \frac{\text{d}\lambda }{\text{d}t}={\lambda }_0-\lambda -{\left|z\right|}^2.(2)$

In this formulation, node activity z is a complex variable x + iy, such that Re(z) = x corresponds to the simulated EEG electrode activity. ω is the frequency of the limit cycle of this system, and can be tuned so that this frequency corresponds to that seen during seizure-like activity. λ is the time-dependent node excitability and λ₀ the constant baseline level of excitability. τ is a real constant modulating the rate at which a node transitions from the seizure-like to the non seizure-like state. α reflects the impact of dynamic inputs received by regions of the brain (nodes) that are not explicitly accounted for within the model dynamics. This noise drives transitions to the limit cycle. Under an Euler-Maruyama scheme, dW draws a value from the uniform distribution bounded by [0, $\sqrt{dt}$ ] at each time-step. Typical values for these model parameters can be found within Table 1.

Table 1

Table 1. Parameters of the modified subcritical Hopf model used in this paper (Junges et al (2020)).

This system is deterministic when α = 0. Within the physically permissible region of λ ∈ [0, 1] and $\left|z\right|\ge 0$, the dynamics of z are characterised by three steady-state solutions. The first is a stable limit cycle at ${\left|z\right|}^2=1+\sqrt{\lambda }$. This represents seizure-like behaviour. The second is a fixed point which exists at z = 0 corresponding to background (inter-ictal) behaviour. Separating these two stable solutions is an unstable limit cycle at ${\left|z\right|}^2=1-\sqrt{\lambda }$. These steady state solutions are shown in Figure 2.

Figure 2

Figure 2. (A): bifurcation diagram of |z| in λ. Stable steady states are indicated as solid lines, and unstable steady states as dashed lines. The red point at (1, 0) represents the Hopf Bifurcation. Two limit points are marked in green at (0, 1) and (0, −1). (B): The nullclines of λ at different values of λ₀, which is parabolic with a peak at $(\left|z\right|=0,\lambda ={\lambda }_0)$. Grey dashed lines indicate the maximum and minimum values of λ₀. When $\left|z\right|=0$ in the interictal state, λ converges to λ₀. In the ictal state, $\left|z\right|$ takes on larger values and consequently λ converges to lower values.

Eq. 2 describes the λ dynamics, and has a single stable steady-state solution in which $\lambda ={\lambda }_0-{\left|z\right|}^2$, illustrated in Figure 2. When z is close to its fixed point, in the inter-ictal state, |z|² is very small and λ therefore tends to λ₀; in the case where $\left|z\right|$ is large, λ tends towards values less than λ₀, decreasing the excitability and forcing a natural return of the system towards the fixed point. The time-scale of this return is determined by the slow variable τ. There exists a steady state in both λ and z at z = 0 and λ = λ₀. As λ₀ increases the size of the basin of attraction of the stable limit cycle decreases, such that the distance between the z = 0 and the unstable limit cycle decreases. This makes it easier for the system to undergo transition to the seizure-like state.

The dynamics in the non-deterministic system (α > 0) are shown for a single node in Figure 3. The system stays near the initial condition of the stable fixed point (z = 0, λ = λ₀) unless the influence of the noise is strong enough that the trajectory crosses the unstable ‘separatrix’ solution causing it to approach the ictal state. Increased activity z results in a drop in the excitability λ - the key mechanism which returns the system from the seizure-like state to the inter-ictal state.

Figure 3

Figure 3. Dynamics of the modified subcritical Hopf model for a single node. (A) Trajectory in phase space for a single node. The direction of flow is anti-clockwise. (B) Simulated single-node EEG activity Re(z) over a timescale of T = 200s. (C) Variation in the slow excitability variable λ.

To appropriately reflect brain network activity, we extend the model consider interacting node dynamics with a network. Here we describe network dynamics as a system of coupled time-dependent stochastic differential equations,

$\text{d}{z}_i=\left(z_i\left({\lambda }_i-1+\text{i}\omega +2{\left|z_i\right|}^2-{\left|z_i\right|}^4\right)+\frac{1}{N}{\displaystyle \sum _{j=1}^N}{M}_{ij}g\left(z_i,z_j\right)\right)\text{d}t+\alpha {\text{dW}}_i,(3)$

$\tau \text{d}{\lambda }_i=\left({\lambda }_0-{\lambda }_i-|z_i{|}^2\right)\text{d}t,(4)$

$g\left(z_i,z_j\right)=\beta \left(z_j-z_i\right).(5)$

N is the total number of nodes within the network such that i = 1, 2, … , N. M_ij is the adjacency matrix categorising the edges within a network. We consider here directed networks without self-loops, so M_ij is allowed to be asymmetric with all diagonal entries set to 0. For the purposes of the present study we consider only binary networks such that entries of M_ij are 1 or 0, corresponding to the existence or nonexistence of an edge. g (z_j, z_i) is the diffusive coupling function where β is a real constant specifying the strength of the system coupling. Through this choice of coupling, interactions between nodes can additionally influence transitions from the steady-state to the stable limit cycle and vice versa.

2.1 Quantifying seizure propensity

In order to quantify the likelihood of the system transitioning into the seizure-like state we utilise the concept of Brain Network Ictogenicity (BNI) (Petkov et al., 2014; Lopes et al., 2018). Since its introduction, BNI has been used in many works as a measure of network seizure propensity. This includes as a way to quantify the change in seizure propensity after network nodes are removed or network edges are perturbed (Sinha et al., 2017; Goodfellow et al., 2016; Lopes et al., 2017; Lopes et al., 2019; Lopes et al., 2020). It evaluates the likelihood of simulated networks to transit into a seizure state, calculated as

$BNI={\displaystyle \sum _{i=1}^{T/dt}}\frac{f\left(m\right)}{TN}.(6)$

T is here the total simulated time-span and dt the time-step when system trajectories are evaluated, using a first-order Euler-Maruyama scheme. To ensure consistency between network simulations, the initial state for all nodes is the background (steady) state. m is the number of nodes found within the seizure-like state during a given time step, such that

$f\left(m\right)=\left\{\begin{array}{cc}0\hfill & m<2\hfill \\ m\hfill & \text{otherwise}\hfill \end{array}\right.(7)$

A node is considered as being in a seizure-like state if ${\left|z\right|}^2>0.5$. We consider the system is in a seizure when at least 2 nodes are in the seizure-like state (m ≥ 2). As we are primarily interested in seizure transitions due to synchronization, cases where m < 2 are discarded to minimise the effect of spontaneous transitions into the stable limit cycle. Consequently, if no two or more nodes simultaneously enter the seizure-like state across an entire simulation then BNI = 0, whereas if all network nodes remain in the high activity state throughout the simulation then BNI = 1. In order to exclude the possibility that any network behaviour is unique to a certain range of the coupling parameter β, BNI is averaged over β ∈ [0, 6] such that it includes behaviour from both the strong and weak coupling regimes (Junges et al., 2020). To account for the impact of stochasticity, the BNI is then re-averaged over several realisations of noise. Timescales of T = 500s or T = 2000s were chosen for each simulation with a time-step of dt = 0.0005s with 5 realisations implemented for each calculation.

In order to reproduce the increase in seizure propensity observed in the honeymoon effect, we consider a gradual increase in the baseline excitability of the network, λ₀. It would be possible, in this framework, to consider higher-complexity perturbations to the system which correlate with increased seizure propensity, such as specific alterations to network topology or connection strength, however in this work we focus on excitability, the propensity for a system to attain high activity states represented by λ₀, as a key feature which directly determines seizure propensity. Figure 4 shows the trajectories in BNI of two 20-node networks as a function of λ₀. We consider both extremes of network connectivity - a fully-connected network, in which all possible edges are present, and a completely disconnected network, in which no edges are present and hence all nodes are isolated and independent of each other. For values of λ₀ < 0.6, BNI was uniformly 0 for both networks, such that no nodes concurrently entered the seizure-like state at any given point during the simulated time-span. The completely disconnected network has a much more gradual and less steep increase in BNI than the fully connected network.

Figure 4

Figure 4. Trajectories of BNI with increasing λ₀ for the fully-connected and completely disconnected 20-node networks.

We quantify the increase of network BNI under increasing λ₀ by two metrics. The area under the curve (AUC) quantifies the overall increase in seizure propensity upon increase of λ₀. We define the “quartile distance”,

$QD={\lambda }_0{|}_{\text{BNI}=0.75}-{\lambda }_0{|}_{\text{BNI}=0.25},(8)$

such that it is the increase in λ₀ required for network BNI to increase from 0.25 to 0.75. This measure quantifies the rate of increase of network trajectories in seizure propensity as a function of λ₀.

2.2 Network features

In order to quantify the network characteristics associated to seizure propensity we use a number of features of binary directed networks which may have an effect on the dynamic properties of the model system. We limit our analysis to five features, though the framework proposed may be applied to any network feature that is believed to be relevant to seizure generation. We choose measures that vary in their complexity and are designed to capture a variety of network properties. We initially consider the first transitive component, a measure which identifies strongly-connected regions of the brain which may be considered “drivers” of activity, known to have a strong impact on the genesis of seizure activity (Benjamin et al., 2012). To assess the phenomenon of directionality at a more global level than FTC, we also consider the trophic incoherence, which measures the directionality of flow in the edges of a network. The presence of “driver” and “responder” regions is known to impact how different nodes and subgraphs influence each other towards similar behaviour (e.g., seizure-like or inter-ictal states) (Terry et al., 2012; Junges et al., 2020). We further consider efficiency and mean clustering coefficient, two global measures which capture respectively global and (averaged) local characteristics of small-worldness, known to be a significant property of complex networks in nature, including of functional brain networks (Bassett and Bullmore, 2006). Finally, we draw on degree variance as a measure of heterogeneity in the graph topology. Efficiency, clustering coefficient and degree variance have previously been shown to associate with epilepsy diagnosis (Faiman et al., 2021), while trophic incoherence is a relevant feature of information flow in networks (Johnson et al., 2014), making these suitable candidates for analysis of adaptations in functional brain connectivity and seizure prognosis.

2.2.1 First transitive component

The First Transitive Component (FTC) is a measure of network connectivity introduced in the context of modelling seizures by Benjamin et al. (Benjamin et al., 2012). The FTC is defined to be a set of nodes consisting of all regions of the network which are strongly-connected (there is a directed path in each direction between any pair of nodes) but receive no information from elsewhere in the network: these regions “drive” but do not “respond” to the behaviour of other nodes in the network. More formally, for a directed network of size N with nodes {n₁, n₂, … , n_N}, consider each distinct pair of nodes {n_i, n_j}. We state that n_i ≪ n_j if there exists a directed path from n_i to n_j. Likewise, n_j ≪ n_i if there exists a directed path from n_j to n_i. A node is a member of the FTC if every node from which it is reachable is reachable from it in return. The FTC is therefore the set of nodes n_j such that any n_i for which n_i ≪ n_j also satisfies n_j ≪ n_i. It has been shown that the FTC is a predictor for the network escape times (another measure of seizure propensity which can act as a proxy for the BNI) of small networks (N ≤ 4) (Benjamin et al., 2012). Figure 5 depicts all 13 non-isomorphic connected 3 node networks grouped by the size of FTC, which we denote as n.

Figure 5

Figure 5. All connected directed networks of size 3 and their first transitive components. FTCs of size n = 1,2, and 3 are highlighted in purple, green and blue respectively. The first network for n = 2 is coloured a different green to highlight that the FTC comprises two separate components of the graph, which are not connected by an edge. This leads to distinct behaviour from other networks with FTC of size 2.

2.2.2 Trophic incoherence

Trophic incoherence is a quantification of the extent to which the overall flow of information in a directed network follows a distinct hierarchy or direction (Johnson et al., 2014). A low trophic incoherence describes a network in which information predominantly flows in a single direction, such as in traditional descriptions of food webs. A demonstration of this property is shown in Figure 6. The relative position of a node within the hierarchy of a directed network is quantified as the “trophic level”. In the context of modelling seizure dynamics, high trophic levels may describe “responder” regions which predominantly receive information from the rest of the network, while low trophic levels would correspond to “driver” regions which predominantly send information to other nodes in the network. Trophic incoherence, then, can be described in this context as a measure of the extent to which different regions of the brain share an equal role in the “driver” and “responder” activity in the network. The revised notion of trophic incoherence used here was introduced by MacKay et al. (MacKay et al., 2020). Since the networks under consideration are unweighted, we will consider the definition in the case where all edges are assigned weight 1.

Figure 6

Figure 6. An example of a network with a low trophic incoherence. Information flows in a clear direction, from “source” nodes (shown in red), via intermediates (shown in magenta), to “sink” nodes (shown in blue).

In the uniformly-weighted case, the imbalance of a node is

$v_i=d_i^{in}-d_i^{\mathit{out}}(9)$

and the total weight of the node is

$u_i=d_i^{in}+d_i^{\mathit{out}}(10)$

where $d_i^{in}$ and $d_i^{\mathit{out}}$ refer to the in-degree and out-degree of node i, respectively.

The trophic level is then defined as the vector solution h to the system of equations

$\mathrm{\Lambda }h=v(11)$

where Λ is the graph Laplacian operator

${\left(\mathrm{\Lambda }h\right)}_j=u_j{h}_j-{\displaystyle \sum _{i\in V}}\left(A_{ij}+A_{ji}\right)h_i(12)$

Finally, then, the trophic incoherence for a uniformly-weighted graph is

$F_0=\frac{1}{\left|E\right|}{\displaystyle \sum _{\left(i,j\right)\in E}}{\left(h_i-h_j-1\right)}^2(13)$

2.2.3 Efficiency

The concept of network efficiency quantifies the ability of nodes in a network to communicate with each other. A highly efficient network contains short paths in both directions between any chosen pair of nodes, while a network with low efficiency contains pairs of nodes which can only communicate through long paths across the network, and may include nodes which are not joined by any path of directed edges, at least in one direction. Formally, the (global) efficiency (Latora and Marchiori, 2001) of a directed network is defined as

$\mathbf{E}\left(G\right)=\frac{1}{N\left(N-1\right)}{\displaystyle \sum _{i\ne j}}\frac{1}{d_{ij}}(14)$

where d_ij is the shortest path length from node i to node j, and its inverse thus describes the efficiency of communication between the two nodes. The global efficiency is therefore an average of efficiencies between all (ordered) pairs of nodes in the directed network. Efficiency is an important marker in functional brain connectivity networks, as it represents the ability of regions of the brain to share information, and describes the global behaviour of small-world networks (Latora and Marchiori, 2001; Rubinov and Sporns, 2010).

2.2.4 Mean local clustering coefficient

Clustering is the property of a network in which neighbouring nodes share mutual neighbours, forming densely-connected groups of nodes within the graph structure. For a binary directed network, the notion of mean clustering coefficient is defined as follows (Fagiolo, 2007):

$C\left(G\right)=\frac{1}{N}{\displaystyle \sum _{i\in V}}\frac{{\left[{\left(A+A^T\right)}^3\right]}_{ii}}{2\left[\left(d_i^{in}+d_i^{\mathit{out}}\right)\left(d_i^{in}+d_i^{\mathit{out}}-1\right)-2d_i^{\leftrightarrow }\right]}(15)$

Here, $d_i^{\leftrightarrow }$ refers to the number of nodes to which i is connected by an edge in both directions. A is the adjacency matrix of the network G, such that ${[{(A+A^T)}^3]}_{ii}$ is the number of triangles that node i is included in, irrespective of the directions of the edges.

Clustering measures the prevalence of small, well-connected cliques in the graph, and is a measure of the local behaviour of small-world networks (Watts and Strogatz, 1998).

2.2.5 Degree variance

The variance of the node degrees quantifies the heterogeneity of a graph (Snijders, 1981). That is, whether nodes in a network are all similarly well-connected, resulting in a low degree variance, or there exist some nodes with many neighbours and some nodes with few, resulting in a high variance. For directed networks, the degree in question is taken to be the out-degree, such that the degree variance is given as

$\mathbb{V}=\frac{1}{N}{\displaystyle \sum _{i\in V}}{\left(d_i^{\mathit{out}}-d^{\bar{\mathit{o}}\mathit{ut}}\right)}^2(16)$

where $d^{\overline{\mathit{out}}}$ is the average out-degree of all nodes, equalling $\frac{\left|E\right|}{N}$, where E is the set of all directed edges in the network.

2.3 Selection of networks

A set of 10,000 20-node random binary directed networks with a mean degree of 2.5 were generated for formal analysis within this paper. This network size and mean degree is in line with typical functional connectivity networks obtained from scalp or intracranial EEG recordings (Aminoff, 2012; Sargolzaei et al., 2015; Lopes et al., 2017). Networks were generated using the NetworkX package in Python (Hagberg et al., 2008), which utilises the Erdös-Rényi algorithm (Erdös and Rényi, 1959; Gilbert, 1959) and specified as at least weakly connected. Figure 7 shows the distribution of the number of nodes contained in the FTC, n, for these networks. We observe a bi-modal distribution, where the majority of networks lie in the range of n < 5 or n > 15. An equivalent distribution for 1000 networks generated with mean degrees of 1.5 and 6 is shown for comparison. As the mean degree of generated networks increases, so does the likelihood of the FTC spanning the entire network. Equivalently as networks become more sparse, the likelihood that a large proportion of the network is contained within the FTC decreases. Our choice of a mean degree of 2.5 is a pragmatic choice, such that it contains a large sample of networks with both high and low n. BNI has been shown to depend on the network mean degree (Petkov et al., 2014). In order to avoid such influence, networks were generated with both equal mean degree and number of nodes N.

Figure 7

Figure 7. Sizes of FTCs of 10,000 randomly-generated 20-node networks of mean degree 2.5, and 1000 networks of mean degree 1.5 and 6. In the case of mean degree 2.5, for each n between 6 and 14, there were fewer than 50 networks $(<0.5\%)$ in the sample. The least-represented size of FTC was n =9, which appeared twice in the sample of 10,000. For the sample of 1000 networks with mean degree 1.5, the number of networks for each n ≥ 9 was less than 5, with no networks having n ≥ 18. For the sample of 1000 networks with mean degree 6, there were no networks with 2 ≤ n ≤ 18, and 97.5% of the networks had n = 20.

2.4 Correlations

Nonlinear correlations are used to measure the relationships between network measures and measures of the increase in BNI with increasing λ₀. We calculate Kendall correlation coefficients between AUC/QD and efficiency, mean local clustering coefficient, trophic incoherence and degree variance; both across all sizes of FTC and for individual sizes of FTC. Significance is evaluated at α = 0.0001, with a Bonferroni correction for multiple comparisons.

3 Results

3.1 3 node networks

In Figure 8 we present the BNI for all 13 non-isomorphic 3 node networks (Benjamin et al., 2012) as a function of the baseline excitability λ₀. We observe that network trajectories in λ₀ show a tendency to group together according to the size of their FTC, n. This is consistent with the work of Benjamin et al. (Benjamin et al., 2012) who have shown that network escape times for 3 node networks group according to n. This behaviour cannot be trivially approximated by counting the number of edges within each network, e.g., 3-edge networks with n = 1, 2, or 3 exhibit distinct behaviour. The set of networks with n = 1 exhibit the most gradual increase in BNI as λ₀ increases. The initial incline of networks with n = 2 is intermediate (between n = 1 and n = 3), with a sharp increase in incline occurring in the region of λ₀ = 0.85. Networks with n = 3 exhibit similar behaviour, albeit with a steeper increase in BNI. The exception to this pattern is the 3-node network consisting of a single sink driven by two otherwise isolated nodes (n = 2), indicated in Figures 5, 8.

Figure 8

Figure 8. (A) BNI as a function of increasing baseline excitability λ₀ for all non-isomorphic three-node networks. All parameters follow the standard values stated in Table 1. Networks with an FTC containing 1 node are shown in purple, those with an FTC containing 2 nodes shown in green and those with an FTC containing the whole network are shown in blue. The trajectory for a single-node sink (ref Figure 5) driven by two otherwise isolated nodes (n = 2) is shown in dark green. Errors are calculated as the standard deviation over 5 realisations of the noise co-efficient α. (B) Focus on the region of 0.8 < λ₀ < 0.9 where trajectories intersect.

The behaviour of these network trajectories was quantified using two measures - the area under the curve (AUC) and the quartile distance (QD). The distribution of these measures for n = 1, 2, 3 is shown in Figure 9. When the single-node sink is excluded, it is clear that the distribution of AUC/QD is dissimilar for different values of n.

Figure 9

Figure 9. Boxplots of the distributions of quartile distance and AUC for each size of FTC in 3-node networks. We exclude the anomalous case for n = 2 in which the two nodes in the FTC are not connected by an edge (see Figures 5, 8; this network has an AUC of 0.0904 and QD of 0.0950). We see that there are clear trends exhibited in the remaining: as n increases, AUC increases and QD decreases. Boxes are drawn from the first to third quartiles with the median value marked in blue. Whiskers are drawn at the maximum and minimum data points of the set, excluding any outliers. Outliers are marked by asterisks (*).

3.2 20 node networks

Having characterised the dynamics in small networks (N = 3), we now consider networks of size N = 20. Figure 10 presents the distribution of QD/AUC measures for increasing values of n. For both measures, values tend to remain similar for networks with FTC sizes beyond n = 6.

Figure 10

Figure 10. Boxplots of the distributions of quartile distance and AUC for each size of FTC in 20-node networks. Beyond n = 6 we see very little variation in either measure which can be associated to n. Boxes are drawn from the first to third quartiles with the median value marked in blue. Whiskers are drawn at the maximum and minimum data points of the set, excluding any outliers. Outliers are marked by stars.

As n alone is no longer sufficient to predict how network BNI will evolve as λ₀ increases, we quantified the relationship between additional network measures and AUC/QD. It should be noted that the number of networks in the range 6 ≤ n ≤ 14 was very low ($<50$ for each n out of a total 10,000, while $>2000$ are present for n = 1, 2). The networks within this range were combined with a sample of 50 networks for each remaining value of n, and evaluated by efficiency, mean clustering coefficient, trophic incoherence and degree variance.

Figure 11 displays the relationships of network measures with AUC and QD across all sizes of FTC. Significant positive correlations were observed for AUC with network efficiency (Kendall tau = 0.393, p < 0.0001) and trophic incoherence (Kendall tau = 0.357, p < 0.0001). A significant negative correlation was observed for AUC with degree variance, (Kendall tau = −0.472, p < 0.0001). No strong correlation was found between AUC and mean clustering coefficient (Kendall tau = −0.0827, p = 1.670 × 10⁻³). Significant negative correlations were observed for QD with network efficiency (Kendall tau = −0.347, p < 0.0001) and trophic incoherence (Kendall tau = −0.398, p < 0.0001). A significant positive correlation was observed for QD with degree variance, (Kendall tau = 0.567, p < 0.0001). No significant correlation was found between QD and mean clustering coefficient (Kendall tau = 0.0575, p = 0.03153).

Figure 11

Figure 11. Scatter plots detailing the relationships between AUC and network measures for 20-node networks, coloured by the size of FTC.

We will now observe in closer detail these relationships for each size of FTC. Due to the low numbers of networks generated for the intermediate values of n, further analysis is conducted only on high and low sizes of FTC: 1 ≤ n ≤ 5 and 16 ≤ n ≤ 20. For these, Kendall correlation coefficients between our network metrics, and network AUC/QD are shown in Figure 12.

Figure 12

Figure 12. (A) 2D correlation plot between AUC and network metrics for 20-node networks (efficiency (NE), mean clustering co-efficient (MCC), trophic incoherence (TI) and degree variance (DV)) for n = 1 − 5, 16 − 20 for 20 node randomly generated networks with a mean degree of 2.5. Kendall correlation coefficients are colour-coded under the right axis colorbar with exact values shown in each cell. Significance values are labelled as * where p < 0.05, ** where p < 0.01 and *** where p < 0.0001. (B) Correlation plot between QD and network metrics for 20-node networks.

Figure 12 shows that there is a clear trend for a negative correlation between the network efficiency and QD, and between trophic incoherence and QD. For high values of n, there is a weak positive correlation between degree variance and QD. No consistent trend in correlation is shown between QD or AUC and mean clustering co-efficient for high or low n. For networks with n ≤ 5 there is a trend of moderate positive correlation of AUC with trophic incoherence and network efficiency, respectively.

To illustrate the results shown above, the trajectories of highest and lowest efficiency and trophic incoherence for n = 3, 18 are shown in Figure 13. For both n = 3 and n = 18, the initial incline of the high-efficiency network is lower than that of the low-efficiency network. The maximum slope of the high-efficiency networks is greater, with trajectories intersecting at approximately λ = 0.85. Similar behaviour is shown for the networks with highest and lowest trophic incoherence. Equivalent plots for n = 1, 5, 16, 20 can be found in Supplementary Figures S1–S4 in the Supplementary Material.

Figure 13

Figure 13. Comparison of trajectories for 20-node networks with minimal and maximal values of efficiency and trophic incoherence for n = 3, 18.

It is important to point out that the relationship between network measures and AUC/QD can be influenced by the choice of network mean-degree (see Supplementary Figure S5 in the Supplementary Material). However, regarding practical applications, this is not a limiting factor of this framework as connectivity networks obtained from brain imaging can be thresholded to match the values presented in this work.

In our analysis, we see that our results for AUC are mirrored in the results for QD; and indeed that AUC and QD are significantly negatively correlated (Kendall tau = −0.581, p < 0.0001), as seen in Supplementary Figure S6 in the Supplementary Material. However, AUC and QD quantify two different aspects of the trajectories of BNI in increasing λ₀: AUC expresses the overall increase of the curve across all λ₀, while QD is an expression of the slope of the trajectory. This is demonstrated in Figure 14, which shows the trajectory of BNI in increasing λ₀ for two networks with the same QD of 0.0775 and AUCs of 0.1172 and 0.1353; and two networks with similar AUCs of 0.1269 and 0.1268 but QDs of 0.0350 and 0.0899 respectively. These are indicated on the scatter plot in Supplementary Figure S6.

Figure 14

Figure 14. Illustration of network trajectories with similar AUC (QD) and strongly different QD (AUC).

Amongst the network measures, network efficiency and trophic incoherence are strongly positively correlated (Kendall tau = 0.693, p < 0.0001) and both of these are significantly negatively correlated with degree variance (Kendall tau = −0.309, −0.328 respectively with p-values $<0.0001$). MCC has no significant correlation with network efficiency (Kendall tau = −0.0784, p = 0.00291), trophic incoherence (Kendall tau = −0.0291, p = 0.260) or degree variance (Kendall tau = 0.0905, p = 0.000738). The relationships between these measures are displayed in Supplementary Figure S7.

Size of FTC has a significant positive correlation with network efficiency (Kendall tau = 0.375, p < 0.0001) and trophic incoherence (Kendall tau = 0.302, p < 0.0001) and a significant negative correlation with degree variance (Kendall tau = −0.366, p < 0.0001). No significant correlation is found between the size of FTC and mean clustering coefficient (Kendall tau = −0.0697, p = 0.0104). We see upon observation in Supplementary Figure S8 that as with AUC and QD, where trends exist between FTC and other network measures, they are not in general monotonic.

3.3 128 node networks

We sought to evaluate whether our findings were consistent when networks of sizes comparable to larger EEG generated networks were considered. Networks of size N = 128 were chosen for analysis. The FTC distribution of 128 node networks with a mean degree of 2.5 was strongly skewed towards lower sizes of FTC. Networks were therefore sampled for analysis within a range of n = 4 − −12. Due to this limited range of FTC sizes compared to network size we do not explore the relationships of AUC and QD with size of FTC at this scale.

We see that on a group level the relationships between our chosen network measures and AUC/QD are preserved as network size increases from N = 20 to N = 128. Figure 15 provides a more detailed examination of these relationships for each value of FTC considered. At a group level, a significant positive correlation was observed for AUC with network efficiency (Kendall tau = 0.4789, p < 0.0001) and trophic incoherence (Kendall tau = 0.4792, p < 0.0001). A significant negative correlation was observed for AUC with degree variance (Kendall tau = −0.3515, p < 0.0001). No strong correlation was found between AUC and mean clustering coefficient (Kendall tau = 0.0530, p = 0.0928.) Significant negative correlations were observed for QD with network efficiency (Kendall tau = -0.3912, p < 0.0001) and trophic incoherence (Kendall tau = −0.3277, p < 0.0001). A significant positive correlation was found for QD with degree variance (Kendall tau = 0.3311, p < 0.0001). No significant correlation was found between QD and mean clustering co-efficient (Kendall tau = −0.0281, p = 0.4379). Figure 15 shows that there is a clear trend for positive correlation between AUC and efficiency, trophic incoherence for all values of FTC considered. No consistent trend is shown between QD and any network metrics within each FTC subset.

Figure 15

Figure 15. (A) 2D correlation plot between AUC and network metrics (efficiency (NE), mean clustering co-efficient (MCC), trophic incoherence (TI) and degree variance (DV)) for n = 4 − 12 for 128 node randomly generated networks with a mean degree of 2.5. Kendall correlation coefficients are colour-coded under the right axis colorbar with exact values shown in each cell. Significance values are labelled as * where p < 0.05, ** where p < 0.01 and *** where p < 0.0001. (B) Correlation plot between QD and network metrics.

4 Discussion

In this study, we propose a model framework that can be used to understand changes in treatment response over time. In this context, treatment could be considered as medication, or more invasive options such as surgery or stimulation. This framework considers the interaction of a set of network features: size of the first transitive component (FTC), efficiency, clustering, incoherence, and heterogeneity (degree variance).

Our results show that, for small networks (3 nodes), distinct patterns of BNI increase as a function of the baseline excitability (λ₀) can be observed. These patterns were shown to be well characterised by the number of nodes in the FTC. Networks with FTC of size 1 have the most gradual (less steep) increase in BNI. From the nature of these networks (seen on the first row in Figure 5), it is observed that the dynamics are mainly influenced by a source node, shown in purple in all networks. This node is not being “controlled” (i.e., forced to keep in the fixed point via diffusive coupling) by any other node, therefore the transition into seizure states is mainly dictated by the effect of the smooth increase in λ₀ on the driving node. This configuration is a conceptual representation of a scenario where seizures emerge in a specific brain region and can easily spread to the rest of the brain.

For networks with FTC of size 3, the graph is strongly connected, and all nodes are controlling each other to some extent. This explains why these networks remain with low values of BNI when 0 < λ₀⪅0.85. However, as λ₀ increases, nodes which transition to the limit-cycle encourage other nodes to do the same. The higher levels of synchronization in these networks lead to a steeper increase in BNI. In this scenario isolated abnormal activity cannot arrest other brain regions, and seizures can only emerge if those regions were already close to the seizure threshold.

For networks with FTC of size 2, levels of BNI increase are intermediary when two interconnected nodes are influencing a third node (light green in Figure 5). In this scenario, we observe a similar but less strong effect than that seen for networks with FTC of size 3 (fully connected). However, when two sources are independently influencing a third node (dark green in Figure 5), the profile of BNI increase is much slower than for any other networks. In this case, the competing sources struggle to arrest the third node at any given time, and the values of BNI are comparatively small even for large values of λ₀ (see Figure 8). The results observed here are in line with the analysis of 3- and 4-node networks presented in (Junges et al., 2020). The different behaviour observed when competing effects are present in the network suggests a more complex relationship between BNI and λ₀ for larger networks. Effectively this occurs due to sub-networks that interact in intricate ways.

We see this phenomenon clearly for larger networks (20 nodes), where the size of the FTC is not enough to fully characterize the relationship between BNI and λ₀. Figure 10 shows that for small FTC sizes (n⪅5), a negative (positive) correlation is observed between n and AUC (QD), with the lowest values of AUC (highest QD) being observed for FTC sizes of 4–5. Networks with larger FTC sizes (n > 6) tend to present increased AUC and decreased QD. A more detailed investigation shows that AUC tends to be lower and QD tends to be higher in less efficient, incoherent and more heterogeneous networks (Figure 12). This relationship with AUC is clearer in networks with small FTC. These results suggest that a network structure containing long distances of communication between nodes, a hierarchical flow and nodes with different levels of influence, in addition to a small driving subnetwork, can more effectively help control the emergence and spread of seizures. These properties may be present, for example, in a graph containing strongly-connected subnetworks embedded within a hierarchy. Interestingly, the clustering coefficient does not seem to influence a network’s response to an increase in λ₀.

The results described above suggest that optimal robustness to an increase in baseline excitability are observed for networks with FTC sizes of 4–5. Minimal AUC and maximal QD suggest that increases in the baseline excitability have a smaller and more gradual effect on seizure propensity in these networks. At the same time, a similar robustness effect is seen for networks with lower efficiency, incoherence and higher degree variance. The effect for these three measures is preserved in larger networks of size 128.

Previous works compared functional networks obtained using electroencephalography (EEG) recordings from people with epilepsy and healthy controls. A recent literature review of biomarkers of idiopathic epilepsy from resting-state EEG explored graph-based markers and showed evidence that, in specific frequency bands, networks from people with epilepsy tend to have elevated mean degree, degree variance, and path length (Faiman et al., 2021). The observed effect of epilepsy in these measures are in line with the results presented in this work. However it is important to notice that our study explores the relationship between the network markers and changes in seizure propensity, and not with seizure propensity itself. The relationship between trophic incoherence and epilepsy, to the best f our knowledge, has not been explored. Nevertheless, studies have shown that spreading processes in networks (such as seizures) can be strongly affected by trophic organisation (Klaise and Johnson, 2016).

Once validated using networks obtained from brain imaging modalities (e.g., electroencephalography) and long term post-intervention seizure monitoring, the framework presented in this work can provide prognostic markers of seizure propensity progression. This will support the development of personalized intervention strategies, aiming to achieve long term seizure freedom. In the context of surgical intervention, network-based models of seizure transition have been extensively used to predict optimal surgical strategies via the estimation of seizure propensity after intervention (Goodfellow et al., 2016; Lopes et al., 2017; Sinha et al., 2017; Lopes et al., 2020). However, these methods tend to ignore potential long-term changes in the system, which could lead to seizure recurrence. The framework presented in this study can significantly extend the current potential of these models by quantifying robustness to mechanisms associated to the honeymoon effect. This would ultimately contribute to support pre-surgical planning via the identification of strategies leading to more robust seizure control.

A limitation of this work is the relatively simple description of the mechanisms leading to the increase in seizure propensity after therapeutic interventions (honeymoon effect). Here, we are representing such an effect via a linear and spatially homogeneous increase in baseline excitability. In reality, such effect might be best represented by complex combinations of local and more intricate changes in excitability, as well as in the network structure itself. Additionally, the specific nature of these mechanisms can depend on several factors, like epilepsy type, age at intervention (neurodevelopmental stage), and/or intervention type. However, this work does not aim to propose a tool to comprehensively estimate the effects of increased seizure propensity in every circumstance. Instead, we propose a general framework to represent the honeymoon effect, which can be modified and validated to describe specific scenarios and provide valuable prognostic markers.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author contributions

EGH: Conceptualization, Formal Analysis, Visualization, Writing–original draft, Writing–review and editing. PK: Conceptualization, Formal Analysis, Visualization, Writing–original draft, Writing–review and editing. JT: Conceptualization, Writing–review and editing. WW: Conceptualization, Supervision, Writing–original draft, Writing–review and editing. LJ: Conceptualization, Supervision, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. JT acknowledges the support of the EPSRC via grants EP/T0277031/1, EP/W035030/1. WW acknowledges the support of Epilepsy Research United Kingdom via grant F2002. LJ acknowledges the support of the Waterloo Foundation via grant no. 1970/3346. All authors gratefully acknowledge the support of the University of Birmingham Dynamic Investment Fund.

Acknowledgments

A preprint of this article was made available on bioRxiv in January 2024 (Harrington et al., 2024).

Conflict of interest

WW and JRT are co-founders, directors and shareholders of Neuronostics Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnetp.2024.1308501/full#supplementary-material

References

Abou-Khalil, B., Löscher, W., and Schmidt, D. (2008). “Case 96 - why do some patients seem to develop tolerance to aeds? development of antiepileptic drug tolerance in a patient with temporal lobe epilepsy,” in Puzzling cases of epilepsy. Editors D. Schmidt, and S. C. Schachter Second Edition (San Diego: Academic Press), 420–432. doi:10.1016/B978-0-12-374005-2.00096-X

CrossRef Full Text | Google Scholar

Aminoff, M. J. (2012). Aminoff’s electrodiagnosis in clinical neurology E-book: expert consult-online and print. Elsevier Health Sciences.

Google Scholar

Bassett, D. S., and Bullmore, E. (2006). Small-world brain networks. Neurosci. 12, 512–523. doi:10.1177/1073858406293182

PubMed Abstract | CrossRef Full Text | Google Scholar

Benjamin, O., Fitzgerald, T. H. B., Ashwin, P., Tsaneva-Atanasova, K., Chowdhury, F., Richardson, M. P., et al. (2012). A phenomenological model of seizure initiation suggests network structure may explain seizure frequency in idiopathic generalised epilepsy. J. Math. Neurosci. 2, 1–30. doi:10.1186/2190-8567-2-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Boggs, J. G., Nowack, W. J., and Drinkard, C. R. (2000). Analysis of the “honeymoon effect” in adult epilepsy patients. Epilepsia 41, 222. doi:10.1111/j.1528-1157.2000.tb01727.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Brodie, M. J., Barry, S. J., Bamagous, G. A., Norrie, J. D., and Kwan, P. (2012). Patterns of treatment response in newly diagnosed epilepsy. Neurology 78, 1548–1554. doi:10.1212/WNL.0b013e3182563b19

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Z., Brodie, M. J., Liew, D., and Kwan, P. (2018). Treatment outcomes in patients with newly diagnosed epilepsy treated with established and new antiepileptic drugs: a 30-year longitudinal cohort study. JAMA Neurol. 75, 279–286. doi:10.1001/jamaneurol.2017.3949

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiarion, G., Sparacino, L., Antonacci, Y., Faes, L., and Mesin, L. (2023). Connectivity analysis in eeg data: a tutorial review of the state of the art and emerging trends. Bioengineering 10, 372. doi:10.3390/bioengineering10030372

PubMed Abstract | CrossRef Full Text | Google Scholar

[Dataset] Faiman, I., Smith, S., Hodsoll, J., Young, A. H., and Shotbolt, P. (2021). Resting-state EEG for the diagnosis of idiopathic epilepsy and psychogenic nonepileptic seizures: a systematic review. A Syst. Rev. 121, 108047. doi:10.1016/j.yebeh.2021.108047

CrossRef Full Text | Google Scholar

De Tisi, J., Bell, G. S., Peacock, J. L., McEvoy, A. W., Harkness, W. F., Sander, J. W., et al. (2011). The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: a cohort study. Lancet 378, 1388–1395. doi:10.1016/S0140-6736(11)60890-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Erdös, P., and Rényi, A. (1959). On random graphs i. Publ. Math. Debr. 6, 290–297. doi:10.5486/pmd.1959.6.3-4.12

CrossRef Full Text | Google Scholar

Fagiolo, G. (2007). Clustering in complex directed networks. Phys. Rev. E 76, 026107. doi:10.1103/PhysRevE.76.026107

PubMed Abstract | CrossRef Full Text | Google Scholar

Gerster, M., Taher, H., Škoch, A., Hlinka, J., Guye, M., Bartolomei, F., et al. (2021). Patient-specific network connectivity combined with a next generation neural mass model to test clinical hypothesis of seizure propagation. Front. Syst. Neurosci. 15, 675272. doi:10.3389/fnsys.2021.675272

PubMed Abstract | CrossRef Full Text | Google Scholar

Gilbert, E. N. (1959). Random graphs. Ann. Math. Statistics 30, 1141–1144. doi:10.1214/aoms/1177706098

CrossRef Full Text | Google Scholar

Goellner, E., Bianchin, M. M., Burneo, J. G., Parrent, A. G., and Steven, D. A. (2013). Timing of early and late seizure recurrence after temporal lobe epilepsy surgery. Epilepsia 54, 1933–1941. doi:10.1111/epi.12389

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodfellow, M., Rummel, C., Abela, E., Richardson, M. P., Schindler, K., and Terry, J. R. (2016). Estimation of brain network ictogenicity predicts outcome from epilepsy surgery. Sci. Rep. 6, 29215. doi:10.1038/srep29215

PubMed Abstract | CrossRef Full Text | Google Scholar

Goodfellow, M., Schindler, K., and Baier, G. (2012). Self-organised transients in a neural mass model of epileptogenic tissue dynamics. NeuroImage 59, 2644–2660. doi:10.1016/j.neuroimage.2011.08.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Hagberg, A., Swart, P., and S Chult, D. (2008). Exploring network structure, dynamics, and function using NetworkX. Tech. rep. Los Alamos, NM United States: Los Alamos National Lab.

Google Scholar

Harrington, G., Kissack, P., Terry, J. R., Woldman, W., and Junges, L. (2024). Treatment effects in epilepsy: a mathematical framework for understanding response over time. bioRxiv. doi:10.1101/2024.01.22.576627

CrossRef Full Text | Google Scholar

Hebbink, J., Meijer, H., Huiskamp, G., van Gils, S., and Leijten, F. (2017). Phenomenological network models: lessons for epilepsy surgery. Epilepsia 58, e147–e151. doi:10.1111/epi.13861

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V. K., Proix, T., Perdikis, D., Woodman, M. M., Wang, H., Bernard, C., et al. (2017). The Virtual Epileptic Patient: individualized whole-brain models of epilepsy spread. NeuroImage 145, 377–388. doi:10.1016/j.neuroimage.2016.04.049

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V. K., Stacey, W. C., Quilichini, P. P., Ivanov, A. I., and Bernard, C. (2014). On the nature of seizure dynamics. Brain 137, 2210–2230. doi:10.1093/brain/awu133

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiruska, P., de Curtis, M., Jefferys, J. G. R., Schevon, C. A., Schiff, S. J., and Schindler, K. (2013). Synchronization and desynchronization in epilepsy: controversies and hypotheses. J. Physiology 591, 787–797. doi:10.1113/jphysiol.2012.239590

PubMed Abstract | CrossRef Full Text | Google Scholar

Johnson, S., Domínguez-García, V., Donetti, L., and Muñoz, M. A. (2014). Trophic coherence determines food-web stability. Proc. Natl. Acad. Sci. 111, 17923–17928. doi:10.1073/pnas.1409077111

PubMed Abstract | CrossRef Full Text | Google Scholar

Junges, L., Woldman, W., Benjamin, O. J., and Terry, J. R. (2020). Epilepsy surgery: evaluating robustness using dynamic network models. Chaos Interdiscip. J. Nonlinear Sci. 30, 113106. doi:10.1063/5.0022171

PubMed Abstract | CrossRef Full Text | Google Scholar

Kalitzin, S. N., Velis, D. N., and da Silva, F. H. L. (2010). Stimulation-based anticipation and control of state transitions in the epileptic brain. Epilepsy and Behav. 17, 310–323. doi:10.1016/j.yebeh.2009.12.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Klaise, J., and Johnson, S. (2016). From neurons to epidemics: how trophic coherence affects spreading processes. Chaos 26, 065310. doi:10.1063/1.4953160

PubMed Abstract | CrossRef Full Text | Google Scholar

Kramer, M. A., Truccolo, W., Eden, U. T., Lepage, K. Q., Hochberg, L. R., Eskandar, E. N., et al. (2012). Human seizures self-terminate across spatial scales via a critical transition. Proc. Natl. Acad. Sci. 109, 21116–21121. doi:10.1073/pnas.1210047110

PubMed Abstract | CrossRef Full Text | Google Scholar

Latora, V., and Marchiori, M. (2001). Efficient behavior of small-world networks. Phys. Rev. Lett. 87, 198701. doi:10.1103/PhysRevLett.87.198701

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. A., Junges, L., Tait, L., Terry, J. R., Abela, E., Richardson, M. P., et al. (2020). Computational modelling in source space from scalp eeg to inform presurgical evaluation of epilepsy. Clin. Neurophysiol. 131, 225–234. doi:10.1016/j.clinph.2019.10.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. A., Perani, S., Yaakub, S. N., Richardson, M. P., Goodfellow, M., and Terry, J. R. (2019). Revealing epilepsy type using a computational analysis of interictal eeg. Sci. Rep. 9, 10169. doi:10.1038/s41598-019-46633-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. A., Richardson, M. P., Abela, E., Rummel, C., Schindler, K., Goodfellow, M., et al. (2017). An optimal strategy for epilepsy surgery: disruption of the rich-club? PLoS Comput. Biol. 13, e1005637. doi:10.1371/journal.pcbi.1005637

PubMed Abstract | CrossRef Full Text | Google Scholar

Lopes, M. A., Richardson, M. P., Abela, E., Rummel, C., Schindler, K., Goodfellow, M., et al. (2018). Elevated ictal brain network ictogenicity enables prediction of optimal seizure control. Front. Neurology 9, 98. doi:10.3389/fneur.2018.00098

PubMed Abstract | CrossRef Full Text | Google Scholar

Löscher, W., Rundfeldt, C., Hönack, D., and Ebert, U. (1996). Long-term studies on anticonvulsant tolerance and withdrawal characteristics of benzodiazepine receptor ligands in different seizure models in mice. i. comparison of diazepam, clonazepam, clobazam and abecarnil. J. Pharmacol. Exp. Ther. 279, 561–572.

PubMed Abstract | Google Scholar

Löscher, W., and Schmidt, D. (2006). Experimental and clinical evidence for loss of effect (tolerance) during prolonged treatment with antiepileptic drugs. Epilepsia 47, 1253–1284. doi:10.1111/j.1528-1167.2006.00607.x

PubMed Abstract | CrossRef Full Text | Google Scholar

MacKay, R. S., Johnson, S., and Sansom, B. (2020). How directed is a directed network? R. Soc. open Sci. 7, 201138. doi:10.1098/rsos.201138

PubMed Abstract | CrossRef Full Text | Google Scholar

Meisel, C., Schulze-Bonhage, A., Freestone, D., Cook, M. J., Achermann, P., and Plenz, D. (2015). Intrinsic excitability measures track antiepileptic drug action and uncover increasing/decreasing excitability over the wake/sleep cycle. Proc. Natl. Acad. Sci. 112, 14694–14699. doi:10.1073/pnas.1513716112

PubMed Abstract | CrossRef Full Text | Google Scholar

Olmi, S., Petkoski, S., Guye, M., Bartolomei, F., and Jirsa, V. (2019). Controlling seizure propagation in large-scale brain networks. PLOS Comput. Biol. 15, e1006805–e1006823. doi:10.1371/journal.pcbi.1006805

PubMed Abstract | CrossRef Full Text | Google Scholar

Petkov, G., Goodfellow, M., Richardson, M. P., and Terry, J. R. (2014). A critical role for network structure in seizure onset: a computational modeling approach. Front. neurology 5, 261. doi:10.3389/fneur.2014.00261

PubMed Abstract | CrossRef Full Text | Google Scholar

Petrik, S., San Antonio-Arce, V., Steinhoff, B. J., Syrbe, S., Bast, T., Scheiwe, C., et al. (2021). Epilepsy surgery: late seizure recurrence after initial complete seizure freedom. Epilepsia 62, 1092–1104. doi:10.1111/epi.16893

PubMed Abstract | CrossRef Full Text | Google Scholar

Proix, T., Bartolomei, F., Guye, M., and Jirsa, V. K. (2017). Individual brain structure and modelling predict seizure propagation. Brain 140, 641–654. doi:10.1093/brain/awx004

PubMed Abstract | CrossRef Full Text | Google Scholar

Riss, J., Cloyd, J., Gates, J., and Collins, S. (2008). Benzodiazepines in epilepsy: pharmacology and pharmacokinetics. Acta Neurol. Scand. 118, 69–86. doi:10.1111/j.1600-0404.2008.01004.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52, 1059–1069. doi:10.1016/j.neuroimage.2009.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sargolzaei, S., Cabrerizo, M., Goryawala, M., Eddin, A. S., and Adjouadi, M. (2015). Scalp eeg brain functional connectivity networks in pediatric epilepsy. Comput. Biol. Med. 56, 158–166. doi:10.1016/j.compbiomed.2014.10.018

PubMed Abstract | CrossRef Full Text | Google Scholar

Sinha, N., Dauwels, J., Kaiser, M., Cash, S. S., Brandon Westover, M., Wang, Y., et al. (2017). Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. Brain 140, 319–332. doi:10.1093/brain/aww299

PubMed Abstract | CrossRef Full Text | Google Scholar

Snijders, T. A. (1981). The degree variance: an index of graph heterogeneity. Soc. Netw. 3, 163–174. doi:10.1016/0378-8733(81)90014-9

CrossRef Full Text | Google Scholar

Téllez-Zenteno, J. F., Dhar, R., and Wiebe, S. (2005). Long-term seizure outcomes following epilepsy surgery: a systematic review and meta-analysis. Brain 128, 1188–1198. doi:10.1093/brain/awh449

PubMed Abstract | CrossRef Full Text | Google Scholar

Terry, J. R., Benjamin, O., and Richardson, M. P. (2012). Seizure generation: the role of nodes and networks. Epilepsia 53, e166–e169. doi:10.1111/j.1528-1167.2012.03560.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Touboul, J., Wendling, F., Chauvel, P., and Faugeras, O. (2011). Neural mass activity, bifurcations, and epilepsy. Neural Comput. 23, 3232–3286. doi:10.1162/NECO_a_00206

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H. E., Bénar, C. G., Quilichini, P. P., Friston, K. J., Jirsa, V. K., and Bernard, C. (2014). A systematic framework for functional connectivity measures. Front. Neurosci. 8, 405. doi:10.3389/fnins.2014.00405

PubMed Abstract | CrossRef Full Text | Google Scholar

Watts, D. J., and Strogatz, S. H. (1998). Collective dynamics of ‘small-world’ networks. Nature 393, 440–442. doi:10.1038/30918

PubMed Abstract | CrossRef Full Text | Google Scholar

Woldman, W., Cook, M. J., and Terry, J. R. (2019). Evolving dynamic networks: an underlying mechanism of drug resistance in epilepsy? Epilepsy and Behav. 94, 264–268. doi:10.1016/j.yebeh.2019.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar