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

Abstract

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

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:

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, ] at each time-step. Typical values for these model parameters can be found within Table 1.

This system is deterministic when α = 0. Within the physically permissible region of λ ∈ [0, 1] and , the dynamics of z are characterised by three steady-state solutions. The first is a stable limit cycle at . 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 . These steady state solutions are shown in Figure 2.

FIGURE 2

Eq. 2 describes the λ dynamics, and has a single stable steady-state solution in which , 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 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

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,

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 asT 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 thatA node is considered as being in a seizure-like state if . 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

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”,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

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

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

and the total weight of the node is

where and 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

where Λ is the graph Laplacian operator

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

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

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):

Here, 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 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

where is the average out-degree of all nodes, equalling , 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

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

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

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

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 ( for each n out of a total 10,000, while 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

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

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

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

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.

References

• 1

  Abou-KhalilB.LöscherW.SchmidtD. (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 SchmidtD.SchachterS. C.Second Edition (San Diego: Academic Press), 420–432. 10.1016/B978-0-12-374005-2.00096-X

  • CrossRef
  • Google Scholar

• 2

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

  • Google Scholar

• 3

  BassettD. S.BullmoreE. (2006). Small-world brain networks. Neurosci.12, 512–523. 10.1177/1073858406293182

  • CrossRef
  • Google Scholar

• 4

  BenjaminO.FitzgeraldT. H. B.AshwinP.Tsaneva-AtanasovaK.ChowdhuryF.RichardsonM. 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. 10.1186/2190-8567-2-1

  • CrossRef
  • Google Scholar

• 5

  BoggsJ. G.NowackW. J.DrinkardC. R. (2000). Analysis of the “honeymoon effect” in adult epilepsy patients. Epilepsia41, 222. 10.1111/j.1528-1157.2000.tb01727.x

  • CrossRef
  • Google Scholar

• 6

  BrodieM. J.BarryS. J.BamagousG. A.NorrieJ. D.KwanP. (2012). Patterns of treatment response in newly diagnosed epilepsy. Neurology78, 1548–1554. 10.1212/WNL.0b013e3182563b19

  • CrossRef
  • Google Scholar

• 7

  ChenZ.BrodieM. J.LiewD.KwanP. (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. 10.1001/jamaneurol.2017.3949

  • CrossRef
  • Google Scholar

• 8

  ChiarionG.SparacinoL.AntonacciY.FaesL.MesinL. (2023). Connectivity analysis in eeg data: a tutorial review of the state of the art and emerging trends. Bioengineering10, 372. 10.3390/bioengineering10030372

  • CrossRef
  • Google Scholar

• 9

  [Dataset] FaimanI.SmithS.HodsollJ.YoungA. H.ShotboltP. (2021). Resting-state EEG for the diagnosis of idiopathic epilepsy and psychogenic nonepileptic seizures: a systematic review. A Syst. Rev.121, 108047. 10.1016/j.yebeh.2021.108047

  • CrossRef
  • Google Scholar

• 10

  De TisiJ.BellG. S.PeacockJ. L.McEvoyA. W.HarknessW. F.SanderJ. W.et al (2011). The long-term outcome of adult epilepsy surgery, patterns of seizure remission, and relapse: a cohort study. Lancet378, 1388–1395. 10.1016/S0140-6736(11)60890-8

  • CrossRef
  • Google Scholar

• 11

  ErdösP.RényiA. (1959). On random graphs i. Publ. Math. Debr.6, 290–297. 10.5486/pmd.1959.6.3-4.12

  • CrossRef
  • Google Scholar

• 12

  FagioloG. (2007). Clustering in complex directed networks. Phys. Rev. E76, 026107. 10.1103/PhysRevE.76.026107

  • CrossRef
  • Google Scholar

• 13

  GersterM.TaherH.ŠkochA.HlinkaJ.GuyeM.BartolomeiF.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. 10.3389/fnsys.2021.675272

  • CrossRef
  • Google Scholar

• 14

  GilbertE. N. (1959). Random graphs. Ann. Math. Statistics30, 1141–1144. 10.1214/aoms/1177706098

  • CrossRef
  • Google Scholar

• 15

  GoellnerE.BianchinM. M.BurneoJ. G.ParrentA. G.StevenD. A. (2013). Timing of early and late seizure recurrence after temporal lobe epilepsy surgery. Epilepsia54, 1933–1941. 10.1111/epi.12389

  • CrossRef
  • Google Scholar

• 16

  GoodfellowM.RummelC.AbelaE.RichardsonM. P.SchindlerK.TerryJ. R. (2016). Estimation of brain network ictogenicity predicts outcome from epilepsy surgery. Sci. Rep.6, 29215. 10.1038/srep29215

  • CrossRef
  • Google Scholar

• 17

  GoodfellowM.SchindlerK.BaierG. (2012). Self-organised transients in a neural mass model of epileptogenic tissue dynamics. NeuroImage59, 2644–2660. 10.1016/j.neuroimage.2011.08.060

  • CrossRef
  • Google Scholar

• 18

  HagbergA.SwartP.S ChultD. (2008). Exploring network structure, dynamics, and function using NetworkX. Tech. rep. Los Alamos, NM United States: Los Alamos National Lab.

  • Google Scholar

• 19

  HarringtonG.KissackP.TerryJ. R.WoldmanW.JungesL. (2024). Treatment effects in epilepsy: a mathematical framework for understanding response over time. bioRxiv. 10.1101/2024.01.22.576627

  • CrossRef
  • Google Scholar

• 20

  HebbinkJ.MeijerH.HuiskampG.van GilsS.LeijtenF. (2017). Phenomenological network models: lessons for epilepsy surgery. Epilepsia58, e147–e151. 10.1111/epi.13861

  • CrossRef
  • Google Scholar

• 21

  JirsaV. K.ProixT.PerdikisD.WoodmanM. M.WangH.BernardC.et al (2017). The Virtual Epileptic Patient: individualized whole-brain models of epilepsy spread. NeuroImage145, 377–388. 10.1016/j.neuroimage.2016.04.049

  • CrossRef
  • Google Scholar

• 22

  JirsaV. K.StaceyW. C.QuilichiniP. P.IvanovA. I.BernardC. (2014). On the nature of seizure dynamics. Brain137, 2210–2230. 10.1093/brain/awu133

  • CrossRef
  • Google Scholar

• 23

  JiruskaP.de CurtisM.JefferysJ. G. R.SchevonC. A.SchiffS. J.SchindlerK. (2013). Synchronization and desynchronization in epilepsy: controversies and hypotheses. J. Physiology591, 787–797. 10.1113/jphysiol.2012.239590

  • CrossRef
  • Google Scholar

• 24

  JohnsonS.Domínguez-GarcíaV.DonettiL.MuñozM. A. (2014). Trophic coherence determines food-web stability. Proc. Natl. Acad. Sci.111, 17923–17928. 10.1073/pnas.1409077111

  • CrossRef
  • Google Scholar

• 25

  JungesL.WoldmanW.BenjaminO. J.TerryJ. R. (2020). Epilepsy surgery: evaluating robustness using dynamic network models. Chaos Interdiscip. J. Nonlinear Sci.30, 113106. 10.1063/5.0022171

  • CrossRef
  • Google Scholar

• 26

  KalitzinS. N.VelisD. N.da SilvaF. H. L. (2010). Stimulation-based anticipation and control of state transitions in the epileptic brain. Epilepsy and Behav.17, 310–323. 10.1016/j.yebeh.2009.12.023

  • CrossRef
  • Google Scholar

• 27

  KlaiseJ.JohnsonS. (2016). From neurons to epidemics: how trophic coherence affects spreading processes. Chaos26, 065310. 10.1063/1.4953160

  • CrossRef
  • Google Scholar

• 28

  KramerM. A.TruccoloW.EdenU. T.LepageK. Q.HochbergL. R.EskandarE. N.et al (2012). Human seizures self-terminate across spatial scales via a critical transition. Proc. Natl. Acad. Sci.109, 21116–21121. 10.1073/pnas.1210047110

  • CrossRef
  • Google Scholar

• 29

  LatoraV.MarchioriM. (2001). Efficient behavior of small-world networks. Phys. Rev. Lett.87, 198701. 10.1103/PhysRevLett.87.198701

  • CrossRef
  • Google Scholar

• 30

  LopesM. A.JungesL.TaitL.TerryJ. R.AbelaE.RichardsonM. P.et al (2020). Computational modelling in source space from scalp eeg to inform presurgical evaluation of epilepsy. Clin. Neurophysiol.131, 225–234. 10.1016/j.clinph.2019.10.027

  • CrossRef
  • Google Scholar

• 31

  LopesM. A.PeraniS.YaakubS. N.RichardsonM. P.GoodfellowM.TerryJ. R. (2019). Revealing epilepsy type using a computational analysis of interictal eeg. Sci. Rep.9, 10169. 10.1038/s41598-019-46633-7

  • CrossRef
  • Google Scholar

• 32

  LopesM. A.RichardsonM. P.AbelaE.RummelC.SchindlerK.GoodfellowM.et al (2017). An optimal strategy for epilepsy surgery: disruption of the rich-club?PLoS Comput. Biol.13, e1005637. 10.1371/journal.pcbi.1005637

  • CrossRef
  • Google Scholar

• 33

  LopesM. A.RichardsonM. P.AbelaE.RummelC.SchindlerK.GoodfellowM.et al (2018). Elevated ictal brain network ictogenicity enables prediction of optimal seizure control. Front. Neurology9, 98. 10.3389/fneur.2018.00098

  • CrossRef
  • Google Scholar

• 34

  LöscherW.RundfeldtC.HönackD.EbertU. (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.

  • Google Scholar

• 35

  LöscherW.SchmidtD. (2006). Experimental and clinical evidence for loss of effect (tolerance) during prolonged treatment with antiepileptic drugs. Epilepsia47, 1253–1284. 10.1111/j.1528-1167.2006.00607.x

  • CrossRef
  • Google Scholar

• 36

  MacKayR. S.JohnsonS.SansomB. (2020). How directed is a directed network?R. Soc. open Sci.7, 201138. 10.1098/rsos.201138

  • CrossRef
  • Google Scholar

• 37

  MeiselC.Schulze-BonhageA.FreestoneD.CookM. J.AchermannP.PlenzD. (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. 10.1073/pnas.1513716112

  • CrossRef
  • Google Scholar

• 38

  OlmiS.PetkoskiS.GuyeM.BartolomeiF.JirsaV. (2019). Controlling seizure propagation in large-scale brain networks. PLOS Comput. Biol.15, e1006805–e1006823. 10.1371/journal.pcbi.1006805

  • CrossRef
  • Google Scholar

• 39

  PetkovG.GoodfellowM.RichardsonM. P.TerryJ. R. (2014). A critical role for network structure in seizure onset: a computational modeling approach. Front. neurology5, 261. 10.3389/fneur.2014.00261

  • CrossRef
  • Google Scholar

• 40

  PetrikS.San Antonio-ArceV.SteinhoffB. J.SyrbeS.BastT.ScheiweC.et al (2021). Epilepsy surgery: late seizure recurrence after initial complete seizure freedom. Epilepsia62, 1092–1104. 10.1111/epi.16893

  • CrossRef
  • Google Scholar

• 41

  ProixT.BartolomeiF.GuyeM.JirsaV. K. (2017). Individual brain structure and modelling predict seizure propagation. Brain140, 641–654. 10.1093/brain/awx004

  • CrossRef
  • Google Scholar

• 42

  RissJ.CloydJ.GatesJ.CollinsS. (2008). Benzodiazepines in epilepsy: pharmacology and pharmacokinetics. Acta Neurol. Scand.118, 69–86. 10.1111/j.1600-0404.2008.01004.x

  • CrossRef
  • Google Scholar

• 43

  RubinovM.SpornsO. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage52, 1059–1069. 10.1016/j.neuroimage.2009.10.003

  • CrossRef
  • Google Scholar

• 44

  SargolzaeiS.CabrerizoM.GoryawalaM.EddinA. S.AdjouadiM. (2015). Scalp eeg brain functional connectivity networks in pediatric epilepsy. Comput. Biol. Med.56, 158–166. 10.1016/j.compbiomed.2014.10.018

  • CrossRef
  • Google Scholar

• 45

  SinhaN.DauwelsJ.KaiserM.CashS. S.Brandon WestoverM.WangY.et al (2017). Predicting neurosurgical outcomes in focal epilepsy patients using computational modelling. Brain140, 319–332. 10.1093/brain/aww299

  • CrossRef
  • Google Scholar

• 46

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

  • CrossRef
  • Google Scholar

• 47

  Téllez-ZentenoJ. F.DharR.WiebeS. (2005). Long-term seizure outcomes following epilepsy surgery: a systematic review and meta-analysis. Brain128, 1188–1198. 10.1093/brain/awh449

  • CrossRef
  • Google Scholar

• 48

  TerryJ. R.BenjaminO.RichardsonM. P. (2012). Seizure generation: the role of nodes and networks. Epilepsia53, e166–e169. 10.1111/j.1528-1167.2012.03560.x

  • CrossRef
  • Google Scholar

• 49

  TouboulJ.WendlingF.ChauvelP.FaugerasO. (2011). Neural mass activity, bifurcations, and epilepsy. Neural Comput.23, 3232–3286. 10.1162/NECO_a_00206

  • CrossRef
  • Google Scholar

• 50

  WangH. E.BénarC. G.QuilichiniP. P.FristonK. J.JirsaV. K.BernardC. (2014). A systematic framework for functional connectivity measures. Front. Neurosci.8, 405. 10.3389/fnins.2014.00405

  • CrossRef
  • Google Scholar

• 51

  WattsD. J.StrogatzS. H. (1998). Collective dynamics of ‘small-world’ networks. Nature393, 440–442. 10.1038/30918

  • CrossRef
  • Google Scholar

• 52

  WoldmanW.CookM. J.TerryJ. R. (2019). Evolving dynamic networks: an underlying mechanism of drug resistance in epilepsy?Epilepsy and Behav.94, 264–268. 10.1016/j.yebeh.2019.03.003

  • CrossRef
  • Google Scholar