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.


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 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 pretreatment 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.
success of treatment with a later return of seizures suggest some mechanism by which the brain adapts towards increased seizure risk (Scenario C).
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.

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 λ 0 the constant baseline level of excitability.τ is a real constant modulating the rate at which a node transitions from the seizurelike 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, dt √ ] 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 |z| ≥ 0, the dynamics of z are characterised by three steady-state solutions.The first is a stable limit cycle at |z| 2 1 + λ √ .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 |z| 2 1 − λ √ .These steady state solutions are shown in Figure 2.
Eq. 2 describes the λ dynamics, and has a single stable steadystate solution in which λ λ 0 − |z| 2 , illustrated in Figure 2. When z is close to its fixed point, in the inter-ictal state, |z| 2 is very small and λ therefore tends to λ 0 ; in the case where |z| is large, λ tends towards values less than λ 0 , 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 λ = λ 0 .As λ 0 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, λ = λ 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 interictal state.
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.

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 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 A node is considered as being in a seizure-like state if |z| 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, λ 0 .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 λ 0 , 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 λ 0 .We consider both extremes of network connectivity -a fullyconnected 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 < 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.
We quantify the increase of network BNI under increasing λ 0 by two metrics.The area under the curve (AUC) quantifies the overall increase in seizure propensity upon increase of λ 0 .We define the "quartile distance", such that it is the increase in λ 0 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 λ 0 .

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.Trajectories of BNI with increasing λ 0 for the fully-connected and completely disconnected 20-node networks.

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 stronglyconnected (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 1 , n 2 , . . ., 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.

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.
In the uniformly-weighted case, the imbalance of a node is and the total weight of the node is where d in i and d out i 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

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

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, d ↔ i 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 smallworld networks (Watts and Strogatz, 1998).

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 d out is the average out-degree of all nodes, equalling |E| N , where E is the set of all directed edges in the network.

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

Correlations
Nonlinear correlations are used to measure the relationships between network measures and measures of the increase in BNI with increasing λ 0 .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.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 leastrepresented 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.

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 λ 0 .We observe that network trajectories in λ 0 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 λ 0 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 = 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.
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.

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.
As n alone is no longer sufficient to predict how network BNI will evolve as λ 0 increases, we quantified the relationship between (A) BNI as a function of increasing baseline excitability λ 0 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 < 0.9 where trajectories intersect.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 (*).
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.
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 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

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 Comparison of trajectories for 20-node networks with minimal and maximal values of efficiency and trophic incoherence for n = 3, 18. 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.

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 (λ 0 ) 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 λ 0 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 )0.85.However, as λ 0 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 λ 0 (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 λ 0 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 λ 0 .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 λ 0 .
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 postintervention 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.

FIGURE 1
FIGURE 1 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 λ 0 , which is parabolic with a peak at (|z| 0, λ λ 0 ).Grey dashed lines indicate the maximum and minimum values of λ 0 .When |z| 0 in the interictal state, λ converges to λ 0 .In the ictal state, |z| takes on larger values and consequently λ converges to lower values.

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 anticlockwise.(B) Simulated single-node EEG activity Re(z) over a timescale of T = 200s.(C) Variation in the slow excitability variable λ.

FIGURE 14
FIGURE 14Illustration of network trajectories with similar AUC (QD) and strongly different QD (AUC).
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.