A mathematical model of neuroimmune interactions in epileptogenesis for discovering treatment strategies.

The development of epilepsy (epileptogenesis) involves a complex interplay of neuronal and immune processes. Here, we present a first-of-its-kind mathematical model to better understand the relationships among these processes. Our model describes the interaction between neuroinflammation, blood-brain barrier disruption, neuronal loss, circuit remodeling, and seizures. Formulated as a system of nonlinear differential equations, the model reproduces the available data from three animal models. The model successfully describes characteristic features of epileptogenesis such as its paradoxically long timescales (up to decades) despite short and transient injuries or the existence of qualitatively different outcomes for varying injury intensity. In line with the concept of degeneracy, our simulations reveal multiple routes toward epilepsy with neuronal loss as a sufficient but non-necessary component. Finally, we show that our model allows for in silico predictions of therapeutic strategies, revealing injury-specific therapeutic targets and optimal time windows for intervention.

Introduction

Epilepsy is a common neurological disorder that affects numerous physiological mechanisms in the CNS (Lytton, 2008; Devinsky et al., 2018). Several processes play an important role in the development of epilepsy. These include the activation of innate and adaptive immune responses (Ravizza et al., 2008; Bauer et al., 2017), disruption of the blood-brain barrier (BBB) integrity (Marchi et al., 2012; Löscher and Friedman, 2020), neuronal loss (Tasch et al., 1999; Dingledine et al., 2014), and remodeling of neural circuits (Liao et al., 2011; Bertram, 2013; Jo et al., 2019). Among the common causes of acquired epilepsy are traumatic brain injury (Annegers and Coan, 2000; Pitkänen and Immonen, 2014), stroke (Olsen, 2001; Pitkänen et al., 2016), central neural system infection (Van Baalen et al., 2010; Ramantani and Holthausen, 2017), and new-onset status epilepticus (SE) (Holtkamp et al., 2005; Gaspard et al., 2018). All four lead to neuroinflammation, which has been widely studied in the context of epilepsy (Barker-Haliski et al., 2017; Rana and Musto, 2018; Vezzani et al., 2019). In particular, signaling pathways of neuroinflammation have been shown to modulate neuronal excitability (Devinsky et al., 2013; Vezzani and Viviani, 2015), seizure threshold (Heida et al., 2004; Galic et al., 2008), and severity of the seizure burden (Auvin et al., 2010a, 2010b; Tan et al., 2015; Rana and Musto, 2018). For example, Patel et al. (2017) have shown that tumor necrosis factor (TNF) knockout mice have a significantly reduced seizure burden. A similar effect was observed in animals with a deletion of TNF receptor 1 (TNFR1), while the ablation of TNF receptor 2 (TNFR2) led to an increase in seizure burden. Moreover, it has been hypothesized that noxious stimuli causing neuroinflammation may not only be associated with an acute injury but also with abnormal neuronal activity, resulting in so-called neurogenic neuroinflammation (Xanthos and Sandkühler, 2014).

The complexity of the involved processes interacting on various timescales makes understanding epileptogenesis (EPG) a great challenge. In such a complex system with several nonlinear interacting processes, mathematical modeling is a useful tool for a better understanding of the system’s dynamics. In addition, modeling helps to systematize and explain the great body of observations obtained in clinical and animal model studies, and to provide valuable predictions for further experimental studies. The mathematical modeling approach has already been successfully applied to studying the dynamics of ictogenesis: initiation, spreading, and termination of epileptic seizures (Fröhlich et al., 2008; Jirsa et al., 2014, 2017; Proix et al., 2017; Gonzalez et al., 2018; González et al., 2019). However, no computational model of EPG that implements multiple major pathomechanisms at their relevant time scales (including neuroinflammation) has yet been developed.

Therefore, in this article, we present a first-of-its-kind mathematical model that describes and simulates the course of acquired epilepsy development and accounts for the epileptogenic role of neuroimmune interactions. Reproducing the data from three animal models, our unified computational framework allows for the simulation of EPG caused by the most common types of neural injuries using a single parameter set. It suggests an explanation for, among other things, how different causes (injury types) lead to similar outcomes (epilepsies) and why the development of epilepsy may take months or years after the triggering injury in certain conditions. Our mathematical analysis reveals that these long time scales result from the underlying dynamical system slowing down and lingering in the vicinity of an unstable fixed point. In addition to the reproduction and explanation of injury-specific features of EPG, our model provides insights into the general principles of EPG and allows us to generate testable predictions for different therapeutic strategies.

Results

Model design

The model describes interactions between neuroinflammation (I), BBB disruption (B), neuronal loss (D), circuit remodeling (R), and epileptic seizures (S) upon neurological injury (Figure 1A). The probability of spontaneous recurrent seizure occurrence is assumed to depend on two seizure-promoting factors: intensity of neuroinflammation and degree of pathological circuit remodeling (Figure 1A, arrows I→S and R→S). Pathological remodeling in circuits (R) may be constituted by, among others, loss of inhibitory neurons (Sloviter, 1987; Knopp et al., 2008), abnormal excitatory synaptogenesis (Weissberg et al., 2015; Kim et al., 2017a, 2017b), or increase of recurrency in neural circuits owing to mossy fiber sprouting (Tauck and Nadler, 1985; Buckmaster, 2014). Such pathological remodeling leads to an increase of excitability in neuronal circuits and rising in the chance of critical synchronization that results in the occurrence of spontaneous seizures (arrow R→S). Facilitation of a neuroinflammatory response was shown to modulate neuronal excitability (Devinsky et al., 2013; Vezzani and Viviani, 2015) and, consequently, lower seizure threshold (Heida et al., 2004; Galic et al., 2008) (arrow I→S). In agreement with a recent finding on the inhibition of neuronal activity by microglia (Badimon et al., 2020), we assume that in conditions of profound neuroinflammation this physiological function may be impaired owing to the adoption of proinflammatory phenotypes by microglia.

Figure 1

Model overview

(A) Interactions between variables in our model.

(B) State-space plot for the rate model with three steady states: “healthy” stable steady state (black), unstable fixed point (white), and “epileptic” stable steady state corresponding to state of progressed EPG (gray). This plot illustrates the dynamical landscape of the system, which consists of two attracting points (stable steady states). The dashed black line going through the unstable fixed point is a separatrix, which separates the basin of attraction of the “healthy” steady state (shaded area) from the basin of attraction of the “epileptic” steady state. Color of arrows indicates the velocity of state changes in corresponding areas of the state-space. Red dashed line corresponds to the glial neurotoxicity threshold Θ under conditions of timescales separation I ≈ B. Model parameters are specified in supplemental item S2.

Epileptic seizures have been shown to cause metabolic stress on the CNS owing to increased energy demands associated with excessive neural activity (Zhang et al., 2015; Prager et al., 2019). Data from human and animal models also suggest that seizures induce leakiness of the BBB and that the leakiness is negatively correlated with time since the last seizure (Van Vliet et al., 2007; Rüber et al., 2018). In our model, we account for this effect of seizures (Figure 1A, arrow S→B), as well as neuroinflammation that can be caused by exposure of the parenchyma to cells or soluble factors infiltrating through the leaky BBB (Farrell et al., 2017; Löscher and Friedman, 2020) (arrow B→I). Neuroinflammation itself can cause leakiness of the BBB via the activation of cells of the neurovascular unit (Obermeier et al., 2013), which results in a positive feedback loop in our model (Figure 1A, arrow I→B). Excessive neuroinflammation may lead to neurotoxicity (Block et al., 2007; Biber et al., 2014), which is also accounted for in our model (arrow I→D). It was shown that severe and repetitive seizures kill neurons (Dingledine et al., 2014), while the same can not be said surely regarding short and isolated seizures. The model design is conservative in this regard: no direct link between seizures and neuronal death is assumed (Figure 1A). Neural loss, in turn, leads to the remodeling of neural circuits. Some remodeling may aim at maintaining functional properties. Here we only consider the kind of remodeling that leads to increased seizure susceptibility (arrow D→R). Furthermore, we also account for pathological circuit remodeling that is independent of neural loss. Examples of such remodeling are the development of chronic inhibition deficits (Kim et al., 2017a, 2017b) and excessive excitatory synaptogenesis (Weissberg et al., 2015) owing to albumin extravasation on BBB dysfunction (Figure 1A, arrow B→R).

The mathematical description of model design, implementation, and analysis are available in the STAR Methods. The rationale behind the parameter choice, the impact of changing the parameter values and specificity of results to chosen parameter set are discussed in detail in STAR Methods. The dynamical landscape of the system describing its stability is illustrated in Figure 1B.

Using a single set of parameters, the mathematical model allows for the simulation of EPG progression caused by various types of neurological injuries. Modeling of different types of injuries is carried out by the application of time sequences of perturbations mimicking pathological effects of the respective injury (Figures 2A–2C). We present simulations of three injury types and associated animal models: Theiler’s murine encephalomyelitis virus (TMEV) mouse model; chemically induced (pilocarpine) SE rodent model; BBB disruption rodent model. These models are distinct in methods of EPG induction and time course of disease pathology. The rationale behind simulation approaches for each simulated animal model is available in the STAR Methods. The simulation protocols are available in Table S2.

Figure 2

Simulation schematics for three animal models of epileptogenesis

(A) TMEV mouse model.

(B) Pilocarpine SE rodent model.

(C) BBB disruption rodent model.

Dependence of epileptogenesis on injury intensity and emergence of long timescales

The risk of epilepsy development and severity of seizure burden have been shown to depend on the intensity of the neurological injury. Patients suffering from mild traumatic brain injury (TBI) have a 2.1% cumulative probability of seizure development over 30 years, while for severe TBI it rises up to 16.7% (Annegers et al., 1998). The dose-dependence is not specific to TBI and also translates to animal models of epilepsy. In infection-induced epilepsy in mice, the increasing viral dose from 3×103 plaque-forming units (PFU) to 3 × 106 PFU leads to an increase in the fraction of animals with seizures from 25% to 80% (Libbey et al., 2011). In SE models, the duration of the acute epileptiform activity determines the incidence of epilepsy in animals (Löscher, 2015).

We tested whether our model captures the observed spectrum of dose-dependence effects of the risk of EPG (and its characteristic features) on the intensity of neurological injury. Figure 3 illustrates the results of the simulation of EPG induced by BBB disruption. The latent period duration (time period between injury and arrival of first seizure, Figure 3A) and seizure burden (Figure 3B) in simulations are in agreement with data reported in the animal model study (Weissberg et al., 2015). Accordingly, seizure manifestation was observed in our simulations either during the infusion of albumin (7 days after pump implantation), or after albumin pump removal (e.g. simulations #6, #7, #8 in Figure 3C). Similar to the animal model data, no neuronal death was observed within the observation period because neuroinflammation did not reach the neurotoxicity threshold for neuronal death (Figure 3D).

Figure 3

Dose-dependence of EPG on the intensity of neurological injury

(A) Comparison of the latent period duration in BBB disruption animal model data from Weissberg et al. (2015) and simulations with the intensity of the injury: matched to Weissberg et al. (2015) (B = 0.25, Toff = 7 days); decreased via lowering albumin concentration (B×50%); decreased via shortening the time window of albumin infusion (Toff×50%); increased via prolongation of the time window of albumin infusion (Toff×150%). Data are represented as mean ± SEM. The unpaired two-group Mann-Whitney U test was performed: ∗p < 0.05; ∗∗p < 0.01; ∗∗∗p < 0.001; ∗∗∗∗p < 0.0001.

(B) Comparison of seizure burden on the first month after injury onset in BBB disruption animal model data from Weissberg et al. (2015) and simulations (annotation identical to caption in A.).

(C) Time sequences of seizure occurrence in individual simulations. Orange bar corresponds to the time window of injury induction.

(D) Time course of neuroinflammation in individual simulations (N = 30). Red dashed line corresponds to the neurotoxicity threshold Θ. Light orange area corresponds to the time window of injury induction.

(E) EPG in response to injuries of four different intensities illustrated with seizure rate development over time post-pump implantation. Simulation results obtained with the rate model. The injury intensity control was implemented by the modification of the duration of the time window of albumin infusion (Toff).

(F) EPG in response to injuries of four different intensities (annotation identical to caption E.) illustrated over the state-space plot. The plot was obtained by performing the timescale separation procedure (see STAR Methods). The dynamical landscape of the system consists of three steady states: “healthy” (black), unstable (white) and “epileptic” (gray). Dashed black line (separatrix) separates basins of attraction of two stable steady states. Red dashed line corresponds to neurotoxicity threshold Θ. Distance between circle markers on EPG traces corresponds to time intervals of 30 days. The state of the system is being pushed away from the origin by the effect of the injury. The exact trajectory of state changes depends on the particular point, from which it takes off after the injury offset. The take-off point may be located in either of two basins of attraction, which would determine the progression or recovery after the neurological injury. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

Our model predicts that with a 50% decrease in the injury intensity the latent period duration has prolonged (5.57 ± 0.34 days vs 7.23 ± 0.47 days, Figure 3A), and the seizure burden has dropped significantly (1.24 ± 0.07 seizures/day vs 0.62 ± 0.04 seizures/day, Figure 3B). The manipulation of injury intensity used in this simulation, in the experimental setting, corresponds to lowering of the albumin concentration in the infused artificial cerebrospinal fluid. An alternative approach to the manipulation of injury intensity is prolongation or shortening of the time window of artificial cerebrospinal fluid infusion. Shortening by 50% did not have a significant effect on latent period duration (Figure 3A), but led to a significant drop in the seizure burden (1.24 ± 0.07 seizures/day vs 0.58 ± 0.03 , Figure 3B). On the other hand, the increase of injury intensity simulated by 50% prolongation of infusion time window led to a significant rise in seizure burden (1.24 ± 0.07 vs 2.13 ± 0.17 , Figure 3B).

In addition to the latent period duration and seizure burden, the intensity of the injury has been also shown to affect the risk of epilepsy development itself. The injury of low intensity does not lead to progressive EPG owing to the restoration of the healthy state (Figures 3E and 3F). Mathematically, the dose-dependence originates from the fact that injuries of low intensity fail to push the system state across the separatrix into the basin of attraction of the "epileptic” steady-state (Figure 3F). Instead, the model recovers without progressive EPG. From the biological and clinical perspectives, this would correspond either to the recovery of the subject without seizure manifestation at all, or to a gradual reduction of seizure burden and increase of interseizure intervals.

Surprisingly, but in line with clinical observations, remarkably long timescales of EPG (up to decades) emerge in our mathematical model, despite the “slowest” variables operating on relatively fast timescales of weeks. From the mathematical point of view, slowing down of dynamics is occurring in our model when the state of the system approaches the unstable steady state (Figure 1B). From the biological and clinical points of view, this would correspond either to the prolongation of the latent period duration, or, in case of slowing down of dynamics after seizure manifestation, to an absence of state changes in a subject (e.g. steady seizure rate, absence of progressive neuronal loss, and so forth). In order to visualize this property, we have performed a simulation of the injury with lowered intensity, which leads to slower EPG (Figure 3E). In this case, the progression of the pathology takes longer owing to slowing down of state changes around the unstable fixed point (Figure 3F). On the other hand, the EPG will be facilitated when caused by injury of increased intensity. In sum, consistent with clinical data, our model can capture the effect of the intensity of the injury on the latent period duration, seizure burden, and the risk of EPG.

Variability of epileptogenesis risk and pathology severity

In addition to the dose-dependent effects of injury intensity, the variability of EPG risk and the severity of pathology is evident even in animals exposed to identical injury. Even in highly standardized conditions of animal experiments, the fraction of animals that develop seizures, and the seizure burden in seizing animals are varying despite identical parameters of induced injury. For example, according to the data from a study by Polascheck et al. (2010), in 12 rats treated with pilocarpine, two have not shown any seizures, while seizure frequency ranged from 1 to 72 in the remaining animals at the time point of 8 weeks after the SE. The variability in EPG outcome is likely to depend on, among others, variability in genetic and epigenetic features of animals in experiment, and variability in various factors when conducting the experimental procedures. Moreover, the brain is an intrinsically stochastic complex system (Deco et al., 2009). Therefore, we added this intrinsic stochasticity to our model allowing for the variability in outcomes of EPG even in identical animals exposed to identical injuries.

To test whether our model is able to account for the experimentally observed variability of EPG outcomes, we simulated infection-induced EPG in the TMEV model (Figure 4). Our simulation reproduces the characteristic temporal pattern where seizures manifest after second day post-infection (day 2.83 ± 0.13) and profoundly drop in frequency during the second week post-infection (Figures 4A and 4F). Moreover, the computational model captures the characteristic time course of neuroinflammation (Figure 4B) as well as neuronal death, which is characterized by the occurrence of macroscopically measurable neuronal damage on day 4 post-infection and its further progression and saturation from the second week post-infection (Figure 4C). Our results suggest that this characteristic plateau of neuronal loss progression originates from the attenuation of the neuroinflammatory response during week two post-infection (Figure 4D). Kirkman et al. (2010) have shown that neuronal death was significantly more abundant in animals that developed seizures versus those that did not. Consistent with this observation, neuronal loss in the model is correlated with seizure burden during the acute post-infection stage (Figure 4E).

Figure 4

Model captures key features of infection-induced EPG

(A) Comparison of characteristic seizure occurrence patterns from TMEV animal model data (left, Patel et al. (2017), and simulation (right). Patel et al. (2017) reported the total number of seizures per day aggregated over N = 11 animals together. For simulation, data are represented as mean ± SEM.

(B) Comparison of neuroinflammation time courses from the TMEV model (left, Patel et al. (2017), and simulation (right). Data are represented as mean ± SEM.

(C) Comparison of neuronal loss score progression from TMEV model (left, Kirkman et al. (2010), and simulation (right). Neuronal loss score for the simulation was computed using the masking procedure from (Kirkman et al., 2010). Data are represented as mean ± SEM. Masking procedure and its effect of “masking out” variability in the simulation results are explained in detail in supplemental item S4.

(D) Neuroinflammation course in individual simulations (N = 30). Red dashed line corresponds to the neurotoxicity threshold Θ. Light red area corresponds to the time window of injury induction.

(E) Neuronal loss one-month post-infection (day 35) is correlated with the severity of seizure burden in the acute phase (week one post-infection). Blue dots correspond to individual simulations. Blue line corresponds to linear regression fit with the coefficient of determination R2 = 0.65.

(F) Time sequences of seizure occurrence in individual simulations. Red bar corresponds to the time window of injury induction. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

Evaluation of seizure occurrence on a follow-up period of one-year post-infection (Figures 4F, 5A, and 5B) illustrates the variability of EPG outcomes in simulated animals. Among 30 simulations, nine did not exhibit any seizures within one week, while the remaining 21 exhibited seizure burdens of various severities. Our model suggests that even for (hypothetically) identical simulated animals exposed to identical injury, the EPG outcome is variable owing to the stochastic nature of spontaneous recurrent seizures. Mathematically, the stochastic nature of seizure generation induces noise in the EPG (Figure 5C), affecting disease progression and outcome (Figures 5A and 5B).

Figure 5

Variability of EPG outcomes in identical simulated animals exposed to identical injury originates from the stochastic nature of spontaneous seizures

(A) Distribution of seizure rate one year after infection for 30 simulated TMEV animals.

(B) Examples of seizure rate development in time for three simulations with different seizure burden outcomes on one-year post-infection (line color code is consistent in all subfigures). The raster plot on top illustrates the occurrence of seizures in time for corresponding simulations.

(C) EPG course for three simulations with different seizure burden outcomes on one-year post-infection illustrated in B-R domain. Distance between circle markers on EPG traces correspond to time intervals of 7 days. Overall visualization period is 1-year post-infection. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

Computational model accounts for complex and injury-specific features of epileptogenesis

EPG is often conceptualized as a two-stage process comprising a clinically silent latent period before the occurrence of a first spontaneous seizure and a subsequent seizure period. However, a growing body of evidence suggests that this view may be overly simplistic (Pitkänen et al., 2015). For example, in a perforant path stimulation rat epilepsy model, epileptiform activity is already observed just after the injury induction before the first seizure (Bumanglag and Sloviter, 2018), questioning the existence of a latent period. Indeed, EPG appears to exhibit various injury-specific features. In the following, we show how our model can account for qualitatively different types of EPG. We illustrate this for different types of EPG-inducing injuries including BBB leakage, infection, and SE (owing to pilocarpine administration) in Figures 3, 4, and 6, respectively. In the case of EPG induced by BBB leakage, the latent period approaches one week in duration (4.97 ± 0.33 days, Figures 3A and 3C). On the contrary, in the TMEV infection model, the occurrence of first spontaneous seizures takes place already after the second day after the viral infection (2.83 ± 0.13 days). Then, this early onset of seizures is followed by a period of profound decrease in seizure activity starting during the second week post-infection (Figures 4A and 4F).

Figure 6

Mathematical model captures injury-specific features of EPG in pilocarpine-induced SE animal model of epilepsy

(A) Comparison of microglial activation progression from the animal model (left, Brackhan et al. (2016), with neuroinflammation course in simulations (right). Data are shown with mean values (red dots) and SEM error bars. Gray area corresponds to the time window of injury induction.

(B) Comparison of neuronal loss progression from the animal model (left, Patel et al. (2017), with neuronal loss in simulations (right). Data are shown with mean values (blue dots) and SEM error bars. Gray area corresponds to the time window of injury induction.

(C) Comparison of neuronal loss progression from animal model data (left, Zhang et al. (2015), with neuronal loss in simulations (right) illustrated for three time points. For simulation results, data are shown with mean values (blue bars) and SEM error bars.

(D) Simulation results illustrate processes underlying the rise of seizure rate after injury despite the relative recovery of the BBB integrity. Orange, light blue, and black thin lines correspond respectively to extent of BBB disruption, degree of circuit remodeling, and seizure rate in individual simulations (N = 30). Solid lines correspond to prediction from the rate model. Gray area corresponds to the time window of injury induction.

(E) Neuroinflammation development in time with the indication of presumed phase of relative recovery characterized by absence of neurotoxicity in rate model prediction. Thin lines correspond to individual simulations (N = 30). Solid lines correspond to prediction from the rate model. Red dashed line corresponds to the neurotoxicity threshold Θ. Gray area corresponds to the time window of injury induction.

(F) Extent of neuronal loss development in time with the indication of presumed phase of relative recovery characterized by the absence of neurotoxicity in rate model prediction (annotation identical to caption in E.). Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

In the pilocarpine-induced SE model of EPG, gliosis and neuronal death are progressing rapidly during the first week after SE (Figures 6A and 6B). Data (Brackhan et al., 2016) suggest that the progression is slowing down during the second week after injury, reaching a plateau as indicated by a comparison of day 5 and day 14 post-SE (Figure 6). However, when looking at a longer timescale, a comparison of the immunoreactivity of neuron-specific nuclear protein (NeuN) in the hippocampus (Zhang et al., 2015) between days 7 and 60 suggests a profound progression of neuronal loss (Figure 6C). Our model captures this temporal pattern of the pathological development of gliosis and neuronal loss with a relative recovery after the acute neurological injury (Figures 6A and 6B) and further progression of pathology in the chronic phase (Figure 6C).

The modeling results suggest that despite the relative recovery of the BBB permeability, seizure burden grows owing to the gradual increase in the degree of circuit remodeling (Figure 6D). The pathological changes in circuity are happening in reaction to the remaining BBB leakage and ensuing damage to the neural population. The slowing of neuronal loss progression (Figure 6B) and later progression during the chronic phase (Figure 6C) are explained by the leveling off of the neuroinflammation after the initial injury and subsequent growth of neurotoxicity associated with the growth of seizure burden (Figure 6D). Thus, seizures take over the propellent role in pathology development from the initial neurological injury. This results in the emergence of the complex temporal pattern of pathology progression in the pilocarpine animal model.

Multicausality and degeneracy in epileptogenesis: Neuronal loss is sufficient but not necessary for inducing epileptogenesis

The role of neuronal loss in epilepsy development has been extensively debated. Massive neurodegeneration in the hippocampus, known as hippocampal sclerosis, is a common pathology of temporal lobe epilepsy and other epilepsy syndromes (Thom, 2014). Moreover, the extent of neuronal loss has been shown to be positively correlated with seizure frequency (Lopim et al., 2016). However, it is still an open question whether neuronal loss is a primary cause of EPG, its consequence, or both (Tasch et al., 1999; Kapur, 2003; Sendrowski and Sobaniec, 2013). In a study by Weissberg et al. (2015), EPG with recurrent seizures was triggered in mice by the induction of BBB disruption without evidence of neuronal loss. This indicates that neuronal loss may not be necessary for EPG. As discussed above, this phenomenon is readily captured by our model, which can produce EPG without cell death based on inflammation and BBB disruption alone (Figure 3F). Given these results, we wondered if cell death, while not being necessary for EPG induction, may still be sufficient for it.

Indeed, rigorous stability analysis of the model and performed simulations predict that neuronal loss alone can trigger the induction of EPG (Figure 7A). Specifically, neuronal loss causes slow remodeling of neural circuits, which gradually lowers the seizure threshold and increases the seizure rate. Mathematically, the presence of neuronal loss modifies the locations of the steady states (Figure 7B): the “healthy” steady-state and the unstable fixed point move toward each other, resulting in a non-zero seizure rate even when the system is resting in the “healthy” steady state. Further increase of neuronal loss leads to a bifurcation where the “healthy” steady-state collides with the unstable fixed point at a certain value of neuronal cell loss Dcritical≈ 0.41 (Figure 7B, for derivation see STAR Methods). For values of neuronal cell loss greater or equal than Dcritical, the development of progressive EPG toward the “epileptic” steady-state is inevitable (Figure 7A). From the biological and clinical perspectives, the subcritical extent of neuronal loss increases the chance of arrival of spontaneous seizure but does not necessarily lead to epileptogenesis onset. On the other hand, the supracritical extent of neuronal loss would cause inevitable epileptogenesis, which, however, may evolve on timescales of up to decades.

Figure 7

Model reveals that neuronal loss is sufficient but not necessary for EPG

(A) EPG progression in five simulations with supercritical (D = 1.0 > Dcritical≈ 0.41) and subcritical (D = 0.3 < Dcritical≈ 0.41) extents of neuronal loss. The raster plots above seizure rate traces indicate seizure times of each simulation.

(B) The effect of neuronal loss on the system stability is illustrated with state-space plots for rising extent of neuronal loss from left to right panels (D = 0; 0.3; ≈0.41; 1.0). The state-space plots were obtained by performing the timescale separation procedure (see STAR Methods). Filled circles correspond to “healthy” (black) and “epileptic” (gray) steady states. Empty circles correspond to unstable fixed (saddle) points. Crossed empty circles correspond to semistable (one of the eigenvalues is zero) fixed points. The subcritical extent of neuronal loss leads to a rise in the chance of spontaneous seizure arrival, while in the case of supracritical extent of neuronal loss gradual epileptogenesis becomes inevitable. For detailed analysis see STAR Methods. Red dashed lines correspond to the neurotoxicity threshold Θ.

(C) Average time until inevitable progressive EPG for different extents of neuronal loss obtained with rate model. The time of progressive EPG was heuristically calculated as the time from the start of the simulation to the time point of the neuroinflammation I(t) reaching 90% of the value corresponding to the “epileptic” steady state. Black dashed line corresponds to the critical extent of neuronal loss Dcritical≈0.41. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

The exact extent of neuronal loss determines the average time until the development of progressive EPG (Figure 7C). However, progressive EPG is also possible in simulated animals with a subcritical extent of neuronal loss owing to stochasticity (D < Dcritical, see simulations with D = 0.3 in Figure 7A). The finding that neuronal loss is sufficient for EPG initiation and progression, while not being necessary for EPG in other types of injuries (Figure 3), highlights the multicausal nature of EPG, where distinct processes may drive the process in isolation or in a convergent fashion. This is in line with a recent proposal that in degenerate systems (Edelman and Gally, 2001) multiple different pathological changes are sufficient but not necessary to cause hyperexcitability (Ratté and Prescott, 2016).

Simulation of therapeutic interventions reveals injury-specific targets and optimal time windows for treatment

Neuroimmune interactions are potential targets in the search for efficient treatments for pharmacoresistant epilepsy. However, not only the selection of prominent targets for intervention but also the timing and duration of interventions seem to matter. For example, the application of rapamycin, which is suspected to have an antiepileptogenic effect via restoration and strengthening of the BBB, over the period of 6 weeks after induced SE was efficient in the reduction of the number of animals developing seizures, seizure frequency, and extent of neuronal loss (van Vliet et al., 2012). In contrast, the treatment that continued only for 2 weeks had no positive effect over a 6 weeks observation period (Sliwa et al., 2012). Our model provides an opportunity for simulating various intervention strategies, allowing for the selection of target and time window and exploring the effects of multi-target interventions.

Specifically, our model shows that a permanent suppression of the effect of seizures on BBB integrity (intervention I in Figure 8A) prevents EPG in a simulated pilocarpine rodent model of epilepsy. The long-term effect on seizure rate is illustrated in Figure 8C. The extent of neuronal loss over the first half a week after injury does not differ among simulations with and without interventions (Figure 8B). Figure 8D illustrates the impact of the suppression of the effect of seizures on BBB integrity. Under this suppression (right plot), only a single attractor corresponding to the “healthy” steady state remains in the state-space of the system, which makes gradual recovery of the system to the healthy state inevitable.

Figure 8

Modeling therapeutic intervention with suppression of seizure effect on BBB integrity reveals injury-specific time window for intervention in pilocarpine rodent model of epilepsy

(A) Time windows for six various interventions and a reference simulation without intervention. Gray area corresponds to the time window of injury induction. The suppression of seizure effect on BBB integrity is simulated with the 100-fold decrease of respective model variable K↓. Simulations are performed using the rate model.

(B) Neuronal loss progression in simulated animals exposed to various types of intervention. Gray area corresponds to the time window of injury induction.

(C) Seizure rate development in simulated animals exposed to various types of intervention.

(D) State-space plots illustrating the state of the system after the injury offset (D = 0.2) without (left) and under the effect of intervention (K↓, right). The state-space plots were obtained by performing the timescale separation procedure (see STAR Methods). Filled circles correspond to “healthy” and “epileptic” steady states. Empty circle corresponds to the unstable fixed point. Red dashed lines indicate the neurotoxicity threshold Θ. Without suppression, the system has two attractor states, while under the effect of the intervention, only a single “healthy” attractor remains.

(E) Response of the system to injury in simulated animals exposed to various types of intervention illustrated in the B-R domain. Red lines correspond to the reference simulation without intervention. Solid black lines starting with “x”-symbol and ending with a star correspond to the time interval of injury induction. Solid color lines starting with “x”-symbol and ending with a star correspond to time windows of intervention. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

Furthermore, we also investigated transient suppression of the effects of seizures on BBB integrity. Suppression during a long time window of 10 weeks (intervention II in Figure 8A) is sufficient to prevent EPG (Figures 8C and 8E), while a shorter window of 2 weeks (intervention III in Figure 8A) does not suffice (Figures 8C and 8E).

Interestingly, not only the duration of an intervention but also its precise timing relative to the injury is crucial for the successful prevention of EPG. For example, in the simulated pilocarpine rodent model of epilepsy, suppression of the effects of seizures on BBB integrity for 5 weeks can prevent EPG when applied 2 weeks after the injury (intervention V in Figure 8A). However, identical interventions starting at 0 or 5 weeks delays (interventions IV and VI in Figure 8A) are inefficient (Figures 8C and 8E).

Moreover, the time windows, during which interventions should be applied, have injury-specific durations and timings. For example, in the simulated TMEV infection rodent model, a 1-week long suppression of the effects of seizures on BBB integrity (intervention II in Figure 9A) is sufficient to prevent EPG (Figure 9C). Interestingly, in this model, the intervention time window has to overlap with the time window of the injury effects. Interventions that are applied with 1 and 2 weeks delays (interventions III and IV in Figure 9A) do not prevent EPG (Figure 9C). In contrast to the simulation of the pilocarpine rodent model (Figure 8B), in the TMEV infection model interventions that are applied during and after the period of injury effects lead to different levels of neuronal loss after the injury offset (Figure 9B). This results in different structures of the systems’ state spaces: in case of interventions I and II, which overlap with the injury time window (Figure 9A), a lower extent on neuronal loss leads to the preservation of a larger basin of attraction of the “healthy” steady state (Figure 9D). This is in contrast to interventions III and IV (Figure 9A), for which the “healthy” steady state and the separatrix are in closeness, proximity (Figure 9E). Therefore, interventions that are not applied during the first week after injury onset are insufficient for the prevention of EPG (Figure 9C).

Figure 9

Modeling therapeutic interventions suppressing the effects of seizures on the BBB integrity and activation of glia by factors infiltrating the parenchyma reveals injury-specific optimal time windows for intervention in the TMEV infection rodent model of epilepsy

(A) Time windows of four interventions suppressing the effects of seizures on BBB integrity. REF indicates reference without intervention. Light red area corresponds to the period of injury induction. The suppression of the effect of seizures on BBB integrity is simulated with a 100-fold decrease of the respective model variable K. Simulations are performed using the rate model.

(B) Neuronal loss progression for the different time windows of K reduction. Light red area corresponds to the period of injury induction.

(C) Seizure rate development for different intervention time windows.

(D) State-space plots illustrating the state of the system after the injury offset (D = 0.34) for interventions that coincided with injury time window without (left) and under the effect of intervention (right). The state-space plots were obtained by performing the timescale separation procedure (see STAR Methods). Filled circles correspond to “healthy” and “epileptic” steady states. Empty circle corresponds to the unstable fixed point. Red dashed lines correspond to the neurotoxicity threshold Θ. Without suppression, the system has two attractor states, while under the effect of the intervention, only a single “healthy” attractor remains.

(E) Same as D but for a higher value of D = 0.38.

(F) Time windows of four interventions with suppression of the activation of glia by factors infiltrating the parenchyma. Light red area corresponds to the period of injury induction. The suppression of activation of glia by factors infiltrating the parenchyma is simulated with 100-fold decrease of the respective model variable κ. Simulations are performed using the rate model.

(G) Neuronal loss progression for different time windows of κ reduction. Light red area corresponds to the time window of injury induction.

(H) Seizure rate development for different time windows of κ reduction. Model and simulation parameters are specified in supplemental items S2 and S3, respectively.

In the next step, we were interested in investigating the effects of different types of interventions in the TMEV infection model of epileptogenesis. We have simulated an intervention that suppresses glial ability to be activated by infiltrating blood factors with a 100-fold decrease of the respective model variable κ↓ (Figure 9F). This type of intervention requires a much longer time window of application but also requires application during the first week after injury onset (Figures 9F–9H). This highlights the injury-specificity of preferred intervention timing. In sum, our model can be used as a framework for the simulation of intervention strategies. It provides means for studying the efficiency of various therapeutic targets, intensities, and time intervals of intervention in an injury-specific manner.

Discussion

The pathophysiology of EPG is associated with the activation of innate and adaptive immune responses, disruption of BBB integrity, neuronal loss, circuit remodeling, and various other processes acting over a range of different timescales. Furthermore, the injury-specific time courses of clinical markers point to a bewildering complexity of the disease. This makes understanding EPG and the development of effective treatments a formidable challenge. Computational and mathematical modeling can be a valuable tool in understanding such complex systems. Indeed, other epilepsy-associated phenomena, such as ictogenesis, have already been successfully studied with mathematical modeling methods (Fröhlich et al., 2008; Jirsa et al., 2014, 2017; Proix et al., 2017; Gonzalez et al., 2018; González et al., 2019). Here, we have presented the first-of-its-kind mathematical model of EPG in the context of acquired epilepsy. Our model accounts for a wide range of EPG phenomena and is a tool for testing different interventions in silico, while generating testable predictions regarding their effectiveness. The model describes the interaction between neuroinflammation, BBB disruption, neuronal loss, circuit remodeling, and seizures in response to neurological injury. Mathematically, the model consists of a system of coupled stochastic nonlinear ordinary differential equations. Our formal analysis of the model has revealed the existence of two stable fixed points, representing the healthy state and the state of developed epilepsy.

Our model captures EPG being triggered by very different types of neurological injuries. Here, we have focused on three such injuries: infection as represented by a TMEV rodent model; chemical intoxication as represented by pilocarpine SE rodent model; and BBB leakage as represented by BBB disruption rodent model. We have found our model to be in good agreement with the experimental data from these animal models using a single set of parameters for all simulations.

The model captures injury-specific characteristics of EPG such as temporal patterns of seizure occurrence, the progression of neuronal loss, neuroinflammation, and BBB disruption. Interestingly, the model suggests an explanation for the emergence of long timescales (years and decades) of disease development despite time-limited injuries that directly affect the CNS for much shorter durations (days). Mathematically, these unexpectedly slow timescales of EPG are explained by a slowing of the system’s dynamics in the vicinity of an unstable fixed point. This resembles the emergence of slow dynamics in, e.g., wound healing after injury with paradoxically long scar formation (Adler et al., 2020). Moreover, our model describes the dependence of the latent period duration, the seizure burden, and the risk of EPG on the intensity of an injury — the dose-dependence effects of injury intensity observed in various human and animal models of EPG.

Our model also captures the multicausal nature of epilepsy. For example, our model describes how, on the one hand, neuronal loss alone may be sufficient to induce EPG, but, on the other hand, it is not at all necessary for EPG. This is in agreement with a recent observation that, in neuronal systems with degenerate mechanisms, several distinct pathologies are sufficient but not necessary to account for hyperexcitability (Ratté and Prescott, 2016). Furthermore, our model suggests that the variability of EPG outcomes originates in part from the stochastic nature of epileptic seizures, which can push the system from the basin of attraction of the “healthy” steady state to that of developed epilepsy.

In order to demonstrate the utility of the model for generating testable predictions of therapeutic interventions, we performed simulations with different intervention targets and time windows for different initial injuries. Our results suggest that therapeutic interventions applied during only a short but critical time window may be just as effective as long-term interventions. Moreover, the optimal time windows for interventions are injury-specific. For example, in the case of a TMEV infection model, the intervention has to be applied during the first week after injury onset and without a delay. In the pilocarpine model, in contrast, a 5-week intervention starting after 2 weeks prevented EPG, while earlier or later interventions were ineffective.

Limitations of the study

Owing to its simplicity, our model also has a number of important limitations. For example, the complex process of neuroinflammation is described by just a single “coarse-grained” variable. This aids mathematical analysis, but it complicates the interpretation of simulation results. Furthermore, the model does not distinguish different seizure types (e.g. focal vs generalized) nor does it allow for an evolution of seizure severity and duration throughout EPG. Besides processes described by the model, phenomena such as channelopathies, neurogenesis, gene transcription, epigenetic modifications, and others are also associated with EPG, but yet to be accounted for in our modeling framework.

Also, the model does not include any positive antiepileptic aspects of neuroinflammation, i.e. processes aimed at the maintenance of healthy CNS function. Future extensions of the model should take into account such protective aspects of neuroimmune interactions. This will allow for computational modeling of paradoxical phenomena such as inflammatory preconditioning and epileptic tolerance. These, together with further exploration of injury-specific targets and time windows for therapeutic interventions, are promising directions for future work. Finally, while we view it as a strength of the model that it accounts for the data from different animal models with a single set of parameters, we acknowledge that there are always inter-species and inter-individual variabilities of physiologic parameters. It will, therefore, be interesting to investigate how parameter variations change the individual susceptibility to EPG and its trajectory. Such an understanding will facilitate the development of individualized interventions in the spirit of precision medicine.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Danylo Batulin (batulin@fias.uni-frankfurt.de).

Materials availability

This study did not generate new unique reagents.

Method details

All simulations, analysis and figure design were performed with Jupyter notebooks in Python version 3.7.4 (for details, see Key resources table).

Mathematical description of the model

The interactions between the processes specified in Figure 1A are modeled with a system of stochastic nonlinear ordinary differential equations:where I(t) is neuroinflammation intensity; B(t) is the extent of BBB disruption; D(t) is the extent of neuronal loss; R(t) is the degree of circuit remodeling. All variables are assumed to be reversible, except for neural loss D(t), which is motivated by the impossibility of recovery of dead neurons and thereby excludes the possibility of neurogenesis in the adult nervous system. Despite the degree of circuit remodeling R(t) being a reversible variable, the presence of an irreversible neuronal loss (D(t) > 0) would cause an irreversible circuit remodeling. This can be shown considering the steady state of Equation 1: R˙ = 0, so R = κB+κD; if D > 0 then consecutively R > 0. κ,κ,κ,κ,κ are parameters for coupling strengths of the respective variables. The described processes are assumed to operate on 3 timescales: fast (seconds-minutes) for epileptic seizures; intermediate (hours-days) for neuroinflammatory reaction (τ); and slow (days-weeks) for permeability of the BBB (τ), neuronal loss (τ) and circuit remodeling (τ).

Mild neuroinflammation, which is a physiological process aiming to maintain tissue homeostasis (Yong et al., 2019), is assumed to have no neurotoxic effects. Thus, neurotoxicity leading to neuronal loss requires excessively activated glia (I(t)>Θ), where Θ is a neurotoxicity threshold and Dmax is a maximum possible extent of neuronal loss.

The term S(I,R) in Equation 1 describes an effect of seizure activity on permeability of the BBB:where κ=K/λmaxTseizure describes the effect of a single spontaneous recurrent seizure on the permeability of the BBB; λmax is a maximum possible amount of seizures per day; K is a parameter defining the maximum possible burden of seizure activity; Tseizure is the seizure duration.

The occurrence of spontaneous recurrent seizures in time is modeled with a Poisson process. This approach has been chosen as the simplest way to account for stochasticity in the process of spontaneous recurrent seizure generation. The seizure rate (the probability of a seizure occurring per unit time) is assumed to be monotonically increasing with the intensity of neuroinflammation I(t) and the extent of circuit remodeling R(t) according to:where κ and κ are parameters scaling the seizure-promoting contribution of, respectively, neuroinflammation and circuit remodeling. The sigmoid shape of the function (supplemental item S5) reflects the saturation effect of maximum possible seizure burden that nervous system may be exposed to within a finite time interval due to metabolic constraints. The assumption of a quadratic dependence of the seizure rate on the intensity of neuroinflammation minimizes seizure-promoting effects of mild neuroinflammation I(t)≳0, which per se is a physiological reaction to noxious stimuli.

Time scale separation and rate model

In addition to the model simulating stochastically occurring seizures (Equations 1 and 2), we developed a rate model, where the Poisson process is approximated with a seizure burden function:

The rate of seizures dictated by Equation 4 does not allow for tracking the occurrence of individual seizures, but it provides a means for more intuitive explanation of the dynamics of the system. Moreover, we can perform the time scale separation procedure for the equation describing neuroinflammation since its timescale is smaller than the timescales of the other processes ().

We assume that the process of neuroinflammatory reaction evolves faster than blood-brain barrier (BBB) disruption and recovery of its integrity, neuronal loss, and circuit remodeling: τ≪τ,τ,τ. Under the condition of an absent neuroinflammatory external input I, we can thus perform a time scale separation. At equilibrium, this will result in I≈κB and the system described in Equation 1 becomes:

From Equation 5, we can obtain a system of equations for fixed values of neuronal loss extent D = Dconst, where 0 ≤ Dconst ≤ Dmax. The resulting system of equations describes the system in the dynamical regimes characterized by the absence of neurotoxicity I≈κB<Θ:

By substituting I = κB in the equation for the extent of BBB disruption, we obtain the system described in B−R dimensions. It is used for analysis and visualization of dynamics with state-space plots for variables B and R:where I = κB<Θ and D = Dconst.

After performing the time scale separation , the fast evolution of the neuroinflammation variable can be approximated by the dynamics of the BBB disruption variable according to Equation 1: I(t)≈B(t). Thus, we can obtain the state-space representations of the model in the B−R domain for particular values of the monotonically rising variable D(t). Stability analysis (STAR Methods) shows that in the absence of neuronal loss, the system is bistable, having 3 steady states: a ‘healthy’ steady state in the origin; an unstable fixed point; and a stable steady state corresponding to the state of progressed EPG (Figure 1B). A separatrix, illustrated with a black dashed line, passes through the unstable fixed point and separates the basins of attraction of the two stable steady states. The neurotoxicity threshold (for neuronal death), illustrated with the red dashed line, divides the state-space into two areas: to the left are states in which no neuronal loss is being induced, and to the right are states in which the neuronal population experiences neurotoxic effects due to glial overactivation.

Stability analysis

In a steady state, B˙ = 0 and R˙ = 0. From Equation 7 we obtain:

Substituting S(I,R) from Equation 4:

The steady states (fixed points) in the system are the result of an intersection of and , which results in (inserting the parameter values from supplemental item S2) two equations:

The intersection of the latter Equation 11 with the horizontal axis will give all that satisfy and . And corresponding values of can be found using Equation 10.

The Jacobian of the linearized system around each fixed point is:where ; ; ; .

The eigenvalues for each fixed point are the eigenvalues of the Jacobian evaluated at . Analysing the fixed point positions and their corresponding eigenvalues, we can describe the state-space of the system:

when the system has 3 fixed points (Figure S1): a stable steady state (negative eigenvalues) around the origin; a saddle point (one positive and one negative eigenvalue); a stable steady state (negative eigenvalues) distanced from the origin in the first quadrant;

when the system undergoes the saddle node bifurcation and has 2 fixed points (Figure S1): a stable steady state (negative eigenvalues) distanced from the origin in the first quadrant, and a semistable point (one eigenvalue equal to 0) in the position of collision of two fixed points;

when the system has 1 fixed point - a stable steady state (negative eigenvalues) distanced from the origin in the first quadrant (Figure S1).

The code for numerical calculation of the fixed point values of and , and the eigenvalues used in the stability analysis can be found at GitHub (see Key resources table for DOI). For the calculation of bifurcation parameter value (the extent of neuronal loss at which the saddle node bifurcation is occurring) see corresponding section of STAR Methods.

Model parameter choice

The parameter values (Table S1) were, whenever possible, directly extracted from the experimental data. For parameters that are not directly measurable, stability analysis of the system and manual parameter search were carried out.

The parameter, corresponding to the maximum daily seizure frequency, was directly inferred from the available experimental data (Libbey et al., 2008; Stewart et al., 2010; Polascheck et al., 2010; Bar-Klein et al., 2014; Weissberg et al., 2015; Loewen et al., 2016; Kim et al., 2017a, 2017b; Patel et al., 2017). Together with the seizure duration and the burden of a single seizure , these parameters are attributes of the stochastic version of the model, which correspond to a single parameter in the rate version of the model. The parameter, which is scaling seizure burden in the system, is estimated with the stability analysis described later in this section. Having estimated, we assumed an average seizure duration value minutes, which in the current version of the model only determines the minimum time discretization step, and obtained the value of the single seizure burden according to .

are the timescale parameters, which do not affect the stability of the system (see Equation 1), but determine the velocity of the dynamic changes. The characteristic timescales were assigned based on knowledge from biology: neuroinflammatory reaction carried out primarily by microglia is assumed to be faster ( day) than the rest of the processes ( days, ratio of timescales 1:10). This assumption was motivated by the data on microglial motility in extending its processes to injury within hours (Davalos et al., 2005) and showing morphological changes (e.g. enlarged soma and shortening of processes) already 48 hours after the status epilepticus (Andoh and Koyama, 2021). The other described processes (vasculature recovery, neuronal loss and debris removal, circuit remodeling) take longer time. However, we have carried out control simulations with alternative ratios of timescales to illustrate the effect of parameter change and confirm the quality of the model fit to the animal model data (Figure S2).

The maximum possible extent of neuronal loss is a parameter that constraints the parameter space. For convenience, it was chosen equal to unity (or 100%). The neurotoxicity threshold of overactivated glia was manually fit. The quality of the fit is illustrated in the comparison of simulation results with alternative parameter values (Figure S3).

By default, all parameters were set equal to unity. However, the presence of the positive feedback loop between and (Figure 1A) imposes the constraint on the and parameters, which can be derived considering the subsystem:

Applying the stability analysis methodology described in STAR Methods Section, we obtain the inequality . This considered, was set to the value 0.1 and the value of was estimated manually (Figure S4). Considering alternative values of , we have identified partial redundancy in the parameter space: several combinations of parameter values may lead to the emergence of identical stability landscapes (Figure S5).

Changing the parameters for the strength of seizure-promoting effects of neuroinflammation and circuit remodeling leads to modulation of seizure rate for given values of neuroinflammation intensity and extent of circuit remodeling (Figure S6).

The effects of and changes on stability of the system is illustrated in Figures S7A and S7B. The effect of the parameter change on the accuracy of the model fit to the experimental data are shown in Figure S7C.

The scaling parameter of neurotoxic effect of overactivated glia was manually fit to capture the temporal pattern of neuronal loss progression for the assumed neurotoxicity threshold value (Figure S8).

The other attribute of the neuronal-loss-related subsystem (right-hand side of Figure 1A) is the scaling parameter for the effect of neuronal loss on circuit remodeling , which has an effect on the stability of the system. It was fit to (i) keep the system sensitive to the presence of neuronal loss and (ii) allow for possibility of the healthy steady state existence in conditions of nonzero neuronal loss (Figure S9).

Despite the quality check of the fit produced with chosen parameters (Table S1) and comparison to simulations with alternative parameter values (Figures S2–S8), we can not rule out the existence of a parameter set that would produce a quantitatively better fit.

In order to address the question of whether the number of simulations for the stochastic version of the model was chosen adequately, we ran additional simulations (dose-dependence effects discussed in Section ’Dependence of EPG on injury intensity and emergence of long timescales’ and illustrated in Figures 3A and 3B) for alternative number of simulations in the cohort and compared the results with (Figure S10). No major differences were observed for increased number of simulations.

Estimation of critical extent of neuronal loss

In this section, we are going to calculate the critical extent of neuronal loss , at which ’healthy’ steady state collides with the unstable fixed point and only one stable steady state remains for . From equations describing the system state for fixed extent of neuronal loss (Equation 6) and seizure burden function (Equation 4), we derive steady state equation for BBB disruption:and steady state equation for circuit remodeling:from which we derive:which will be referred to as linear R.

Inserting the parameter values (supplemental item S2) into Equation 14:

Defining , we can derive B from Equation 17:

Defining from Equation 18, we obtain:

From Equation 19, we can derive:

Now, we replace f with :

Taking logarithms of both sides:

From Equation 22, we can obtain nonlinear R equation:

The intersections between linear R and nonlinear R will give us all the fixed points of the system. With the parameters defined in supplemental item S2, this system of equations always has at least one fixed point for . In addition to this fixed point, a saddle node bifurcation can emerge when two additional fixed points are generated as a result of a change of parameter in the equations. Assuming that D could play the role of such a bifurcation parameter, we need to find its critical value such that linear R becomes tangential to nonlinear R. Decreasing this critical value would result in the emergence of two fixed points; however, increasing this value beyond would result in no intersection between the nullclines, and hence the system will have only one fixed point which was defined before. To find the value of in linear R which results in a tangent line to the nonlinear curve, first, we need to find the first derivative of nonlinear R with respect to B. Secondly, we should find all those values , where s is equal to the slope of linear R, or in other words, in linear R that is equal to (Equation 16, supplemental item S2). This indicates that we should find in for the nonlinear R. Finding the first derivative of Equation 23:

Equating Equation 24 to , we obtain a polynomial:

Solving this polynomial numerically (code available at GitHub, see Key resources table for DOI), we obtain the following values for  = [-1.259; 0.7448; 0.01439]. Since the BBB disruption variable can not take negative values, the solution [-1.259] is discarded. Using the equation for nonlinear R (Equation 16), we can calculate the corresponding R values for  = [0.7448; 0.0143]:

= [0.4560; 0.0146].

Note that these values of should hold in both nonlinear R and linear R, since they are the result of intersections between the nullclines. From linear R (Equation 16), we can derive the equation for :

For values of and , we can calculate  = [-577.5155; 0.4103]. The extent of neuronal loss can not take negative values, so we have to discard one of the solutions [-577.5155]. Thus, we have found the only critical extent value of neuronal loss  = 0.4103 0.41.

Simulation of various neurological injuries

In this paper, we carried out the simulation of three common animal models of epileptogenesis.

The first animal model is the TMEV mouse model, which is commonly used in epilepsy research for modeling infection-induced epilepsy in humans (Libbey et al., 2008, 2011; Stewart et al., 2010; Gerhauser et al., 2019). During the first week post infection, the acute neuroinflammatory response is developing. It is characterized by the presence of excessive concentrations of proinflammatory molecules (Patel et al., 2017), micro- and astrogliosis (Kirkman et al., 2010), which are followed by relative recovery during subsequent weeks post infection. Taking this dynamics of pathology development into account, we simulate infectious injury by induction of neuroinflammation with onset and offset within the first week post infection (Figure 2A).

The second animal model is a pilocarpine rodent model, which is based on induction of SE (excessively long generalized seizure) by injection of pilocarpine with further pharmacological termination of SE (Polascheck et al., 2010; Zhang et al., 2015; Brackhan et al., 2016; Kim et al., 2017a, 2017b). SE induces neuronal loss, neuroinflammation and profound leakage of the BBB. In our framework, we do not define the inflammatory perturbation for simulation of pilocarpine since it would be indirectly induced by BBB disruption. Thus, we simulate the induction of pilocarpine model injury as a combination of two external perturbations (Figure 2B): neuronal death, which allows us to account for initial SE-associated neuronal loss (Auvin et al., 2010a, 2010b), and BBB leakage, which normalizes after 1–2 days post SE (Bankstahl et al., 2018).

The third animal model is based on the induction of BBB leakage. It is obtained via exposure of the brain tissue to an artificial cerebrospinal fluid containing bile salts (Seiffert et al., 2004; Tomkins et al., 2007). Alternatively, an artificial cerebrospinal fluid may contain serum albumin, which mimics the extravasation of this protein in the brain parenchyma when BBB is dysfunctional (Weissberg et al., 2015). This animal model is simulated by setting B to a high value for a period corresponding to the time of application of bile salts or pumping albumin into the brain of an animal (Figure 2C).

Animal model data

BBB disruption rodent model data used in Figures 3A and 3B:

The following data from Weissberg et al. (2015) were used in this study:

latent period duration of 4.9 1.3 days (mean SEM), N = 10;

spontaneous seizures frequency of 1.16 0.16 seizures per day (mean SEM), N = 10.

Theiler’s murine encephalomyelitis virus (TMEV) mouse model data used in Figures 4A–4C.

The following data from Patel et al., 2017 were used in this study:

number of seizures per day for N = 11 mice was extracted from Figure 2 (Patel et al., 2017). Average seizure frequency per mice was calculated for 3 time intervals: day 1 post infection, days 2–7 post infection and days 8–15 post infection;

TNF protein fold change (relative to PBS-injected control mice) on day 1 post infection (N = 8): 6.9 0.6, day 5 post infection (N = 6): 206.2 14.9, day 14 post infection (N = 5): 34.8 7.1. Data are presented in mean SEM

The following data from Kirkman et al. (2010) were used in this study:

neuronal cell loss score for 2 hippocampi (mean SEM) on days 1–35 post infection from Figure 2 (Kirkman et al., 2010), N = 4–13 per time point group.

Chemically-induced (pilocarpine) SE rodent model data used in Figures 6A–6C:

The following data from Brackhan et al. (2016) were used in this study:

microglial activation score for the hippocampus (mean SEM) on days 0, 2, 5, 14 post SE from Figure 4 (Brackhan et al., 2016), N=3-5 per time point group;

neuronal cell loss score for the hippocampus (mean SEM) on days 0, 2, 5, 14 post SE from Figure 4 (Brackhan et al., 2016), N=3-5 per time point group.

The following data from Zhang et al. (2015) were used in this study:

NeuN-immunoreactive cells count per in the hippocampus of pilocarpine treated animals from Figure 5 (Zhang et al., 2015). Fraction of cells missing (in %) was computed for days 7 and 60 after pilocarpine injection relatively to values for untreated animals.

Quantification and statistical analysis

All experiments were performed as mathematical model simulations. Group allocation of samples is described in Table S2. No data were excluded as outliers. For stochastic model simulations, the sample size of N = 30 was chosen, which is twice larger than a common sample size in animal model experiments (Weissberg et al., 2015; Patel et al., 2017; Kirkman et al., 2010; Brackhan et al., 2016; Zhang et al., 2015). Data are presented both in raw format (e.g. neuroinflammation intensity development over time, time sequences of seizure occurrence), and in format of mean standard error of the mean (SEM), when convenient (e.g. latent period duration, seizure occurrence frequency for the whole sample). Due to bistability of system states (Figure 1B, STAR Methods), data could not be assumed to have Gaussian statistics and non-parametric tests were used. Specifically, an unpaired two-group Mann-Whitney U test was performed for analyses of dose-dependence effects on the intensity of the injury. was chosen as a threshold for statistical significance and exact p values were reported.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Deposited data
Data from BBB disruption rodent model of epileptogenesis	Weissberg et al., 2015	https://doi.org/10.1016/j.nbd.2015.02.029
Data from Theiler’s murine encephalomyelitis virus mouse model of epileptogenesis	Patel et al., 2017	https://doi.org/10.1523/ENEURO.0105-17.2017
Data from Theiler’s murine encephalomyelitis virus mouse model of epileptogenesis	Kirkman et al., 2010	https://doi.org/10.1111/j.1528-1167.2009.02390.x
Data from chemically-induced (pilocarpine) SE rodent model of epileptogenesis	Brackhan et al., 2016	https://doi.org/10.2967/jnumed.116.172494
Data from chemically-induced (pilocarpine) SE rodent model of epileptogenesis	Zhang et al., 2015	https://doi.org/10.7150/ijms.10527
Software and algorithms
Python version 3.7.4	Python Software Foundation	https://www.python.org
Jupyter Notebook	The Jupyter Software	https://jupyter.org/
Data generation, analysis and figure production scripts	This paper	https://doi.org/10.5281/zenodo.6035834https://github.com/danylodanylo/math-model-epileptogenesis