# Neural mass modeling of slow-fast dynamics of seizure initiation and abortion

Elif Köksal Ersöz, Julien Modolo, Fabrice Bartolomei, Fabrice Wendling

► To cite this version:

Elif Köksal Ersöz, Julien Modolo, Fabrice Bartolomei, Fabrice Wendling. Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLoS Computational Biology, 2020, 16 (11), pp.e1008430. 10.1371/journal.pcbi.1008430. hal-02998818v2

HAL Id: hal-02998818

https://hal.science/hal-02998818v2

Submitted on 1 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Neural mass modeling of slow-fast dynamics of seizure initiation and abortion

Elif Köksal Ersöz¹, Julien Modolo¹, Fabrice Bartolomei²,³, Fabrice Wendling¹*

1 University of Rennes, Inserm-U1099, LTSI, Rennes, France, 2 Aix Marseille University, Inserm, INS, Institut de Neurosciences des Systèmes, Marseille, France, 3 APHM, Timone Hospital, Clinical Neurophysiology, Marseille, France

* fabrice.wendling@inserm.fr

<figure><img src="image_2.png" /></figure>

## Abstract

Epilepsy is a dynamic and complex neurological disease affecting about 1% of the world-wide population, among which 30% of the patients are drug-resistant. Epilepsy is characterized by recurrent episodes of paroxysmal neural discharges (the so-called seizures), which manifest themselves through a large-amplitude rhythmic activity observed in depth-EEG recordings, in particular in local field potentials (LFPs). The signature characterizing the transition to seizures involves complex oscillatory patterns, which could serve as a marker to prevent seizure initiation by triggering appropriate therapeutic neurostimulation methods. To investigate such protocols, neurophysiological lumped-parameter models at the mesoscopic scale, namely neural mass models, are powerful tools that not only mimic the LFP signals but also give insights on the neural mechanisms related to different stages of seizures. Here, we analyze the multiple time-scale dynamics of a neural mass model and explain the underlying structure of the complex oscillations observed before seizure initiation. We investigate population-specific effects of the stimulation and the dependence of stimulation parameters on synaptic timescales. In particular, we show that intermediate stimulation frequencies (>20 Hz) can abort seizures if the timescale difference is pronounced. Those results have the potential in the design of therapeutic brain stimulation protocols based on the neurophysiological properties of tissue.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Köksal Ersöz E, Modolo J, Bartolomei F, Wendling F (2020) Neural mass modeling of slow-fast dynamics of seizure initiation and abortion. PLoS Comput Biol 16(11): e1008430. https://doi.org/10.1371/journal.pcbi.1008430

**Editor:** Peter Neal Taylor, Newcastle University, UNITED KINGDOM

**Received:** May 28, 2020

**Accepted:** October 8, 2020

**Published:** November 9, 2020

**Copyright:** © 2020 Ersöz et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting Information files. Simulation files are available from the GitHub database (https://github.com/elifkoksal/NMM_BurstingDynamics).

**Funding:** We acknowledge the funding support from NINDS Grant R01 NS092760-01A1 (EKE). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

## Author summary

Epilepsy is a complex disease affecting 1% of the worldwide population of which 30% of the patients are drug-resistant and seeking for alternative therapeutic methods, such as neurostimulation. Epileptic seizures can be hallmarked by preceding pre-ictal phases which are a possible window of opportunity to trigger electrical stimulation with the objective to prevent seizure initiation. Biophysiological models are an appropriate framework to understand underlying dynamics and transitions between different epileptogenic phases. In this study, we consider a typical pre-ictal regime with complex bursting-type oscillations, which can be accurately reproduced by a neural mass model. By analyzing the multiple time-scaled structure of the model, we identify the key role of the subpopulations of GABAergic interneurons. We show that appropriate brain stimulation targeting
GABAergic interneurons is able to abort pre-ictal bursting, thus preventing seizures to develop.

Introduction

Epilepsy is a severe, multi-causal chronic disease defined by the recurrence of unpredictable seizures that severely affect patients' quality of life. In 30% of patients, antiepileptic drugs [1] remain inefficient to control the occurrence of seizures. In most cases, drug-resistant epilepsies are 'focal' [2], as characterized by an epileptogenic zone (EZ) that is relatively circumscribed in one of the two cerebral hemispheres. There is a large body of evidence supporting that the balance between excitatory and inhibitory processes is modified in the EZ [3] due to multiple, not mutually exclusive, pathological mechanisms resulting from changes occurring at the cellular level (e.g. hyperexcitability caused by potassium and chloride dysregulation, review in [4]), up to the network level (e.g. hyperexcitability caused by altered glutamatergic or GABAergic synaptic transmission, review in [5]). Unfortunately, surgical treatment can only be offered to 15-20% drug-resistant patients [6] in whom the benefit-to-deficit ratio is favorable. Therefore, alternative therapeutic procedures aimed at reducing seizures' frequency are urgently needed.

Among these procedures, direct electrical stimulation of the brain is an increasingly popu-
lar technique of treating epilepsy, as evidenced by both animal and human studies [7]. Stimula-
tion targets have included deep brain structures such as thalamic nuclei, hippocampus or
cortical targets [8]. It has been acknowledged for decades that stimulation of the cortex during
routine brain mapping procedures may induce epileptiform discharges or seizures, but more
recently pulse trains have demonstrated their potential in aborting abnormal epileptiform
activity [9]. Direct stimulation has been shown to be effective in suppressing epileptic activity,
however with inconsistent results among patients. Furthermore, brain stimulation in drug-
refractory patients is recognized to be still largely empirical [10]. A rational definition of stim-
ulation protocols is indeed still missing, as evidenced in randomized controlled trials [11].

In this context, the specific objective of the present study is to exploit neuro-inspired mod-
els to design neurostimulation protocols aiming at aborting seizures at their onset. More spe-
cifically, we investigate a well-defined pattern of interictal-to-ictal transition characterized by
the occurrence of pre-ictal rhythmic large amplitude spikes followed by a fast onset activity, as
observed in stereo-EEG recordings (SEEG, intracerebral electrodes). Although not the unique
one, this commonly encountered pattern has long been considered as a hallmark of the EZ,
especially in mesial temporal lobe epilepsy [12–14]. More particularly, we focus on pre-ictal
bursting characterized by active episodes (fast epileptic spikes), repeated (quasi-) periodically
and separated by quiescent (slow-wave and/or silent) phases. First, we accurately reproduce
human electrophysiological patterns in a neural mass model (NMM) featuring glutamatergic
pyramidal neurons as well as two types of GABAergic interneurons (somatostatin-positive or
SOM+, and parvalbulmin-positive or PV+). After integrating neurostimulation effects in the
model as a parametrizable exogenous membrane perturbation of the main cells and interneu-
rons, we analyze the slow-fast nature of this nonlinear dynamical system in the bursting regime
by using numerical bifurcation analysis and geometric singular perturbation theory (GSPT)
[15,16]. Following this approach, the mechanisms leading to the pre-ictal bursting are deter-
mined, and the perturbation effects are explained geometrically. Here, the perturbation corre-
sponds to a suprathreshold constant current stimulation that is stronger than that used for
neuromodulation. In the following, the word "perturbation" refers to this "strong perturba-
tion". Identification of the model structures to be targeted for bursting abortion highlight the
key role of SOM+ interneurons in suppressing pre-ictal epileptic activity. Overall, the results of this methodological study elucidate mathematically the nature of pre-ictal spike bursting and provide aspiring information to design optimal direct stimulation protocols targeting this specific epileptiform pattern. To the best our knowledge, this is the first study deciphering explicitly the multiple time-scale structure of a neurophysiologically grounded NMM and modelling pre-ictal bursting.

# Model and methods

## Model

Since the early seventies, and following the pioneer works of W. Freeman [17] and F. Lopez da Silva [18], NMMs have been extensively used to study LFPs not only in the context of physiological activity (such as brain rhythms [19], and/or visually evoked potentials [20]), but also pathophysiological activity (epilepsy [21–24], Alzheimer’s disease [25], Parkinson’s disease [26]). NMMs are average descriptions of the temporal dynamics of neuronal assemblies, and therefore model parameters are themselves lumped. To some extent, these models are complementary of models of microcircuits in which neurons are explicitly represented. The important point there is that NMMs are neuro-inspired, since they implement interactions between different subpopulations of neurons (glutamatergic-excitatory- neurons and GABAergic-inhibitory- interneurons of various types) and model parameters are physiologically grounded [27,28]. Typically, the time constants governing average post-synaptic potentials (PSPs) match the rising and decaying times of corresponding experimentally recorded PSPs [27]. Regarding average firing rates, sigmoid functions account for threshold and saturation effects that are major features of neuronal physiology [28] and account for non-linearities in neuronal dynamics. Importantly, NMMs have been used to simulate LFPs which closely resemble those experimentally recorded, and led to predictions which could be experimentally verified [29].

We consider the NMM presented in [22] which includes three interacting neuronal subpopulations: pyramidal neurons (divided into two subpopulations) and inhibitory interneurons (SOM+ and PV+, also called “dendrite-projecting slow” and “soma-projecting fast” interneurons, respectively) (see S1 Fig for the block diagram of the model). The average PSP of each subpopulation is determined by two functions: 1) a ‘pulse-to-wave’ function, $S(v) = 5/(1+\exp(0.56(6-v)))$, transforming the incoming PSPs into a firing rate; and 2) the input firing rate is converted into the mean PSP of the corresponding subpopulation by a linear transformation, that is $h(t) = Wt/\tau_w \exp(-t/\tau_w)$, where W represents the average synaptic gain and $\tau_w$ is the average synaptic time constant, which is related to rise and decay times of PSPs and lumps together the passive membrane time constant and all signal propagation delays. The system reads:

$$ 
\begin{aligned}
\dot{y}_0 &= \frac{A}{\tau_a} S(y_1 - y_2 - y_3) - \frac{2}{\tau_a} \dot{y}_0 - \frac{1}{\tau_a^2} y_0, \\
\dot{y}_1 &= \frac{A}{\tau_a} \{p(t) + C_2 S(C_1 y_0)\} - \frac{2}{\tau_a} \dot{y}_1 - \frac{1}{\tau_a^2} y_1, \\
\dot{y}_2 &= \frac{B}{\tau_b} C_4 S(C_3 y_0) - \frac{2}{\tau_b} \dot{y}_2 - \frac{1}{\tau_b^2} y_2, \\
\dot{y}_3 &= \frac{G}{\tau_g} C_7 S(C_5 y_0 - C_6 y_1) - \frac{2}{\tau_g} \dot{y}_3 - \frac{1}{\tau_g^2} y_3, \\
\dot{y}_4 &= \frac{B}{\tau_b} S(C_3 y_0) - \frac{2}{\tau_b} \dot{y}_4 - \frac{1}{\tau_b^2} y_4.
\end{aligned}
\quad (1) $$
<table><caption>Table 1. Parameter values during simulated background activity.</caption><thead><tr><th>A (mV)</th><th>B (mV)</th><th>G (mV)</th><th>p (Hz)</th><th>C<sub>1</sub></th><th>C<sub>2</sub></th><th>C<sub>3</sub></th><th>C<sub>4</sub></th><th>C<sub>5</sub></th><th>C<sub>6</sub></th><th>C<sub>7</sub></th><th>τ<sub>a</sub>(s)</th><th>τ<sub>b</sub>(s)</th><th>τ<sub>g</sub>(s)</th></tr></thead><tbody><tr><td>5</td><td>40</td><td>35</td><td>90</td><td>135</td><td>108</td><td>35</td><td>25</td><td>450</td><td>121</td><td>121</td><td>0.01</td><td>0.05</td><td>0.003</td></tr></tbody></table>

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

Variables $y_i$ stand for the PSPs generated at the level of pyramidal cells ($y_0$), excitatory inputs on pyramidal cells ($y_1$) (or in other words, the pyramidal subpopulation denoted by $y_0$ excites the pyramidal subpopulation $y_1$), SOM+ interneurons ($y_2$), PV+ interneurons ($y_3$), and inhibitory inputs on PV+ interneurons ($y_4$). Parameters A, B, G are the synaptic gains, $\tau_a$, $\tau_b$, $\tau_g$ are the synaptic time constants, connectivity constants $C_i$ represent the average number of synaptic contacts, and $p(t)$ is the external (noisy) cortical input $p(t) = p+\xi$, where $p$ is the mean of the external input, and $\xi$ is a random variable following a normal distribution with zero mean and standard deviation $\sigma$. Table 1 presents the parameter values used in this manuscript unless otherwise stated. The main difference between the parameter sets of [22] and Table 1 is the connectivity strengths of the circuit involving PV+ interneurons. Note that $\tau_a$ is 3.3 times and $\tau_b$ is 16.6 times greater than $\tau_g$. These differences would introduce multiple time-scale dynamics in the system. Below, we recall primaries of slow-fast analysis before expressing (1) in slow-fast formulation.

## Primaries of slow-fast analysis

Complex oscillations, e.g. spiking, bursting and mixed-mode oscillations, are widely present in electrophysiological signals, and are hallmarked by interacting components running in different timescales (for a recent review see [30]). Mathematical models describing these dynamics involve inevitably variables with different timescales. For example, the interaction of two variables running at two different time-scales (1 fast and 1 slow variables) is sufficient for generating spiking, whereas bursting emerges only in higher dimensional systems with at least two different time-scales (2 fast and 1 slow variables/subsystems). Geometric methods in the singular limits (denoted as slow-fast or multiple timescale analysis) simplify the investigation of the model dynamics by breaking into several reduced systems. For instance, bursting patterns can be classified in a singular limit where the slowest variables are treated as parameters [31,32].

A slow-fast system in the general slow form reads,

$$ \epsilon \dot{x} = f(x, z, \epsilon), $$

$$ \dot{z} = g(x, z, \epsilon), $$

with fast variables $x$ and slow variables $z$ of arbitrary dimensions, time scale parameter $0<\epsilon\ll 1$, and dot represents derivation with respect to time $t$. The dynamics of a slow-fast system can be divided into *fast* and *slow* epochs. Each of these epochs can be investigated with the slow-fast analysis in a hybrid manner and then can be concatenated, so that one can understand the underlying structure giving sharp transitions (excitable responses to external inputs) and complex oscillatory patterns (spiking, bursting and subthreshold oscillations) [33].

An important geometrical object for both the slow and the fast dynamics is the *critical manifold* $C^0$, defined as the nullcline of the fast variable $C^0 = \{(x, z) | f(x, z, 0) = 0\}$, eventually obtained by setting $\epsilon = 0$. For the differential-algebraic system defined for $\epsilon = 0$, the so-called *reduced system* (*slow subsystem*) approximates the slow dynamics of the original system. The critical manifold $C^0$ both defines the phase space of the reduced system and equilibrium points of the layer problem expressed in the fast time-scale, that is

$$ x' = f(x, z, 0), $$

$$ z' = 0, $$
where (') denotes derivative with respect to $\tau = t/\epsilon$. The stability of the layer problem determines the characteristics of the critical manifold. The critical manifold $C^0$ is normally hyperbolic along the set for which $\det(f_x(x,z,0)) \neq 0$, which can be attracting, repelling or saddle type. The Fenichel theory [15] guarantees that these normally hyperbolic points of the critical manifold perturb smoothly in $\epsilon$ and give slow manifolds ($C^\epsilon$) of the original system for small enough $\epsilon>0$. If $C$ is folded, attracting and repelling branches of $C^0$ meet along the fold set $\mathcal{F} = \{\det(f_x(x_{fold}, z_{fold}), 0) = 0\}$, where normal hyperbolicity is lost. Extension of the classical Fenichel theory to non-hyperbolic sets provides a tool to investigate the slow dynamics near $\mathcal{F}$, and one can expect canard solutions in the neighborhood of such sets [34].

## Slow-fast formulation of the model

One can notice that the variable $y_4$ in (1) is equivalent to $y_2$, thus the dimension of (1) can be reduced by multiplying the PSP variables with $C_i$s before the 'pulse-to-wave' conversion. Further, by applying the variable conversion,

$$ \left( \frac{y_0}{\tau_a}, \frac{y_1}{\tau_a}, \frac{y_2}{\tau_b}, \frac{y_3}{\tau_g}, y_5, y_6, y_7, y_8 \right) \rightarrow (v_0, v_1, v_2, v_3, v_4, v_5, v_6, v_7, v_8), $$

system (1) can be written as:

$$
\begin{align}
\tau_g \dot{v}_3 &= y_8, \nonumber \\
\tau_g \dot{v}_8 &= G S(C_5 \tau_a v_0 - C_6 \tau_b v_2) - v_3 - 2y_8, \nonumber \\
\tau_a \dot{v}_0 &= y_5, \nonumber \\
\tau_a \dot{v}_5 &= A S(A \tau_a p + C_2 \tau_a v_1 - C_4 \tau_b v_2 - C_7 \tau_g v_3) - v_0 - 2y_5, \nonumber \\
\tau_a \dot{v}_1 &= y_6, \nonumber \\
\tau_a \dot{v}_6 &= A S(C_1 \tau_a v_0) - v_1 - 2y_6, \nonumber \\
\tau_b \dot{v}_2 &= y_7, \nonumber \\
\tau_b \dot{v}_7 &= B S(C_3 \tau_a v_0) - v_2 - 2y_7. \tag{2}
\end{align}
$$

Intuitively, system (1), hence system (2), are multiple-time-scale systems which can result in complex epileptogenic patterns for appropriate choices of parameters. Thus, understanding the multiple-time-scale structure of (2) is indispensable for designing brain stimulation protocols aiming at aborting the aforementioned oscillatory patterns. In order to proceed a slow-fast analysis of (2), where the time-scaling is not explicit, and use the standard methods of GSPT, we first normalize time $t$ with respect to $\tau_g$, as $\tilde{t} = t/\tau_g$, then define two parameters, $\delta = \tau_g/\tau_a$ and $\varepsilon = \tau_a/\tau_b$. We further assume that $\varepsilon$ and $\delta$ are independent of the synaptic time constants ($\tau_a, \tau_b, \tau_g$), similar to the approach followed in [35], that is to say variations in $(\varepsilon, \delta)$ for theoretical analysis do not require varying $(\tau_a, \tau_b, \tau_g)$. At the end system (2) is expressed in an
explicit slow-fast formulation:

$$
\begin{align}
\frac{dv_3}{dt} &= y_8 := F_3(y_8), \nonumber \\
\frac{dy_8}{dt} &= G S(C_5 \tau_a v_0 - C_6 \tau_b v_2) - v_3 - 2y_8 := F_8(v_0, v_2, v_3, y_8), \nonumber \\
\frac{dv_0}{dt} &= \delta y_5 := \delta F_0(y_5), \nonumber \\
\frac{dy_5}{dt} &= \delta (A S(A \tau_a p + C_2 \tau_a v_1 - C_4 \tau_b v_2 - C_7 \tau_g v_3) - v_0 - 2y_5) := F_5(v_0, v_1, v_2, v_3, y_5), \tag{3} \\
\frac{dv_1}{dt} &= \delta y_6 := \delta F_1(y_6), \nonumber \\
\frac{dy_6}{dt} &= \delta (A S(C_1 \tau_a v_0) - v_1 - 2y_6) := \delta F_6(v_0, v_1, y_6), \nonumber \\
\frac{dv_2}{dt} &= \delta \epsilon y_7 := \delta \epsilon F_2(y_7), \nonumber \\
\frac{dy_7}{dt} &= \delta \epsilon (B S(C_3 \tau_a v_0) - v_2 - 2y_7) := \delta \epsilon F_7(v_0, v_2, y_7). \nonumber
\end{align}
$$

System (3) is a three-time-scale system for small enough values of $(\delta, \epsilon)$ [27–31], with $(v_3, y_8)$ being fast variables, $(v_0, y_5, v_1, y_6)$ slow variables, and $(v_2, y_7)$ superslow variables. System (3) is written using the (fast) time $\tilde{t}$, and called the *fast system*. As can be noticed, rescaling of (1) identifies the small parameters such that GSPT can be applied and expresses the importance of timescales. We follow [36,37] to analyze the three time scaled slow-fast structure of (3). Defining $\tilde{t}_s = \delta \tilde{t}$ gives the *slow system*:

$$
\begin{align}
\delta \frac{dv_3}{dt_s} &= F_3(y_8) \nonumber \\
\delta \frac{dy_8}{dt_s} &= F_8(v_0, v_2, v_3, y_8), \nonumber \\
\frac{dv_0}{dt_s} &= F_0(y_5), \nonumber \\
\frac{dy_5}{dt_s} &= F_5(v_0, v_1, v_2, v_3, y_5), \nonumber \\
\frac{dv_1}{dt_s} &= F_1(y_6), \nonumber \\
\frac{dy_6}{dt_s} &= F_6(v_0, v_1, y_6), \nonumber \\
\frac{dv_2}{dt_s} &= \epsilon F_2(y_7), \nonumber \\
\frac{dy_7}{dt_s} &= \epsilon F_7(v_0, v_2, y_7). \tag{4}
\end{align}
$$

where $F_8$ are as defined for (3), with $i$ representing the system variables' indices on the left-hand side. S2A Fig presents the bifurcation diagram of the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem of (4) for $\epsilon = 0$, where $v_2$ acts as a parameter, and a periodic bursting orbit of (3) for $\epsilon = 0.01$. We see that the orbit agrees with the bifurcation diagram when $\epsilon$ is decreased. Details of the bursting behavior are explained in Sec. *Bursting analysis*.
Defining $\tilde{t}_{ss} = \varepsilon\tilde{t}_s = \varepsilon\delta\tilde{t}$ gives the superslow system,

$$
\begin{align}
\varepsilon\delta \frac{dv_3}{d\tilde{t}_{ss}} &= F_3(y_8), \nonumber \\
\varepsilon\delta \frac{dy_8}{d\tilde{t}_{ss}} &= F_8(v_0, v_2, v_3, y_8), \nonumber \\
\varepsilon \frac{dv_0}{d\tilde{t}_{ss}} &= F_0(y_5), \nonumber \\
\varepsilon \frac{dy_5}{d\tilde{t}_{ss}} &= F_5(v_0, v_1, v_2, v_3, y_5), \nonumber \\
\varepsilon \frac{dv_1}{d\tilde{t}_{ss}} &= F_1(y_6), \nonumber \\
\varepsilon \frac{dy_6}{d\tilde{t}_{ss}} &= F_6(v_0, v_1, y_6), \nonumber \\
\frac{dv_2}{d\tilde{t}_{ss}} &= F_2(y_7), \nonumber \\
\frac{dy_7}{d\tilde{t}_{ss}} &= F_7(v_0, v_1, y_7). \tag{5}
\end{align}
$$

Systems (3), (4) and (5) are equivalent if $\varepsilon \neq 0$ and $\delta \neq 0$, but they give nonequivalent dynamics in the singular limits $\varepsilon \to 0$ and/or $\delta \to 0$. The limit $\delta \to 0$ in the fast system (3) eliminates the slow and superslow dynamics and yields the fast layer problem,

$$
\begin{align}
\frac{dv_3}{d\tilde{t}} &= F_3(y_8), \nonumber \\
\frac{dy_8}{d\tilde{t}} &= F_8(v_0, v_2, v_3, y_8), \tag{6}
\end{align}
$$

which describes the dynamics of the fast variables $(v_3, y_8)$ for fixed values of $(v_0, v_2)$, $(v_0^0, v_2^0)$ for instance. The critical manifold is defined by the four-dimensional set of equilibria of the fast layer problem (6), which reads,

$$S^0 = \{(v_3, y_8, v_0, v_2) | F_3(y_8) = 0 \cap F_8(v_0, v_2, v_3, y_8) = 0\},$$

and $S^0$ is eventually in the $(y_8 = 0)$-space. The stability of $S^0$ is determined by deriving the Jacobian of $S^0$ with respect to the fast variables, that is,

$$Jac(S^0_{v_3, y_8}) = \begin{bmatrix} 0 & 1 \\ -1 & -2 \end{bmatrix}.$$

Since $\det(Jac(S^0_{v_3, y_8})) \neq 0$, and the eigenvalues are $\lambda_{1,2} = -1$, the $S^0$ is normally hyperbolic and stable. Hence, $S^0$ is perturbed to local invariant slow manifolds for sufficiently small $\delta>0$.
Another singular limit is obtained by letting $\delta \to 0$ in the slow system (4) gives the algebraic-differential slow reduced problem,

$$
\begin{align*}
0 &= F_3(y_8), \\
0 &= F_8(v_0, v_2, v_3, y_8), \\
\frac{dv_0}{dt_s} &= F_0(y_5), \\
\frac{dy_5}{dt_s} &= F_5(v_0, v_1, v_2, v_3, y_5), \\
\frac{dv_1}{dt_s} &= F_1(y_6), \\
\frac{dy_6}{dt_s} &= F_6(v_0, v_1, y_6), \\
\frac{dv_2}{dt_s} &= \epsilon F_2(y_7), \\
\frac{dy_7}{dt_s} &= \epsilon F_7(v_0, v_2, y_7),
\end{align*}
\tag{7}
$$

which describes the dynamics on $S^0$. System (7) is a two-time-scale problem for $\epsilon$ sufficiently small and it gives the slow layer problem in the $\epsilon \to 0$ limit,

$$
\begin{align*}
0 &= F_3(y_8), \\
0 &= F_8(v_0, v_2, v_3, y_8), \\
\frac{dv_0}{dt_s} &= F_0(y_5), \\
\frac{dy_5}{dt_s} &= F_5(v_0, v_1, v_2, v_3, y_5), \\
\frac{dv_1}{dt_s} &= F_1(y_6), \\
\frac{dy_6}{dt_s} &= F_6(v_0, v_1, y_6),
\end{align*}
\tag{8}
$$

where $v_2$ appears as a parameter. A periodic orbit of (7) for $\epsilon = 0.01$ and the bifurcation diagram of (8) as a function of $v_2$ is projected on the $(v_0, v_2)$-plane in S2B Fig.

In the slow layer problem (8), the slow variables $(v_0, y_5, v_1, y_6)$ evolve along fibers defined by $(v_3, y_8, v_0, y_5, v_1, y_6, v_2, y_7) = (v_3, y_8, v_0, y_5, v_1, y_6, v_2, y_7)$, where $(v_2^0, y_7^0)$ are constant and $(v_3, y_8, v_0, y_5, v_1, y_6, v_2^0, y_7^0)$ restricted to $S^0$. The equilibria of (8) defines the superslow manifold $L^0$

$$
\begin{align*}
F_3(y_8) &= 0 \cap F_8(v_0, v_2, v_3, y_8) = 0 \cap \\
L^0 &= \{(v_3, y_8, v_0, y_5, v_1, y_6, v_2, y_7) \in S^0 | \\
& \phantom{=} F_0(y_5) = 0 \cap F_5(v_0, v_1, v_2, v_3, y_5) = 0 \cap \}, \\
F_1(y_6) &= 0 \cap F_6(v_0, v_1, y_6) = 0
\end{align*}
$$

which is a subset of $S^0$. The superslow manifold $L^0$ is reduced to

$$
\begin{align*}
L^0 &= \{(v_3, y_8, v_0, y_5, v_1, y_6, v_2, y_7) \in S^0 | A S(A \tau_a p + C_2 \tau_a A S(C_1 \tau_a v_0) - C_4 \tau_b v_2 - C_7 \tau_g v_3) - v_0 = 0\}, \\
\text{where } v_3 &= G S(C_5 \tau_a v_0 - C_6 \tau_b v_2) \text{ is on } S^0. \text{ The curve } L^0 \text{ perturbs to locally slow invariant}
\end{align*}
$$
manifolds for $\epsilon>0$ along the hyperbolic branches of $L^0$, while the dynamics of near the non-hyperbolic fold points should be investigated using GSPT. Finally, the *superslow reduced problem* obtained by setting $\epsilon\to 0$ in (5) reads

$$ 
\begin{aligned}
\frac{dv_2}{d\tilde{t}_s} &= F_2(y_7), \\
\frac{dy_7}{d\tilde{t}_s} &= F_7(v_0, v_2, y_7).
\end{aligned}
\qquad (9) $$

This algebraic-differential system determines the superslow dynamics restricted to $L^0$, and eventually to $S^0$.

## Strongly perturbed system

Electrical (through direct stimulation) and optical (through optogenetics, using light pulses in genetically modified animals) perturbations can alter action potential firing through modification of the mean membrane potential of the considered neural subpopulation. In the following, we use the word “perturbation” to refer to a strong suprathreshold input. We assumed an additive model for the stimulation effect onto the mean membrane PSP [38]. Thus, the external input $I_{ext}(t)$ is included in the ‘pulse-to-wave’ functions of the NMM in (3), and the system reads:

$$ 
\begin{aligned}
\tau_g \dot{v}_3 &= y_8, \\
\tau_g \dot{y}_8 &= G S(C_5 \tau_a v_0 - C_6 \tau_b v_2 + k_{PV} I_{ext}) - v_3 - 2y_8, \\
\tau_a \dot{v}_0 &= y_5, \\
\tau_a \dot{y}_5 &= A S(A \tau_a p + C_2 \tau_a v_1 - C_4 \tau_b v_2 - C_7 \tau_g v_3 + k_P I_{ext}) - v_0 - 2y_5, \\
\tau_a \dot{v}_1 &= y_6, \\
\tau_a \dot{y}_6 &= A S(C_1 \tau_a v_0 + k_P I_{ext}) - v_1 - 2y_6, \\
\tau_b \dot{v}_2 &= y_7, \\
\tau_b \dot{y}_7 &= B S(C_3 \tau_a v_0 + k_{SOM} I_{ext}) - v_2 - 2y_7.
\end{aligned}
\qquad (10) $$

where $k_i$ with $i = \{P,SOM,PV\}$ is the coupling coefficient between the stimulation and the considered subpopulation, governing the impact of $I_{ext}(t)$ on the subpopulation. Same generalization holds for the slow and superslow reduced problems given by (7) and (9), respectively.

## Clinical data

Clinical data used for the purpose of this study consisted in SEEG signals collected in a patient with drug-resistant focal epilepsy that required invasive EEG exploration. Recordings were performed using intracranial multichannel electrodes (DIXI Medical, 5–18 contacts; length, 2 mm, diameter, 0.8 mm; 1.5 mm apart), 10 of which are used for investigations. The patient-specific position of each electrode results from a clinical decision process which is based on the many assumptions made from non-invasive data (video EEG recording, semiology, neuro-imaging data, clinical examination). Electrodes were implanted according to the stereotactic method of Talairach [39]. SEEG signals were recorded on a Deltamed system on a maximum number of channels equal to 128, and were sampled at 256 Hz and recorded to hard disk (16 bits/sample) using no digital filter. The only filter present in the acquisition procedure was a hardware analog high-pass filter (cut-off frequency of 0.16 Hz) used to remove very slow variations that sometimes contaminate the baseline. In the patient for which data is displayed in the
remainder of the manuscript, a surgical operation was performed 6 month after pre-surgical exploration (cortectomy of the frontal dorsolateral region). The chosen monopolar signal was recorded on one electrode located in the depth of a sulcus located in the dorsolateral frontal cortex. It was chosen because the recorded structure was highly epileptogenic considering the exhibited epileptiform activity (during interictal periods and at the onset of seizures). Histological data revealed the presence of a focal cortical dysplasia (Taylor). After surgery, the patient was seizure free (Engel IA). As a reminder, SEEG is always carried out as part of normal clinical care of patients who give informed consent in the usual way. Patients were informed that their data may be used for research purposes.

## Computational methods

The bifurcation analysis was in done with AUTO-07p [40]. Model equations were implemented in XPPaut [41]. Stochastic differential equations were iterated using Euler-Maruyama method with a step size $dt = 10^{-4}$ second. All simulation files are available from the GitHub database (https://github.com/elifkoksal/NMM_BurstingDynamics).

## Results

### Pre-ictal spiking during interictal to ictal transition

In partial (i.e. focal) epilepsies, the onset of seizures is characterized by the appearance of a rapid discharge, typically in the gamma frequency band ([25, 120] Hz) [27]. This fast onset activity has long been recognized as a hallmark of the epileptogenic zone, and a number of methods have been proposed to make use of this pattern to identify the epileptogenic zone [42,43]. Interestingly, fast onset activity can be preceded by a specific electrophysiological pattern consisting of sustained large amplitude bursts with superimposed faster spikes, which can be observed in various etiologies [44]. Two examples of this pre-ictal pattern, as recorded in two different patients during pre-surgical investigation with depth electrodes, are shown in Fig 1. As depicted, this dynamical regime starts with sporadic bursts, which become periodic to change into a sustained discharge of pre-ictal bursts. In the first recording in Fig 1A, the number of spikes of the bursts gradually decreases during the pre-ictal burst phase, which continues approximately for 14 seconds. In the second example in Fig 1D, pre-ictal bursting continues for approximately 35 seconds. In this example, the amplitude of pre-ictal bursts increases. The pre-ictal burst phases are followed by the fast activity that actually marks the onset of the seizure. Patterns of SEEG signals recorded during interictal-ictal transitions for different seizures of a given patient are typical and reproducible [45,46]. As for the patients whose pre-ictal activity are exemplified in Fig 1, such a transition is persistent with quantitative variations in 1) duration of the pre-ictal bursting period, 2) amplitude of discharges, and 3) number of fast spikes in the bursts.

System (1) represents a physiologically relevant system that includes neuronal subpopulations of excitatory glutamatergic pyramidal neurons, inhibitory SOM+ and PV+ GABAergic interneurons, and makes average PSPs at the level of each subpopulation accessible. This NMM has been extensively explored to establish relationships between model parameters and electrophysiological patterns observed in SEEG recordings [22,47]. For instance, increasing the ratio of the synaptic gain of the excitatory pyramidal cell population and inhibitory SOM+ interneuron population introduces a region in the parameter space where the system can undergo different stages of epileptogenic activity as a function of the synaptic gain of inhibitory SOM+ interneuron population, parameter *B*, and the synaptic gain of inhibitory PV+ interneuron population, parameter *G*.
<figure><img src="image_4.png" /><figcaption>**Fig 1.** SEEG signals recorded in two patients with epilepsy during the interictal to ictal transition and simulated signals. (a) Epileptic seizure recorded in a patient showing the typical pre-ictal spiking pattern with three phases: sporadic spikes, pre-ictal bursts, and fast onset. (b) Zoom into each phase of the actual SEEG signal. The pre-ictal burst type-1 is followed by the pre-ictal burst type-2. (c) Simulated signals corresponding to each phase. (d) Epileptic seizure recorded in a different patient showing the typical pre-ictal spiking pattern with three phases. The region of pre-ictal bursting is zoomed in panel (e).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1008430.g001
The thorough exploration of system (1) led to the identification of three key parameters: *B*, *G*, and the strength of the excitatory synaptic coefficient on the PV+ interneuron subpopulation *C*<sub>5</sub>. Indeed, the tuning of these three parameters enables replicating of the different pre-ictal stages shown in Fig 1A. These results are illustrated by the bifurcation diagram in Fig 2A for the parameter set given in Table 1. As depicted, the decrease of parameter *B* yields a transition from background activity to fast onset activity, though pre-ictal spiking. Fig 2B shows where these activity regions are localized in the (*B*,*G*)-parameter space.

Let us walk through the bifurcation diagram in Fig 2A, starting the equilibrium point at *B* = 0, and increase *B*. The system first undergoes a supercritical Hopf bifurcation (H<sub>1</sub>) at *B*≈1.51, giving stable limit cycles at ≈30 Hz (gamma activity). The amplitude of these solutions increases with *B* and they become unstable *via* a saddle node bifurcation of periodic orbits at *B*≈6.56. The branch of periodic orbits connects to a Hopf bifurcation (H<sub>2</sub>) of subcritical type at *B*≈4.56, while the range of oscillations extends until *B*≈6.56. The branch of equilibrium points undergoes two limit point (LP) bifurcations (LP<sub>1</sub> for *B*≈7.43 and LP<sub>2</sub> for *B*≈4.55) between H<sub>1</sub> and H<sub>2</sub>. Then, we identify a third Hopf bifurcation (H<sub>3</sub>) of supercritical type for *B*≈7.74, where a branch of stable oscillations around ≈6 Hz appears. As *B* increases, this branch connects to stable bursting orbits by passing through several limit folds around *B*≈8.86. At *B*≈10 the system reaches to the maximum number of spikes per burst orbit (11 spikes for this parameter set). Increasing *B* decreases the number of spikes *via* the horizontal up-and-downs in *y*<sub>0</sub> between *B*∈(9.5,22.5). The bursts terminate at *B*≈22.5. The branch holding the unstable equilibrium points forms a Z-shaped curve with two folds at *B*≈21.3 (LP<sub>3</sub>) and *B*≈35.6 (LP<sub>4</sub>), with unstable focus on the upper branch, saddles in the middle and stable nodes on the lower branch after a subcritical Hopf bifurcation (H<sub>4</sub>) which gives unstable limit cycles making a heteroclinic connection with the middle branch. For *B*>35.6, the system only has stable equilibrium points as solutions in the unnumbered white area, which corresponds to physiological background activity.

Continuation of the LP and Hopf bifurcations marked on Fig 2A in the (*B*,*G*)-space is shown in Fig 2B. It can be seen that the locations of LP<sub>3</sub>, LP<sub>4</sub>, H<sub>3</sub>, H<sub>4</sub> points do not depend on *G*, whereas the locations of LP<sub>1</sub>, LP<sub>2</sub>, H<sub>1</sub> and H<sub>2</sub>, which are related to the oscillations at ≈30 Hz, do. The fast onset region does not exist for small values of *B* if *G*<5, for which the system only has equilibrium points for *B*<*B*<sub>H3</sub>. As *G* increases, the system undergoes Hopf bifurcations (H<sub>1</sub> and H<sub>2</sub>) that introduce gamma-band fast oscillations (for *G*>5). The amplitude of these fast oscillations are smaller closer to H<sub>1</sub>, and this difference is more visible in simulated LFPs. Under stochastic inputs, the system can yield fast oscillations for *B*∈(0, *B*<sub>H1</sub>), and a mixture of fast and slow oscillations of ≈6 Hz for *B*∈(*B*<sub>H2</sub>, *B*<sub>H3</sub>). Furthermore, *G* controls the amplitude of the spikes of bursting solutions, which increases with *G*.

Assume that system (1) is initially in the background activity mode, which corresponds to the white region in Fig 2A for $B>35.6$, where unique stable equilibrium points on the bifurcation curve is observed. In the blue region between the two folds $LP_3$ and $LP_4$, the bifurcation curve takes a Z-form with stable nodes on the lower branch, unstable nodes in the middle and saddle-nodes on the upper branch. For $B$ values in this blue region, system (1) under a stochastic input $p(t)$ undergoes sporadic bursts, with an increasing probability as $B$ approaches to the left fold. As $B$ decreases, the system enters into the bursting region (orange region). Note that further increasing $B$ would increase the number of spikes. However, in the recordings we see that the number of spikes decreases in the course of the pre-ictal bursting regime as the system approached the low voltage fast onset (LVFO), transition from type-1 to type-2 bursting. This change is very subtle to be reproduced in the model because, as detailed in the Sec. *Burst analysis*, the number of spikes depends on the presynaptic potential on PV+ interneurons: the lower it is, the more spikes within the burst are obtained. Thus, the number of spikes increases when
<figure><img src="image_5.png" /><figcaption>**Fig 2. Bifurcation diagrams of the system (1) showing the different pre-ictal stages numbered 1–3 in Fig 1. (a)**
Amplitude of the PSP of the pyramidal cell subpopulation is plotted as a function of the synaptic gain of the SOM+ interneuron subpopulation *B*. The other system parameters are given in Table 1 with *G* = 35. Bold and dashed lines correspond to stable and unstable solutions, respectively. Region 1 (blue) corresponds to sporadic bursts, region 2 (orange) to sustained bursts, and region 3 (purple) to low voltage fast onset activity. The system yields large amplitude ≈30 Hz oscillations in the unnumbered green shaded region. The unnumbered yellow shaded area corresponds to high amplitude stable equilibrium points, and white corresponds to physiological background activity. The arrow shows the route from background to low voltage fast onset activity in the parameter space. In order to preserve the readability of the diagram, the bifurcations along the branches of periodic solutions are not shown. (b) Co-dimension 2 diagram of the Hopf (H) and limit point (LP) bifurcations marked on panel (a) in the parameter space of *B* and *G* (synaptic gain of the PV+ interneuron subpopulation). The LP₁ and LP₂ points merge on a cusp (CP) bifurcation, and the H₁ and H₂ merge on a zero-Hopf (ZH) bifurcation.</figcaption></figure>

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

the inhibitory input decreases, or when the excitation onto PV+ interneurons increases. At this stage, transition from the type-1 bursting to type-2 bursting is obtained by keeping *B* constant, but decreasing the excitatory post-synaptic potential (EPSP) on PV+ by decreasing progressively *C*₅ to 300 to reduce the number of spikes; and increase *G* to increase spikes amplitude. Under these variations, the bifurcation diagram in Fig 2A remains qualitatively the same, the most important quantitative change being the location of the Hopf bifurcation point H₁ related to the LVFO. As shown in Fig 2B, increasing *G* moves the H₁ towards right along the *B*-axis, and initiates the LVFO for higher values of *B*≈*B*<sub>H1</sub>. Under stochastic input, one can observe LVFO for *B*≤*B*<sub>H1</sub>, as well (purple regions in Fig 2). The large amplitude fast oscillations that are obtained for higher *B* values in (*B*<sub>H1</sub>, *B*<sub>H2</sub>) can also be observed in SEEG recordings during the transition to epileptic seizures.

## Bursting analysis

We investigate the bursting dynamics of system (1) using system (3), which is a kind of nondimensionalized version of (1) but expressed in an explicit slow-fast formulation. Fig 3A shows a bursting solution of (3) in the (v₀, v₂, v₃)-space, the critical manifold S⁰ and the superslow manifold L⁰ (see Sec. *Slow-fast formulation of the model* for definition). The critical manifold S⁰ is normally hyperbolic (not folded) and stable, and stretches between almost horizontal surfaces (lower and upper) with an almost vertical plane. The superslow manifold L⁰ has branches both on the lower horizontal surface and vertical surface of S⁰. While the part of L⁰ on the vertical surface of S⁰ is stable, the part on the lower horizontal surface of S⁰ is divided into stable and unstable sections at two fold points LP₁ and LP₂. The curve L⁰ is stable along the branch that is almost parallel to the v₂-axis, unstable along the branch between LP₁ and LP₂, and then
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Bursting orbit of system (3). (a) Solution of (3) (blue orbit) and $L^0$ (red curve) on the critical surface $S^0$(green surface) projected on the $(v_0, v_2, v_3)$-space. Single-headed, double-headed and triple-headed arrows indicate the flow direction during superslow, slow and fast time-scales, respectively. LP denotes limit point bifurcation. The $L^0$ curve changes stability at two limit points, $LP_1$ and $LP_2$ (red dots). The middle branch of the $L^0$ curve between these limit points is unstable (dashed). (b) Time course of the variables $(v_3, v_0, v_2)$ of the orbit plotted in panel (a). (c) Solution of (3) projected on the bifurcation diagram (black curve) of (4) for $\varepsilon = 0$ where $v_2$ is threaded as a parameter. Arrows show the direction of the flow with respective time-scales. Bold and dashed lines correspond to stable and unstable solutions, respectively. H donates a Hopf bifurcation, LP a limit point bifurcation. The equilibrium points along the black Z-shaped curve are unstable on the middle branch of the curve, between $LP_1$ at $v_2^{LP1} = 4.778$ and $LP_2$ at $v_2^{LP2} = 20.66$ (black dots), and on the upper branch between $H_1$ at $v_2^{H1} = 0.27$ and $H_2$ at $v_2^{H2} = 14.27$ (green dots). The amplitude of the stable limit cycles is bounded by the green continuous curves connecting the $H_1$ and $H_2$ in the $\varepsilon = 0$ limit.</figcaption></figure>

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

becomes stable again. The stable and unstable branches of $L^0$ are normally hyperbolic, whereas the fold points $LP_1$ and $LP_2$ are not.

The critical manifold $S^0$ and superslow manifold $L^0$ perturb for small enough values of time-scale parameters, hence the three-time-scale dynamics of (3) approximate to $S^0$ and $L^0$. During the superslow time-scale under (9), the bursting orbit follows the stable branch of $L^0$ almost parallel to the $v_2$-axis. Near the fold point $LP_1$, the trajectory bends in the $v_3$-direction along the vertical plane of $S^0$ and enters into the spiking regime, which runs in fast time-scale. The spiking terminates close to $LP_2$ and the trajectory jumps back to the stable branch of $L^0$ almost parallel to the $v_2$-axis in slow time-scale under (8). Fig 3B shows the time series in $\tilde{t}$ of the orbit in Fig 3A.

For a better understanding of the bursting dynamics, we consider system (4) at $\varepsilon = 0$ for which the variables of the slowest subsystem $(v_2, y_7)$ act as parameters of the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem. Since only $v_2$ appears in the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem, its dynamics depend on $v_2$. In Fig 3C, the bursting orbit in Fig 3A is superimposed on the bifurcation diagram of the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem in (4) at $\varepsilon = 0$ as a function of $v_2$. Although the fastest variables of (4) are $(v_3, y_8)$, we chose $v_0$ vs $v_2$ for a clearer visualization (the same trajectory and the bifurcation diagram are given on the $(v_3, v_2)$-plane Fig 4). We see that the corresponding system poses a Z-shaped bifurcation diagram as a function of $v_2$ with two folds, $v_2^{LP1}$ and $v_2^{LP2}$. The equilibrium points are stable on the lower branch of the Z-shaped curve for $v_2 > v_2^{LP1}$, unstable along the middle branch between $v_2^{LP1}$ and $v_2^{LP2}$. The upper branch has two
<figure><img src="image_7.png" /><figcaption>**Fig 4.** Spike number as a function of $C_5$. (a) Solution of (4) with 3 spikes for $C_5 = 500$ projected on the bifurcation diagram of the fast system (5) as a function of $v_2$. Arrows indicate the direction of the flow. (b) Solution of (4) with 1 spike for $C_5 = 300$ projected on the bifurcation diagram of the fast system (5) as a function of $v_2$. Arrows indicate the direction of the flow. (c) Co-dimension 2 diagrams of the Hopf (H) points (green) and the limit points (LP, black) in the $(v_2, C_5)$ parameter space marked on the left and middle panels. As $C_5$ decreases, $H_2$ moves leftwards and eventually the spike number decreases. For $C_5 = 139$, $H_2$ and $LP_1$ are aligned at $v_2 = 4.778$. A further decrease in $C_5$ places $H_2$ on the left of $LP_1$ and leaves no chances for a bursting solution.</figcaption></figure>

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

supercritical Hopf bifurcations, at $v_2^{H1}$ and $v_2^{H2}$, with stable limit cycles in between. Along the upper branch, equilibrium points are stable for $v_2 < v_2^{H1}$ and $v_2^{H2} < v_2 < v_2^{LP2}$. The bursting behavior resulting from this bifurcation structure in the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem is classified as ‘fold/Hopf bursting’ by Izhikevich [32] due to the presence of a ‘fold/Hopf hysteresis in the bifurcation diagram.

System (4) may undergo through these bifurcations in a repetitive manner for $\varepsilon \neq 0$, which results eventually in the bursting solutions for small enough values of $\varepsilon$. As the arrows on Fig 3C and the traces on Fig 3B demonstrate, the trajectory follows the lower stable branch during the quiescence phase of the bursting, which terminates near $v_2 \approx v_2^{LP1}$. Then, it jumps to the region of the stable limit cycles on the upper branch, which initiates the active phase of the bursting. The spiking frequency during the active phase is faster at the beginning than the end due to the fact that the Hopf bifurcation at $v_2^{H1}$ gives limit cycles with $\approx 30$ Hz frequency whereas the Hopf bifurcation at $v_2^{H2}$ gives limit cycles with $\approx 10$ Hz. The spiking terminates at $v_2 \approx v_2^{H2}$, but the active phase continues until the trajectory jumps back to the stable lower branch at $v_2 \approx v_2^{LP2}$. We underline that as $\varepsilon \to 0$, the bursting orbit attaches more and more the bifurcation diagram obtained for $\varepsilon = 0$ (see S2 Fig for an example).

The main difference between the type-1 and type-2 bursting is the number of spikes during the active phase of bursting. In the model, the variations in the number of spikes can be met by changing the excitation level on the PV+ interneurons: as aforementioned, the number of spikes increases with the amount of excitation received by PV+ interneurons. This can be achieved by either decreasing inhibition or by increasing excitation. For instance, decreasing $B$ in region-2 in Fig 2 increases the number of spikes. In (4) at $\varepsilon = 0$, the excitation on PV+ depends on two synaptic coupling coefficients, $C_5$ and $C_6$. The effect of $C_6$ will be similar to the one of $B$, since they both scale the PSP of SOM+ interneurons given by the variable $v_2$ in (4). Below, the role of the excitatory synapses in the $(v_3, y_8, v_0, y_5, v_1, y_6)$-subsystem by changing $C_5$ is investigated.
As displayed by Figs 3C, 4A and 4B, the spikes are bounded by $LP_1$ and $H_2$ in the bifurcation diagram of the ($v_3, y_8, v_0, y_5, v_1, y_6$)-subsystem as a function of $v_2$. The distance between $LP_1$ and $H_2$ in $v_2$ affects the number of spikes; the further they are, the more spikes the burst has. In Fig 4C, $LP$ and $Hopf$ bifurcations are continued in the parameter space of ($v_2, C_5$). While the $LP_1$ and $LP_2$ lie along almost vertical lines, the $Hopf$ bifurcation points form a V-shaped curve along which the left arm locates the $H_1$ points and the right arm the $H_2$ points. The distance between $H_2$ and $LP_1$ increases with $C_5$, hence, the spike number. At $C_5 = 139$, $H_2$ and $LP_1$ are aligned. Further decrease in $C_5$ places $H_2$ on the left of $LP_1$ and leaves no chances for a bursting solution. The system yields only relaxation type of oscillations for $C_5<139$.

Overall, the aforementioned analysis shows that pre-ictal bursting runs in three-time-scales. The system sustains the bursting regime for a certain range of parameter $B$ denoting SOM+ synaptic gain. The complex pre-ictal bursting pattern can be accurately adjusted by tuning parameters $G$, which controls the PV+ synaptic gain, and the connectivity coefficient $C_5$, which controls PV+ excitability. In particular, the number of spikes and their amplitude can be adjusted by tuning $C_5$ and $G$, respectively.

## Strong perturbation analysis

The mean membrane potential of a neuronal subpopulation can be altered by electrical and optical stimulations. Under the assumption of an additive model for the stimulation effect, an external input can be included in the ‘pulse-to-wave’ functions, $S(v)$, of the NMM after being scaled by the subpopulation specific impact coefficients (see Sec. *Strongly perturbed system*).

A pulse input (biphasic or monophasic) changes the PSP of the perturbed subpopulation by shifting it above its base level $S(0)$. We assume that a neural mass block, given by

$$ \dot{A}y = W/\tau_w S(I_{ext}(t)) - 2/\tau_w \dot{y} - 1/\tau_w^2 y, $$ receives biphasic stimulation. The PSP of the neural mass block increases during the anodal pulse (positive perturbation or depolarization of the membrane potential), but decreases (discharges) during the cathodic pulse (negative perturbation or hyperpolarization of the membrane potential) and between the inter-pulse intervals of the biphasic input. Depending on the pulse width, pulse amplitude, and mostly on the synaptic time constant of the neural mass block, this shift may be sustained or not. For instance, the discharge takes longer in a neural mass block with slow synaptic kinetics than the one with fast synaptic kinetics. If the pulse frequency is sufficiently high to stimulate the neural mass before it completely resumes to its base level, then the PSP of the neural mass can oscillate above the base level. As visualized in S3 Fig, the same perturbation shifts the PSP of a neural mass with slow synaptic kinetics, while the neural mass with fast synaptic kinetics decays to its base level during the inter-pulse intervals of the stimulation. Increasing the stimulation frequency can keep the PSP of the neural mass with fast kinetics above the base level, and therefore the firing rate and PSP of the fast neural mass increase with the stimulation frequency.

The bursting solution is driven by the slow oscillations in system (3) (see Sec. *Bursting Analysis*). The slow dynamics of (3) (subsystems representing the pyramidal cell and SOM+ inter-neuron subpopulations) can be approximated by (7), which preserves the burst-driver slow oscillations behavior for the same parameter values yielding bursting oscillations in (3) (see Sec. *Slow-fast formulation* and S2 Fig). Thus, it is sufficient to investigate the response of (7) under perturbation to understand the impact of the perturbation on the burst-driver slow dynamics. The most common signal delivered to brain tissue in the field of deep brain stimulation (DBS) is bi-phasic pulses with balanced anodic/cathodic phases of brief durations (approximately 100 µs). Below, the impact of anodic and cathodic constant external inputs is considered without taking into account their duration, to simply understand how they alternate the phase space of system (7). For this purpose, a constant input ($I_{ext} = 1$) scaled with the
impact coefficients $k_p$ and $k_{SOM}$ is applied to (7) by following the formulism given by (10) (see Sec. *Strongly perturbed system*).

In Fig 5, subpopulations representing pyramidal cells and inhibitory SOM+ interneurons are perturbed. The left panels of Fig 5 show the $y_5$-nullsurface $\Theta$, $v_2$-nullsurface $\Sigma$ and the superslow manifold $L^0$ projected on the ($y_5$, $v_2$, $v_0$)-space. The solution of (7) is visible on the left panels, and the solution of (3) for the same parameters is given on the right panels of Fig 5. Fig 5A shows the case where only the SOM+ interneuron subpopulation described by the ($v_2$, $y_7$)-subsystem in (7) is subject to the constant external input ($k_p = 0$). In the absence of any perturbation ($k_{SOM} = 0$), $\Theta$ and $\Sigma$ intersect for $v_2 > 20$. System (7) has a limit cycle which flows on $\Theta$ and (3) a burst orbit (black solutions Fig 5A and 5A1, respectively). The quiescence phase of the burst corresponds to the slow passage following $L^0$ where $v_0 \approx 0$, and the active phase correspond to the trajectory on the upper sheet of $\Theta$.

A key point in terms of controlling bursting activity through direct stimulation is that an input leading to a bifurcation from the stable limit cycle to an equilibrium point can prevent the system from bursting by keeping the system in the silent phase. This can be achieved by an input that ensures an intersection between $\Sigma$ hyperplane and the lower branch of $L^0$. Indeed, for $k_{SOM} = 1$, (7) possesses a stable equilibrium point near the left fold of $L^0$ which traps the trajectory (green dot in Fig 5A). For the same input, the bursting in (3) is aborted (green solution in Fig 5A1). On the other hand, a negative constant input ($k_{SOM} = -1$), moves $\Sigma$ away from the left fold of $L^0$. Being $\Sigma$ closer to the upper branch of $L^0$ prolongs the active phase of the burst and increases the number of spikes, as seen in Fig 5A2. These observations indicate that increasing the excitation on SOM+ interneurons can abort bursting.

In Fig 5B, only the subsystem representing the pyramidal cells receives the perturbation ($k_{SOM} = 0$). The input on the pyramidal cell subpopulation acts on $\Theta$. While positive constant input ($k_p = 1$) increases the distance between the lower fold of $L^0$ and $\Sigma$, negative constant input ($k_p = -1$) decreases this distance. Both systems (7) and (3) preserve the oscillatory behavior for these values of $k_p$, yet, the oscillation frequency decreases for $k_p = -1$ due to the decreased distance between the lower fold of $L^0$ and $\Sigma$. Thus, hyperpolarization of pyramidal cells by increasing inhibition on them can abort bursting.

As aforementioned above, pulsed stimulation increases the firing rate of a neuronal population. However, a stimulation applied to one specific region might not affect all neural populations in the same manner. This can be due to the relative position of electrodes with respect to neurons, cell specific firing thresholds, or synchronization level within neural subpopulations. However, such features can bring certain advantages in aborting bursting. Fig 5C shows the response of the system when both subpopulations of pyramidal cells and SOM+ interneurons are perturbed, the oscillatory behavior in systems (7) and (3) continues under the same positive constant input ($k_p = k_{SOM} = 1$). With such input, the number of spikes during the active phase is decreased (Fig 5C1). If the subpopulation of SOM+ interneurons is perturbed more strongly than the subpopulation of pyramidal cells ($k_p = 1, k_{SOM} = 2$), the system can bifurcate to the resting state (Fig 5C2).

Although the reduced system (7) does not include the fast dynamics of the PV+ interneurons, the effect of the perturbation on the subpopulation of PV+ interneurons can be understood geometrically. First, let us notice that increasing the inhibition on SOM+ interneurons encourages spiking (Fig 5A and 5A2), while increasing the excitation on SOM+ aborts bursting (Fig 5A and 5A1). Perturbing (stimulating) the subpopulation of PV+ interneurons increases the PSP from PV+ interneurons to pyramidal cells, reduces the PSP from pyramidal onto SOM+, and in turn favors bursting. Another way to illustrate the impact of perturbing the PV+ on bursting is to examine the diagram in Fig 4. Anodic pulses can shorten the quiescent phase in the ‘fold/Hopf’ hysteresis loop by kicking the trajectory to the region of stable
<figure><img src="image_8.png" /><figcaption>**Fig 5.** Geometrical analyses of a constant input. Constant input is applied to SOM+ interneurons (a), to pyramidal cells (b) and to both SOM+ interneurons and pyramidal cells (c). Left panels show the projection of the nullsurfaces, critical slow manifold and the orbit of the reduced model (7). Right panels show the LFP signal of the full system (3) subject to the constant inputs analyzed on the left. All parameters are as given in Table 1, except $B = 15$. (a) The $y_5$-nullsurface $\Theta$ (blue surface for $k_p = 0$), and $y_7$-nullsurface $\Sigma$ (red surface for $k_{SOM} = -1$, black surface for $k_{SOM} = 0$, green surface for $k_{SOM} = 1$) are projected on the ($v_2, y_5, v_0$)-space. The blue curve $L^0$ (stable on the bold, unstable on the dashed) is on the intersection between $\Theta$ and the $\{y_5 = 0\}$-hyperplane. The black and red orbits are the solutions of the system for $k_{SOM} = 0$ and $k_{SOM} = -1$, respectively. For $k_{SOM} = 1$, the solution approaches to the green stable equilibrium point on the intersection between $\Sigma(k_{SOM} = 1)$ and $L^0$. Panel (a1) shows the time series for $k_{SOM} = \{0, 1\}$, and panel (a2)</figcaption></figure>
for $k_{SOM} = -1$. (b) The $y_5$-nullsurface $\Theta$ (red surface for $k_P = 1$, green surface for $k_P = -1$), and $y_7$-nullsurface $\Sigma$ (black surface for $k_{SOM} = 0$) are projected on the $(v_2, y_5, v_0)$-space. The red curve $L^0$ (stable on the bold, unstable on the dashed) is on the intersection between $\Theta(k_P = 1)$ and the $\{y_5 = 0\}$-hyperplane. The green curve $L^0$ (stable on the bold, unstable on the dashed) is on the intersection between $\Theta(k_P = -1)$ and the $y_5 = 0$ hyperplane. The green and red orbits are the solutions of the system for $k_P = 1$ and $k_P = -1$, respectively. Panel (b1) shows time series for $k_P = 1$, and panel (b2) for $k_P = -1$. (c) The $y_5$-nullsurface $\Theta$ (red surface for $k_P = 1$) and $y_7$-nullsurface $\Sigma$ (green surface for $k_{SOM} = 1$, blue surface for $k_{SOM} = 2$) are projected on the $(v_2, y_5, v_0)$-space. The red curve $L^0$ (stable on the bold, unstable on the dashed) is on the intersection between $\Theta(k_P = 1)$ and the $y_5 = 0$ hyperplane. The green curve $L^0$ (stable on the bold, unstable on the dashed) is on the intersection between $\Theta(k_P = -1)$ and the $\{y_5 = 0\}$-hyperplane. The green orbit is the solution of the system for $(k_P, k_{SOM}) = (1,1)$. For $(k_P, k_{SOM}) = (1,2)$ the solution approaches to the cyan stable equilibrium point on the intersection between $\Sigma(k_{SOM} = 2)$ and $L^0$. Panel (c1) shows time series for $(k_P, k_{SOM}) = (1,1)$, and panel (c2) for $(k_P, k_{SOM}) = (1,2).$

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

limit cycles between the two Hopf bifurcations $H_1$ and $H_2$. Hence, such pulses applied periodically can increase the bursting frequency by shortening the quiescent phase. On the other hand, cathodic pulses can lengthen the quiescent phase by hooking the trajectory near the down state of the hysteresis loop. One can also think of inhibitory effect of PV+ subpopulation on the pyramidal cell population which could abort bursting. Indeed, it would be possible if an input was capable of depolarizing the PV+ subpopulation during the quiescence phase of the bursting where it is inhibited by the SOM+ subpopulation. To overcome the inhibition from the SOM+ subpopulation to the PV+ subpopulation, such an input should be very high in amplitude (order of 20 if a constant input on the PV+ subpopulation is considered), which could induce undesired discharges.

Overall, this geometric perturbation analysis helps to clarify the role of hyperpolarizing and depolarizing inputs on ongoing bursting activity. In particular, depolarization of the subpopulation of SOM+ interneurons or hyperpolarization of the subpopulation of pyramidal cells can abort bursting by keeping the sum of PSPs at low levels. Depolarization of the subpopulation of PV+ interneurons contributes to bursting.

## Stimulation applied during the pre-ictal burst regime

The analysis in Sec. *Perturbation Analyses* has shown that a positive constant input applied to the subpopulation of SOM+ interneurons can bifurcate the limit cycle (oscillating epilepsy-like activity) to an equilibrium point (background activity), while a positive constant input on the subpopulations of pyramidal cells and PV+ interneurons preserve bursting and high frequency oscillations (Fig 5). Hence, an appropriate strategy for pre-ictal bursting abortion consists in the excitation of the SOM+ interneuron subpopulation.

In this section, the results obtained from the mathematical analysis are translated into an *in silico* set-up mimicking experimental conditions. Typically, charge-balanced bi-phasic pulses (pulse width = 0.5 ms and total duration 1 ms) with an arbitrary unit (arb. unit) amplitude are applied during the pre-ictal bursting/spiking regime in the presence of a stochastic input. In order to test our predictions on the role of different neural populations, only SOM+ interneurons are perturbed in Fig 6A ($k_{SOM} = 1$, $k_P = 0$, $k_{PV} = 0$), whereas in Fig 6B all neural subpopulations are perturbed homogenously (coupling coefficients $k_{SOM} = k_P = k_{PV} = 1$).

Results indicate that pre-ictal bursts frequency decreases when the stimulation is switched on at, typically at the instant $t = 5s$ in both cases. The bursting regime can be aborted if the stimulation frequency and amplitude are sufficiently high. The minimum values of the stimulation frequency and amplitude to abort bursting depend on which neuronal subpopulation receives the stimulation. When only the SOM+ interneuron subpopulation is stimulated, the minimum stimulation frequency and amplitude required to abort bursting are lower than the case where all neural subpopulations are stimulated homogenously. As exemplified in Fig 6C,
<figure><img src="image_9.png" /><figcaption>**Fig 6.** System (1) under stimulation. Biphasic stimulation with a 0.5 ms pulse width (total pulse duration is 1 ms) is applied to the system in type-1 bursting. Panels (a) and (b) show the energy map of the simulated LFP signal that is lower in the blue region than the yellow region (see the color bar on the right). The energy of the LFP signal was computed by using $E = \sum_{n=1}^{n=N} |x(n)|^2$, where *n* stands for the index of and *N* for the size of the discrete signal *x(n)*. (a) Only the SOM+ interneurons receive the biphasic perturbation. (b) The pyramidal cell, SOM+ interneurons and PV+ interneurons receive the same biphasic perturbation. Panels (c), (d) and (e) show the time course of the marked stimulation on panels (a) and (b). (c) 15 Hz biphasic stimulation with 10 arb. unit amplitude is applied to the SOM+ interneurons ($k_{SOM} = 1, k_p = k_{PV} = 0$). (d) 15 Hz biphasic stimulation with 10 arb. unit amplitude is homogenously applied to all subpopulations ($k_{SOM} = k_p = k_{PV} = 1$). (e) 25 Hz biphasic stimulation with 10 arb. unit amplitude is homogenously applied to all subpopulations ($k_{SOM} = k_p = k_{PV} = 1$).</figcaption></figure>

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

bursting is suppressed at $f = 15$ Hz for an amplitude of 10 arb. unit when only SOM+ interneurons are stimulated. While the same stimulation can considerably decreases the frequency of bursting events (Fig 6D) when all subpopulations are impacted, the stimulation frequency
should be increased to 25 Hz for a complete bursting suppression (Fig 6E). For a sufficiently high stimulation frequency for which bursting is aborted, increasing the stimulation amplitude can depolarize the PV+ subpopulation. However, such an input induces high frequency components to the burstless response observed in the LFP, as opposed to the low amplitude response which is closer to the physiological background activity (data now shown). While the purpose of the stimulation is not only aborting bursting but also keeping the system as close as possible to background activity, low amplitude stimulations result in more beneficial outputs for clinical applications.

The difference between type-1 and type-2 bursting is the number of spikes during the active phase that is related to the excitatory input onto PV+ interneurons. In particular, the EPSP is larger in the former case. Despite this difference, the bursting mechanisms in both types are the same; i.e. slow oscillations in the SOM+ interneurons drive sequentially and periodically the same type of bifurcations in the subsystem of pyramidal cells and PV+ interneurons. Hence, the strategy for aborting bursting relying on aborting oscillations in the SOM+ subsystem does not depend on the bursting type. The estimations on the stimulation parameters (in terms of frequency and amplitude) given in Fig 6C, which are for type-1 bursting, are capable of aborting type-2 bursting and sporadic bursting, as well, because both of the regimes are less excited than type-1 bursting.

Discussion

Epilepsy is a dynamic and complex disease running on different time-scales [48–50]. Epileptic activity is characterized by long interictal periods (outside seizures), during which the brain behaves mostly as a normal brain, then marked by brief ictal episodes (seizures). The seizure onset, i.e. the transition from interictal to ictal activity, has a wide repertoire in human focal epilepsies [44,51]. In this study, we focused on a specific electrophysiological pattern generally referred to as “pre-ictal spikes” or “pre-ictal discharges”, which has been particularly described in mesial temporal lobe seizures [12–14] but that may also be observed as a seizure onset pattern in neocortical seizures from various origins [44,52]. This complex pattern is signed by large amplitude fast spikes followed by a slow discharge, thus holding the properties of a bursting and is called “pre-ictal bursting” in this paper.

We successfully reproduced the complex pre-ictal bursting pattern in a NMM featuring three subsets of neurons (subpopulations of pyramidal neurons, SOM+ and PV+ interneurons) in [22]. The slow-fast formulation of the model unveiled its three-time-scale structure and the following analysis detailed the mechanisms responsible for the pre-ictal bursting. In particular, the bursting process in the model arose from a high level of excitation among pyramidal neurons as well as onto the PV+ interneuron subpopulation. In the bursting regime, the slow oscillations mediated by the SOM+ interneurons are the drivers of bursting solutions, and the number of spikes during an active phase of a burst depends on the level of excitation on the PV+ interneurons. Ultimately, we showed that a perturbation that was able to stop the slow oscillations in the SOM+ interneuron subpopulation would be sufficient to stop pre-ictal bursting activity.

These model predictions corroborate some experimental findings. Indeed, *in vitro* data from human specimen suggested that a glutamate-dependent cellular and/or synaptic plasticity process underlies the occurrence of pre-ictal discharges during the transition to seizure. Pre-ictal discharges would initiate changes in glutamatergic and GABAergic signaling in groups of neurons larger than those involved in interictal discharges. Repeated discharges would result from a dynamic process that ultimately leads to ictal events [53]. Along the same line, as extensively reviewed in [54], both excitatory and inhibitory networks are involved in
epileptogenesis and seizure generation. In particular, GABAergic-mediated mechanisms con-
tribute to synchronizing neuronal networks in epileptic brain structures. Notably, interneuro-
nal activity is enhanced and synchronized during sustained epileptic spikes [55,56].

This viewpoint is particularly interesting if the role of the GABAergic system in the suppres-
sion of epileptiform pre-ictal activity is considered when direct brain stimulation applied dur-
ing the interictal period. For instance, optogenetic stimulation of the CA3 region of
hippocampus leads to considerable reduction of seizures in the CA3 region by entrainment of
GABAergic interneurons targeting GABA<sub>A</sub> receptors [57,58]. Low-frequency stimulation of
fiber tracts during the inter-ictal period has also been shown to reduce seizures through activa-
tion of the GABA<sub>B</sub> signaling in animal models of temporal lobe epilepsy activity [59–61], as
well as with the application of an electrical field [62]. The success of low-frequency stimulation
of fiber tracts in focal cortical seizures has also been linked to GABAergic effects [63,64].

Our results are in line with the above reported data, and indicate that an abortive stimula-
tion of the epileptic activity during the pre-ictal bursting regime should primarily target the
GABAergic system (mostly on interneurons with slow synaptic kinetics). Stimulating the
GABAergic system yielded more pronounced effect as compared with the stimulation pyrami-
dal neurons. The stimulation frequency required to change the PSPs of neural subpopulations
was directly linked with their kinetics: the slower they are, the lower stimulation frequency
needs to be. At this point, SOM+ interneurons were impacted more than other subpopula-
tions, since SOM+ interneurons have the slowest synaptic kinetics among the considered neu-
ronal types in the model. Besides, the model structure privileges the SOM+ subpopulation to
be depolarized by an external stimulation since the SOM+ subpopulation does not receive any
IPSP from other subpopulations. Increasing the stimulation frequency would depolarize the
SOM+ subpopulation further, thus reinforcing the global inhibition generated by the SOM+
subpopulation. In addition, it has been estimated that a single GABAergic cell may affect more
than a thousand pyramidal cells [65,66], which may explain how the activation of GABAergic
neurons may become predominant and exert powerful anti-epileptic effects.

Another prediction of this study is the contributing role of PV+ interneurons stimulation on pre-ictal bursting. More specifically, depolarizing the subpopulation of PV+ interneurons contributes to bursting by increasing the number of spikes during the active phase. Also, as it was discussed above, anodal pulses on PV+ interneurons can prompt the active phase and increase the frequency of pre-ictal bursts. Such observation is in agreement with a previous study by Assaf and Schiller [67], in which optogenetic activation of PV+ interneurons in the ictal regime had an anti-epileptic effect, but a pro-epileptic effect when they were activated in the inter-ictal regime. More recently, it was discussed that paradoxical effects of PV+ activation shown in [68] could be related to the timing of the neurostimulation [69]. Therefore, our results support that a precise, on-demand (closed-loop) stimulation system is required to deliver stimulation at an optimal timing, and avoid the promotion of epileptiform activity. The properties of epileptic discharges during the ictal phase are different than the ones of the pre-ictal phase considered in this study. Thus, seizure abortion during the ictal phase may require different stimulation protocols. Investigation of suitable targets and signal waveforms can be considered as a possible avenue for future work.

DBS for epileptic patients is an ongoing research topic, and unfortunately, the lack of ran-
domized control trials comparing different stimulation protocols hampers obtaining definite
results on optimal stimulation protocols [8,70]. Low-frequency electrical and optical stimula-
tion (< 5Hz) applied during interictal phases has been shown to reduce the frequency of inter-
ictal spikes and seizure initiation in animal and human studies [58,71]. High-frequency
electrical stimulation (>100 Hz) applied during ictal phases has also been shown to terminate
seizures [72–74]. Here, we considered the pre-ictal phase, which is between the interictal and
ictal phases. We showed that stimulation with a frequency greater than 20 Hz can abort pre-
ictal oscillations and keep the system close to background activity by depolarizing the subpop-
ulation of SOM+ interneurons. From our results, the suggested frequency range lies between
the ranges of low- and high-frequency stimulations and beyond. This can be due to the fact
that the considered epileptogenic phase (pre-ictal) is “in-between” the phases where low-fre-
quency (interictal phase) and high-frequency (ictal phase) stimulations are successful. Further-
more, a single-lumped NMM considers the local synchronized activity at the mesoscopic scale,
whereas epilepsy is a large-scale network disease where heterogeneity and spatial dynamics
can be crucial for epileptic activity. Therefore, spatio-temporal effects of ongoing neural
dynamics, synaptic plasticity or differences in the activation functions or neural subpopula-
tions could not be studied in our model. Nevertheless, our study explains the complex locally
observed pre-ictal patterns. Furthermore, our results suggest an alternative stimulation proto-
col in terms of frequency and timing of stimulation delivery.

Overall, our study adds to the recent literature on computational studies of seizure abor-
tion. In the context of absence seizures, Taylor et al. [75] have proposed application of single
pulse stimulation in a model known to have bistable properties. While close in the level of
description to [75], our model differs in two aspects. First, our model of focal seizures is differ-
ent than of absence seizure due to the difference between underlying neural circuits. Second,
our mathematical analysis is based on the slow-fast features of the dynamical system. There-
fore, transitions from epileptic discharges to normal background activity are obtained by stim-
ulations which differ from a single pulse not only in terms of duration, but also in terms of
waveform.

It has long been reported that pre-ictal spiking/bursting is an emerging feature of the inter-
ictal to ictal transition and is specific to epileptogenic regions. From a mathematical viewpoint,
both spiking (a single bump followed by a quiescent phase) and bursting (a sequence of spikes
(bumps) followed by a quiescent phase) oscillations in a neural context result from the interac-
tion between the slow and fast variables of a multiple time-scale system. While the type of the
oscillation depends on the bifurcation structure of the fast subsystem, it is always the slow sub-
system that drives the recurrent transitions between the quiescence and active (spiking) phases
[32], here the subpopulation of SOM+ interneurons. Since the essence of spiking/bursting is
the same in general sense, stimulation protocols mainly affecting slow oscillations during the
pre-ictal phase would abort pre-ictal spiking/bursting activity. In other words, the burst-abor-
tion strategy presented in this paper would also be appropriate to abort spiking. Yet, it is essen-
tial to identify the neuronal subpopulations of the brain region under consideration, the
connections between these neuronal subpopulations and their roles in such slow-fast regimes
to optimize the stimulation frequency, since as shown in this manuscript, subpopulations with
slower kinetics are more responsive to pulsed stimulations. For instance, a pre-ictal spiking
regime mediated by GABA<sub>B</sub> interneurons may be aborted by using lower stimulation frequen-
cies than a pre-ictal spiking regime mediated by GABA<sub>A</sub> interneurons, since GABA<sub>B</sub> interneu-
rons have slower kinetics than GABA<sub>A</sub> interneurons. Depending on the neuroanatomy and
neurophysiology of a specific brain region (type of subpopulations and connections between
them), activation of specific types of interneurons can be achieved *via* the modulation of differ-
ent neural targets [57, 58, 60, 76–78].

LFPs recorded by SEEG electrodes can provide hints on timescale separations, since there is
a tight link between the form of the observed complex oscillations and timescale separations.
For instance, the pre-ictal bursts examined in this work involve fast spikes followed by a slower
large amplitude discharge. While fast spikes indicate an interaction of fast and slow compo-
nents, the slower large amplitude discharge spikes indicate an interaction of slow and super-
slow components. The natural presence of three distinct timescales in the model, which are
introduced by the kinetics of glutamatergic receptors (NMPA and AMPA) and GABAergic receptors of two distinct types, enables modeling these complex oscillations. However, computation of the synaptic time scales from the LFP is an ill posed problem since the LFP is a complex mixture of excitatory and inhibitory post-synaptic potentials at the level of pyramidal neurons.

Slow-fast analysis of the mathematical models of neural systems with complex oscillatory patterns has contributed to discover the roles of biological ingredients [79–88], unveil the fine structures (e.g. excitability thresholds, spike adding mechanisms and subthreshold oscillations ...etc.) [30,35,89–96], and design controllers [97,98]. Slow-fast thinking has also been insightful for constructing phenomenological models of epileptic dynamics and analyzing their dynamics [99–101]. Recently, response types of brief electrical pulses in coupled NMMs have been investigated using some elements of slow-fast analysis [102]. In [103] a regime of canard solutions has been reported in sleep/wake transitions in a NMM, also in an extended NMM formulation in [84]. As opposed to [82–84], here we reformulated a widely studied NMM in an explicit slow-fast form and unveiled its three-time-scale structure. Thus, the system is an example of three-time-scale systems that are beginning to be explored. Besides, the structure of the model is widely used in engineering studies, in particular in systems with feedback [33]. Hence, the methodology used is this paper could be beneficial in many other research areas.

During our investigations we also observed canard solutions organizing the transition from slow-wave (≈6 Hz) to bursting oscillations through a spike-adding mechanism in between. We did not further explore this interesting mechanism since the main purpose of this paper was to understand the perturbation effect on pre-ictal bursting solutions, which are away from the canard regime in the parameter space. Further analysis concerning the classification of slow dynamics near the fold points, canard solutions and spike-adding mechanisms in the line of [36,37,104] are among the possible extensions of this work.

Another possible avenue to extend this work would be to consider the possibility to perform patient-specific bifurcation analyses of epileptiform patterns to propose patient-specific stimulation parameters (most critically, stimulation frequency) that would result in the abortion of the said epileptiform patterns. Current direct brain stimulation protocols in epilepsy use indeed relatively generic parameters, without consideration for the type, localization or extent of the epileptogenic network; a possible factor to explain the lack of consistency for this therapy so far for drug-refractory epilepsy. For our methods to be applicable an adaptive closed-loop detection system, such as a brain responsive neurostimulation system, can be taken into account, which could detect pre-ictal discharges, and then intervene with an appropriate stimulation.

Finally, we should emphasize that the NMM considered here was initially proposed a model for hippocampal activity. As shown in this study, this NMM can reproduce complex oscillatory patterns at the macroscopic level resulting from interaction of district subpopulations with different kinetics. More recently, the model was shown to have a more general scope as the embedded circuitry is valid that of most of the regions at macroscopic level (see [105] and references there in). Thus, appropriately formulated NMMs and the tools presented here can be used to study the complex dynamics observed in other cortical areas and to investigate effects of external perturbations.

**Supporting information**

S1 Fig. Block diagram of the neural mass model (1). The model features three types of neuronal subpopulations, namely pyramidal neurons (PYR and PYR'), GABAergic SOM+ interneurons (SOM+) and GABAergic PV+ interneurons (PV+). Average PSP at the level of each
subpopulation (denoted by $y$ variables) is determined by a pulse-to-wave function $S(v)$ and a linear dynamic transfer function $h(t)$. Properties of $h(t)$ are determined by synaptic gains ($A$, $B$, $G$) and synaptic time constants ($1/a$, $1/b$, $1/g$). Parameters $C_i$ denote average synaptic contacts. Cortical input is denoted by $p(t)$.

(TIF)

**S2 Fig. Periodic orbits of (4) and (8) for $\varepsilon = 0.01$.** The system is put in the bursting regime by taking $B = 18$ and the other parameters are as given in Table 1. (a) Solution of (3) for $\varepsilon = 0.01$ projected on the bifurcation diagram (black curve) of (4) for $\varepsilon = 0$ where $v_2$ is treated as a parameter. Stable and unstable solutions are indicated with bold and dashed lines, respectively. The equilibrium points along the black Z-shaped curve are unstable on the middle branch of the curve, between the limit points (LP) $LP_1$ and $LP_2$ (black dots), and on the upper branch between the supercritical Hopf (H) bifurcation points $H_1$ and $H_2$ (green dots). The amplitude of the stable limit cycles is bounded by the green continuous curves connecting the $H_1$ and $H_2$ points in the $\varepsilon = 0$ limit. Arrows show the direction of the flow. (b) Solution of (7) for $\varepsilon = 0.01$ projected on the bifurcation diagram of (8) (red curve) where $v_2$ is treated as a parameter. Stable and unstable solutions are indicated with bold and dashed lines, respectively. The equilibrium points along the black Z-shaped curve are unstable on the middle branch of the curve, between the $LP_1$ and $LP_2$ limit points (red dots). Arrows show the direction of the flow.

(TIF)

**S3 Fig. Response of a neural mass block to biphasic balanced pulses.** A neural mass (NM) block, that is $\dot{y} = M/\tau S(I_{ext}(t)) - 2/\tau \dot{y} - 1/\tau^2 y$ with $\tau = \{0.05, 0.01, 0.003\}$ and $M\tau = 1$, receives biphasic pulses $I_{ext}(t)$ at different pulse frequency, amplitude and width. Although $M\tau$ is constant across the trials, the amplitude of the response varies due to the difference between the synaptic kinetics. (a) Amplitude of the steady state oscillations of NM evoked by $I_{ext}(t)$ of pulse width 0.05 ms, amplitude 1 (arb. unit) and frequency in $f = [1, 100] Hz$. Red dashed line is the base level (denoted by $S(0)$) of the NM in the absence of any inputs. The NM with slow kinetics ($\tau = 0.05$) detaches from the base level for a smaller stimulation frequency than the NM with fast kinetics ($\tau = 0.003$). Amplitude of the steady state oscillations decreases with frequency. Panels (a1) and (a2) show the responses at $f = 10$ Hz and $f = 100$ Hz, respectively.

Same color codes are used in panels (a), (a1) and (a2). (b) Amplitude of the steady state oscillations of NM evoked by $I_{ext}(t)$ of pulse width 0.05 ms, frequency $f = 10$ Hz and amplitude $amp = [1,100]$ (arb.unit). Red dashed line is the base level (denoted by $S(0)$) of the NM in the absence of any inputs. The difference between the base line and min ($y(t)$) is larger for the NM with slow kinetics ($\tau = 0.05$) than the NM with fast kinetics ($\tau = 0.003$). Amplitude of the steady state oscillations increases with amplitude then does not change further for $amp > 20$ (arb.unit). Panels (b1) and (b2) show the time course of the responses at $amp = 1$ (arb.unit) and $amp = 20$ (arb.unit), respectively. Same color codes are used in panels (b), (b1) and (b2).

(c) Amplitude of the steady state oscillations of NM evoked by $I_{ext}(t)$ of frequency $f = 10$ Hz, pulse width $width = [0.05, 50]$ ms with amplitude $amp = 0.05 width^{-1}$ (arb.unit). Red dashed line is the base level (denoted by $S(0)$) of the NM in the absence of any inputs. The NM responds with a lower shoot to increasing the pulse width which appears for narrower pulses for the NM with fast kinetics ($\tau = 0.003$) than the NM with slow kinetics ($\tau = 0.03$). Amplitude of the steady state oscillations increases with amplitude then decreases with pulse width. Panels (c1) and (c2) show the time course of the responses at $width = 10$ ms and $width = 50$ ms, respectively. Same color codes are used in panels (c), (c1) and (c2).

(TIF)
## Author Contributions

**Conceptualization:** Elif Köksal Ersöz, Fabrice Wendling.

**Data curation:** Elif Köksal Ersöz, Fabrice Bartolomei, Fabrice Wendling.

**Formal analysis:** Elif Köksal Ersöz.

**Funding acquisition:** Fabrice Wendling.

**Investigation:** Elif Köksal Ersöz, Julien Modolo, Fabrice Bartolomei, Fabrice Wendling.

**Methodology:** Elif Köksal Ersöz.

**Project administration:** Elif Köksal Ersöz, Fabrice Wendling.

**Resources:** Fabrice Bartolomei.

**Software:** Elif Köksal Ersöz.

**Supervision:** Elif Köksal Ersöz, Julien Modolo, Fabrice Wendling.

**Validation:** Elif Köksal Ersöz, Fabrice Bartolomei.

**Visualization:** Elif Köksal Ersöz.

**Writing – original draft:** Elif Köksal Ersöz, Fabrice Wendling.

**Writing – review & editing:** Elif Köksal Ersöz, Julien Modolo, Fabrice Bartolomei, Fabrice Wendling.

## References

1. Billakota S, Devinsky O, Kim K-W. Why we urgently need improved epilepsy therapies for adult patients. Neuropharmacology. 2019;107855. https://doi.org/10.1016/j.neuropharm.2019.107855 PMID: 31751547

2. Kwan P, Brodie MJ. Early Identification of Refractory Epilepsy. N Engl J Med. 2000; 342 (5):314–319 https://doi.org/10.1056/NEJM200002033420503 PMID: 10660394

3. Avoli M, de Curtis M, Gnatkovsky V, Gotman J, Köhling R, Lévesque M, et al. Specific imbalance of excitatory/inhibitory signaling establishes seizure onset pattern in temporal lobe epilepsy. J Neurophysiol. 2016; 115 (6):3229–3237 https://doi.org/10.1152/jn.01128.2015 PMID: 27075542

4. Kahle KT, Khanna AR, Duan J, Staley KJ, Delpire E, Poduri A. The KCC2 Cotransporter and Human Epilepsy: Getting Excited About Inhibition. Neuroscientist. 2016; 22 (6):555–562 https://doi.org/10.1177/1073858416645087 PMID: 27130838

5. Dudek FE, Staley KJ. The Time Course and Circuit Mechanisms of Acquired Epileptogenesis. In: Noe-bels JL, Avoli M, Rogawski MA, Olsen RW, Delgado-Escueta AV, editors. Jasper's Basic Mechanisms of the Epilepsies. 4th ed. Bethesda (MD): National Center for Biotechnology Information (US); 2012.

6. Morris H, Najm I, Kahane P. Epilepsy surgery: patient selection. In: Lüders HO, editor. Textbook of Epilepsy Surgery. London: Informa Healthcare; 2008. p. 230–7.

7. Fisher RS, Velasco AL. Electrical brain stimulation for epilepsy. Nat Rev Neurol. 2014; 10 (5):261–270. https://doi.org/10.1038/nrneurol.2014.59 PMID: 24709892

8. Sprengers M, Vonck K, Carrette E, Marson AG, Boon P. Deep brain and cortical stimulation for epilepsy. Cochrane Database Syst Rev. 2017; 7:CD008497. https://doi.org/10.1002/14651858. CD008497.pub3 PMID: 28718878

9. Lesser RP, Lee HW, Webber WRS, Prince B, Crone NE, Miglioretti DL. Short-term variations in response distribution to cortical stimulation. Brain. 2008; 131 (6):1528–1539. https://doi.org/10.1093/brain/awn044 PMID: 18337272

10. Laxpati NG, Kasoff WS, Gross RE. Deep Brain Stimulation for the Treatment of Epilepsy: Circuits, Targets, and Trials. Neurotherapeutics. 2014; 11 (3):508–526. https://doi.org/10.1007/s13311-014-0279-9 PMID: 24957200

11. Heck CN, King-Stephens D, Massey AD, Nair DR, Jobst BC, Barkley GL, et al. Two-year seizure reduction in adults with medically intractable partial onset epilepsy treated with responsive
neurostimulation: Final results of the RNS System Pivotal trial. Epilepsia. 2014; 55(3):432-441.
https://doi.org/10.1111/epi.12534 PMID: 24621228

12. Lee S-A, Spencer DD, Spencer SS. Intracranial EEG Seizure-Onset Patterns in Neocortical Epilepsy. Epilepsia. 2000; 41 (3):297-307. https://doi.org/10.1111/j.1528-1157.2000.tb00159.x PMID: 10714401

13. Bartolomei F, Wendling F, Régis J, Gavaret M, Guye M, Chauvel P. Pre-ictal synchronicity in limbic networks of mesial temporal lobe epilepsy. Epilepsy Res. 2004; 61(1):89-104. https://doi.org/10.1016/j.elepsyres.2004.06.006 PMID: 15451011

14. Velascol AL, Wilson CL, Babb TL, Engel J Jr. Functional and Anatomic Correlates of Two Frequently Observed Temporal Lobe Seizure-Onset Patterns. Neural Plast. 2000; 7:201878.

15. Fenichel N. Geometric singular perturbation theory for ordinary differential equations. J Differ Equ. 1979; 31(1):53-98.

16. Jones CKRT. Geometric singular perturbation theory. In: Arnold L, Jones CKRT, Mischaikow K, Rau-gel G, Johnson R, editors. Dynamical Systems: Lectures Given at the 2nd Session of the Centro Inter-nazionale Matematico Estivo (CIME) held in Montecatini Terme, Italy, June 13-22, 1994. Berlin, Heidelberg: Springer; 1995. p. 44-118. (Lecture Notes in Mathematics).

17. Freeman WJ. Simulation of chaotic EEG patterns with a dynamic model of the olfactory system. Biol Cybern. 1987; 56(2-3):139-150.

18. Lopes da Silva FH, Hoeks A, Smits H, Zetterberg LH. Model of brain rhythmic activity: The alpha-rhythm of the thalamus. Kybernetik. 1974; 15(1):27-37. https://doi.org/10.1007/BF00270757 PMID: 4853232

19. Ursino M, Cona F, Zavaglia M. The generation of rhythms within a cortical region: Analysis of a neural mass model. NeurolImage. 2010; 52(3):1080-1094. https://doi.org/10.1016/j.neuroimage.2009.12.084 PMID: 20045071

20. Jansen BH, Rit VG. Electroencephalogram and visual evoked potential generation in a mathematical model of coupled cortical columns. Biol Cybern. 1995; 73: 357-366.

21. Suffczynski P, Kalitzin S. Lopes Da Silva FH. Dynamics of non-convulsive epileptic phenomena mod- eled by a bistable neuronal network. Neuroscience. 2004; 126(2):467-484. https://doi.org/10.1016/j.neuroscience.2004.03.014 PMID: 15207365

22. Wendling F, Bartolomei F, Bellanger JJ, Chauvel P. Epileptic fast activity can be explained by a model of impaired GABAergic dendritic inhibition. Eur J Neurosci. 2002; 15(9):1499-1508. https://doi.org/10.1046/j.1460-9568.2002.01985.x PMID: 12028360

23. Rodrigues S, Barton D, Szalai R, Benjamin O, Richardson MP, Terry JR. Transitions to spike-wave oscillations and epileptic dynamics in a human cortico-thalamic mean-field model. J Comput Neurosci. 2009; 27(3):507-526. https://doi.org/10.1007/s10827-009-0166-2 PMID: 19499316

24. Goodfellow M, Schindler K, Baier G. Self-organised transients in a neural mass model of epileptogenic tissue dynamics. NeurolImage. 2012; 59(3):2644-2660. https://doi.org/10.1016/j.neuroimage.2011.08.060 PMID: 21945465

25. Bhattacharya BS, Coyle D, Maguire LP. A thalamo-cortico-thalamic neural mass model to study alpha rhythms in Alzheimer's disease. Neural Netw. 2011; 24(6):631-645. https://doi.org/10.1016/j.neunet.2011.02.009 PMID: 21435838

26. Liu F, Wang J, Liu C, Li H, Deng B, Fietkiewicz C, et al. A neural mass model of basal ganglia nuclei simulates pathological beta rhythm in Parkinson's disease. Chaos Interdiscip J Nonlinear Sci. 2016; 26(12):123113. https://doi.org/10.1063/1.4972200 PMID: 28039987

27. Molaee-Ardekani B, Benquet P, Bartolomei F, Wendling F. Computational modeling of high-frequency oscillations at the onset of neocortical partial seizures: From 'altered structure' to 'dysfunction'. Neuro-Image. 2010; 52(3):1109-1122. https://doi.org/10.1016/j.neuroimage.2009.12.049 PMID: 20034581

28. Lopes da Silva F. Neural mechanisms underlying brain waves: from neural membranes to networks. Electroencephalogr Clin Neurophysiol. 1991; 79(2):81-93. https://doi.org/10.1016/0013-4694(91)90044-5 PMID: 1713832

29. Wendling F, Gerber U, Cosandier-Rimele D, Nica A, De Montigny J, Raineteau O, et al. Brain (Hyper) Excitability Revealed by Optimal Electrical Stimulation of GABAergic Interneurons. Brain Stimulat. 2016; 9(6):919-932. https://doi.org/10.1016/j.brs.2016.07.001 PMID: 27576186

30. Bertram R, Rubin JE. Multi-timescale systems and fast-slow analysis. Math Biosci. 2017; 287:105-121.

31. Rinzel J. A Formal Classification of Bursting Mechanisms in Excitable Systems. In: Teramoto E, Yumaguti M, editors. Mathematical Topics in Population Biology, Morphogenesis and Neurosciences. Lecture Notes in Biomathematics, vol 71. Berlin, Heidelberg: Springer Berlin Heidelberg; 1987. p. 267-281. https://doi.org/10.1007/978-3-642-93360-8_26
32. Izhikevich EM. Neural excitability, spiking and bursting. Int J Bifurc Chaos. 2000; 10(06):1171-1266.

33. Kuehn C. Multiple Time Scale Dynamics. Applied Mathematical Sciences. Springer International Publishing; 2015.

34. Krupa M, Szmolyan P. Extending Geometric Singular Perturbation Theory to Nonhyperbolic Points—Fold and Canard Points in Two Dimensions. SIAM J Math Anal. 2001; 33(2):286-314.

35. Rubin J, Wechselberger M. Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model. Biol Cybern. 2007; 97(1):5-32. https://doi.org/10.1007/s00422-007-0153-5 PMID: 17458557

36. Desroches M, Kirk V. Spike-Adding in a Canonical Three-Time-Scale Model: Superslow Explosion and Folded-Saddle Canards. SIAM J Appl Dyn Syst. 2018; 17(3):1989-2017.

37. Nan P, Wang Y, Kirk V, Rubin JE. Understanding and Distinguishing Three-Time-Scale Oscillations: Case Study in a Coupled Morris-Lecar System. SIAM J Appl Dyn Syst. 2015; 14(3):1518-1557.

38. Suffczynski P, Kalitzin S, da Silva FL, Parra J, Velis D, Wendling F. Active paradigms of seizure anticipation: Computer model evidence for necessity of stimulation. Phys Rev E. 2008; 78(5).

39. Talairach J, Bancaud J, Szikla G, Bonis A, Geier S, Vedrenne C. New approach to the neurosurgery of epilepsy. Stereotaxic methodology and therapeutic results. 1. Introduction and history. Neurochirurgie. 1974; 20(Suppl 1):1-240. PMID: 4603712

40. Doedel EJ, Champneys A, Fairgrieve TF, Yu AB, Kuznetsov AP, Oldeman BE, et al. Auto-07p: Continuation and bifurcation software for ordinary differential equations. 2007. https://doi.org/10.1103/PhysRevE.75.016213 PMID: 17358241

41. Ermentrout B. Simulating, Analyzing, and Animating Dynamical Systems [Internet]. Software, Environments and Tools. Society for Industrial and Applied Mathematics; 2002.

42. Bartolomei F, Chauvel P, Wendling F. Epileptogenicity of brain structures in human temporal lobe epilepsy: a quantified study from intracerebral EEG. Brain. 2008; 131 (7):1818-1830. https://doi.org/10.1093/brain/awn111 PMID: 18556663

43. Andrzejak RG, David O, Gnatkovsky V, Wendling F, Bartolomei F, Francione S, et al. Localization of Epileptogenic Zone on Pre-surgical Intracranial EEG Recordings: Toward a Validation of Quantitative Signal Analysis Approaches. Brain Topogr. 2015; 28(6):832-837. https://doi.org/10.1007/s10548-014-0380-8 PMID: 24929558

44. Lagarde S, Buzori S, Trebuchon A, Carron R, Scavarda D, Milh M, et al. The repertoire of seizure onset patterns in human focal epilepsies: Determinants and prognostic values. Epilepsia. 2019; 60 (1):85-95. https://doi.org/10.1111/epi.14604 PMID: 30426477

45. Wendling F, Bellanger J-J, Badier J-M, Coatrieux J-L. Extraction of spatio-temporal signatures from depth EEG seizure signals based on objective matching in warped vectorial observations. IEEE Trans Biomed Eng. 1996; 43(10):990-1000. https://doi.org/10.1109/10.536900 PMID: 9214816

46. Wendling F, Badier JM, Chauvel P, Coatrieux JL. A method to quantify invariant information in depth-recorded epileptic seizures. Electroencephalogr Clin Neurophysiol. 1997; 102(6):472-485. https://doi.org/10.1016/s0013-4694(96)96633-3 PMID: 9216480

47. Wendling F, Benquet P, Bartolomei F, Jirsa V. Computational models of epileptiform activity. J Neurosci Methods. 2016; 260:233-251. https://doi.org/10.1016/j.jneumeth.2015.03.027 PMID: 25843066

48. da Silva FHL, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Dynamical diseases of brain systems: different routes to epileptic seizures. IEEE Trans Biomed Eng. 2003; 50 (5):540-548. https://doi.org/10.1109/TBME.2003.810703 PMID: 12769430

49. Silva FLD, Blanes W, Kalitzin SN, Parra J, Suffczynski P, Velis DN. Epilepsies as Dynamical Diseases of Brain Systems: Basic Models of the Transition Between Normal and Epileptic Activity. Epilepsia. 2003; 44(s12):72-83.

50. Meisel C, Kuehn C. Scaling Effects and Spatio-Temporal Multilevel Dynamics in Epileptic Seizures. PLoS One. 2012; 7(2):e30371. https://doi.org/10.1371/journal.pone.0030371 PMID: 22363431

51. Perucca P, Dubeau F, Gotman J. Intracranial electroencephalographic seizure-onset patterns: effect of underlying pathology. Brain J Neurol. 2014; 137(Pt 1):183-196. https://doi.org/10.1093/brain/awt299 PMID: 24176980

52. Lagarde S, Bonini F, McGonigal A, Chauvel P, Gavaret M, Scavarda D, et al. Seizure-onset patterns in focal cortical dysplasia and neurodevelopmental tumors: Relationship with surgical prognosis and neuropathologic subtypes. Epilepsia. 2016; 57(9):1426-1435. https://doi.org/10.1111/epi.13464 PMID: 27406939

53. Huberfeld G, Menendez de la Prida L, Pallud J, Cohen I, Le Van Quyen M, Adam C, et al. Glutamatergic pre-ictal discharges emerge at the transition to seizure in human epilepsy. Nat Neurosci. 2011; 14 (5):627-634. https://doi.org/10.1038/nn.2790 PMID: 21460834
54. Avoli M, de Curtis M. GABAergic synchronization in the limbic system and its role in the generation of epileptiform activity. Prog Neurobiol. 2011; 95 (2):104–132. https://doi.org/10.1016/j.pneurobio.2011. 07.003 PMID: 21802488

55. Lasztóczi B, Nyitrai G, Héja L, Kardos J. Synchronization of GABAergic Inputs to CA3 Pyramidal Cells Precedes Seizure-Like Event Onset in Juvenile Rat Hippocampal Slices. J Neurophysiol. 2009; 102 (4):2538–2553. https://doi.org/10.1152/jn.91318.2008 PMID: 19675286

56. Muldoon SF, Pasqualetti F, Gu S, Cieslak M, Grafton ST, Vettel JM, et al. Stimulation-Based Control of Dynamic Brain Networks. PLoS Comput Biol. 2016 Sep9; 12(9). https://doi.org/10.1371/journal. pcbi.1005076 PMID: 27611328

57. Chiang C-C, Ladas TP, Gonzalez-Reyes LE, Durand DM. Seizure Suppression by High Frequency Optogenetic Stimulation Using In Vitro and In Vivo Animal Models of Epilepsy. Brain Stimulat. 2014; 7 (6):890–899. https://doi.org/10.1016/j.brs.2014.07.034 PMID: 25108607

58. Ladas TP, Chiang C-C, Gonzalez-Reyes LE, Nowak T, Durand DM. Seizure reduction through inter-
neuron-mediated entrainment using low frequency optical stimulation. Exp Neurol. 2015; 269:120-
132. https://doi.org/10.1016/j.expneurol.2015.04.001 PMID: 25863022

59. Rashid S, Pho G, Czigler M, Werz MA, Durand DM. Low frequency stimulation of ventral hippocampal commissures reduces seizures in a rat model of chronic temporal lobe epilepsy. Epilepsia. 2012; 53 (1):147–156. https://doi.org/10.1111/j.1528-1167.2011.03348.x PMID: 22150779

60. Toprani S, Durand DM. Long-lasting hyperpolarization underlies seizure reduction by low frequency deep brain electrical stimulation. J Physiol. 2013; 591(22):5765–5790. https://doi.org/10.1113/ jphysiol.2013.253757 PMID: 23981713

61. Kano T, Inaba Y, D'Antuono M, Biagini G, Levésque M, Avoli M. Blockade of in vitro ctogenesis by low-frequency stimulation coincides with increased epileptiform response latency. J Neurophysiol. 2015; 114(1):21–28. https://doi.org/10.1152/jn.00248.2015 PMID: 25925325

62. Svejgaard B, Andreasen M, Nedergaard S. Role of GABAB receptors in proepileptic and antiepileptic effects of an applied electric field in rat hippocampus in vitro. Brain Res. 1710; 2019:157–162.

63. Couturier NH, Durand DM. Corpus callosum low-frequency stimulation suppresses seizures in an acute rat model of focal cortical seizures. Epilepsia. 2018; 59(12):2219–2230. https://doi.org/10.1111/ epi.14595 PMID: 30426470

64. Couturier NH, Durand DM. Comparison of fiber tract low frequency stimulation to focal and ANT stimu-
lation in an acute rat model of focal cortical seizures. Brain Stimulat. 2019; 13(2):499–506.

65. Li X-G, Somogyi P, Ylinen A, Buzsáki G. The hippocampal CA3 network: An in vivo intracellular label-
ing study. J Comp Neurol. 1994; 339(2):181–208. https://doi.org/10.1002/cne.903390204 PMID:
8300905

66. Sik A, Penttonen M, Ylinen A, Buzsaki G. Hippocampal CA1 interneurons: an in vivo intracellular label-
ing study. J Neurosci. 1995; 15(10):6651–6665. https://doi.org/10.1523/JNEUROSCI.15-10-06651.
1995 PMID: 7472426

67. Assaf F, Schiller Y. The antiepileptic and ctogenic effects of optogenetic neurostimulation of PV-
expressing interneurons. J Neurophysiol. 2016; 116(4):1694–1704. https://doi.org/10.1152/jn.00744.
2015 PMID: 27486107

68. Lévesque M, Chen L-Y, Etter G, Shiri Z, Wang S, Williams S, et al. Paradoxical effects of optogenetic stimulation in mesial temporal lobe epilepsy. Ann Neurol. 2019; 86(5):714–728. https://doi.org/10. 1002/ana.25572 PMID: 31393618

69. Krook-Magnuson E. The Devil's in the Details: How to Harness Inhibition for Seizure Control. Epilepsy Curr. 2020;1535759720902078. https://doi.org/10.1177/1535759720902078 PMID: 32313505

70. Zangiabadi N, Ladino LD, Sina F, Orozco-Hernández JP, Carter A, Téllez-Zenteno JF. Deep Brain Stimulation and Drug-Resistant Epilepsy: A Review of the Literature. Front Neurol. 2019; 10. https://doi.org/10.3389/fneur.2019.00601 PMID: 31244761

71. Yamamoto J, Ikeda A, Kinoshita M, Matsumoto R, Satow T, Takeshita K, et al. Low-frequency electric cortical stimulation decreases interictal and ictal activity in human epilepsy. Seizure. 2006; 15(7):520–527. https://doi.org/10.1016/j.seizure.2006.06.004 PMID: 16908203

72. Schiller Y, Bankirer Y. Cellular Mechanisms Underlying Antiepileptic Effects of Low- and High-Fre-
quency Electrical Stimulation in Acute Epilepsy in Neocortical Brain Slices In Vitro. J Neurophysiol.
2007; 97(3):1887–1902. https://doi.org/10.1152/jn.00514.2006 PMID: 17151229

73. Chiang C-C, Lin C-CK JM-S, Durand DM. High frequency stimulation can suppress globally seizures induced by 4-AP in the rat hippocampus: An acute in vivo study. Brain Stimulat. 2013; 6(2):180–189.

74. Velasco AL, Velasco M, Velasco F, Menes D, Gordon F, Rocha L, et al. Subacute and Chronic Electri-
cal Stimulation of the Hippocampus on Intractable Temporal Lobe Seizures: Preliminary Report. Arch
Med Res. 2000;13. https://doi.org/10.1016/s0188-4409(00)00064-3 PMID: 11036183
75. Taylor PN, Wang Y, Goodfellow M, Dauwels J, Moeller F, Stephani U, et al. A Computational Study of Stimulus Driven Epileptic Seizure Abatement. PLoS One. 2014; 9(12):e114316. https://doi.org/10.1371/journal.pone.0114316 PMID: 25531883

76. Sah N, Sikdar SK. Transition in subicular burst firing neurons from epileptiform activity to suppressed state by feedforward inhibition. Eur J Neurosci. 2013; 38(4):2542-2556. https://doi.org/10.1111/ejn. 12262 PMID: 23725217

77. Xu Z, Wang Y, Chen B, Xu C, Wu X, Wang Y, et al. Entorhinal Principal Neurons Mediate Brain-stimulation Treatments for Epilepsy. EBioMedicine. 2016; 14:148-160. https://doi.org/10.1016/j.ebiom. 2016.11.027 PMID: 27908611

78. Krook-Magnuson E, Armstrong C, Oijala M, Soltesz I. On-demand optogenetic control of spontaneous seizures in temporal lobe epilepsy. Nat Commun. 2013; 4:1376. https://doi.org/10.1038/ncomms2376 PMID: 23340416

79. Bacak BJ, Kim T, Smith JC, Rubin JE, Rybak IA. Mixed-mode oscillations and population bursting in the pre-Bötzinger complex. elife. 2016; 5:e13403. https://doi.org/10.7554/eLife.13403 PMID: 26974345

80. Tabak J, Rinzel J, O'Donovan MJ. The Role of Activity-Dependent Network Depression in the Expression and Self-Regulation of Spontaneous Activity in the Developing Spinal Cord. J Neurosci. 2001; 21 (22):8966-8978. https://doi.org/10.1523/JNEUROSCI.21-22-08966.2001 PMID: 11698607

81. Tabak J, Senn W, O'Donovan MJ, Rinzel J. Modeling of Spontaneous Activity in Developing Spinal Cord Using Activity-Dependent Depression in an Excitatory Network. J Neurosci. 2000; 20(8):3041-3056. https://doi.org/10.1523/JNEUROSCI.20-08-03041.2000 PMID: 10751456

82. Nunemaker CS, Bertram R, Sherman A, Tsaneva-Atanasova K, Daniel CR, Satin LS. Glucose Modulates [Ca2+]i Oscillations in Pancreatic Islets via Ionic and Glycolytic Mechanisms. Biophys J. 2006; 91(6):2082-2096. https://doi.org/10.1529/biophysj.106.087296 PMID: 16815907

83. Lin CC, Segel LA. Mathematics Applied to Deterministic Problems in the Natural Sciences. Society for Industrial and Applied Mathematics; 1988. p 621. (Classics in Applied Mathematics).

84. Duan W, Lee K, Herbison AE, Sneyd J. A mathematical model of adult GnRH neurons in mouse brain and its bifurcation analysis. J Theor Biol. 2011; 276(1):22-34. https://doi.org/10.1016/j.jtbi.2011.01. 035 PMID: 21300070

85. Rubin JE, Hayes JA, Mendenhall JL, Negro CAD. Calcium-activated nonspecific cation current and synaptic depression promote network-dependent burst oscillations. Proc Natl Acad Sci. 2009; 106 (8):2939-2944. https://doi.org/10.1073/pnas.0808776106 PMID: 19196976

86. FitzHugh R. Impulses and Physiological States in Theoretical Models of Nerve Membrane. Biophys J. 1961; 1(6):445-466. https://doi.org/10.1016/s0006-3495(61)86902-6 PMID: 19431309

87. Bertram R, Previte J, Sherman A, Kinard TA, Satin LS. The Phantom Burster Model for Pancreatic β-Cells. Biophys J. 2000; 79(6):2880-2892. https://doi.org/10.1016/S0006-3495(00)76525-8 PMID: 11106596

88. Rotstein HG, Pervouchine DD, Acker CD, Gillies MJ, White JA, Buhl EH, et al. Slow and Fast Inhibition and an H-Current Interact to Create a Theta Rhythm in a Model of CA1 Interneuron Network. J Neurophysiol. 2005; 94(2):1509-1518. https://doi.org/10.1152/jn.00957.2004 PMID: 15857967

89. Harvey E, Kirk V, Wechselberger M, Sneyd J. Multiple Timescales, Mixed Mode Oscillations and Canards in Models of Intracellular Calcium Dynamics. J Nonlinear Sci. 2011; 21(5):639.

90. Krupa M, Popović N, Kopell N, Rotstein HG. Mixed-mode oscillations in a three time-scale model for the dopaminergic neuron. Chaos Interdiscip J Nonlinear Sci. 2008; 18(1):015106.

91. Rotstein HG, Oppermann T, White JA, Kopell N. The dynamic structure underlying subthreshold oscillatory activity and the onset of spikes in a model of medial entorhinal cortex stellate cells. J Comput Neurosci. 2006 Dec; 21(3):271-292. https://doi.org/10.1007/s10827-006-8096-8 PMID: 16927211

92. Rotstein HG, Wechselberger M, Kopell N. Canard Induced Mixed-Mode Oscillations in a Medial Entorhinal Cortex Layer II Stellate Cell Model. SIAM J Appl Dyn Syst. 2008; 7(4):1582-1611.

93. De Maesschalck P, Desroches M. Numerical Continuation Techniques for Planar Slow-Fast Systems. SIAM J Appl Dyn Syst. 2013; 12(3):1159-1180.

94. De Maesschalck P, Wechselberger M. Neural Excitability and Singular Bifurcations. J Math Neurosci JMN. 2015; 5(1):16. https://doi.org/10.1186/s13408-015-0029-2 PMID: 26246435

95. Mitry J, McCarthy M, Kopell N, Wechselberger M. Excitable Neurons, Firing Threshold Manifolds and Canards. J Math Neurosci. 2013; 3(1):12.

96. Köksal Ersöz E, Desroches M, Guillamon A, Rinzel J, Tabak J. Canard-induced complex oscillations in an excitatory network. J Math Biol. 2020; 80:2075-2107.
97. Rubin JE, Terman D. High Frequency Stimulation of the Subthalamic Nucleus Eliminates Pathological Thalamic Rhythmicity in a Computational Model. J Comput Neurosci. 2004; 16(3):211-235.

98. Wilson D, Moehlis J. A Hamilton-Jacobi-Bellman approach for termination of seizure-like bursting. J Comput Neurosci. 2014; 37(2):345-355.

99. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014; 137(8):2210-2230.

100. Proix T, Bartolomei F, Chauvel P, Bernard C, Jirsa VK. Permittivity Coupling across Brain Regions Determines Seizure Recruitment in Partial Epilepsy. J Neurosci. 2014; 34(45):15009-15021. https://doi.org/10.1523/JNEUROSCI.1570-14.2014 PMID: 25378166

101. Saggio ML, Spiegler A, Bernard C, Jirsa VK. Fast-Slow Bursters in the Unfolding of a High Codimension Singularity and the Ultra-slow Transitions of Classes. J Math Neurosci. 2017; 7(1).

102. Hebbink J, van Gils SA, Meijer HGE. On analysis of inputs triggering large nonlinear neural responses Slow-fast dynamics in the Wendling neural mass model. Commun Nonlinear Sci Numer Simul. 2020; 83:105103.

103. Weigenand A, Costa MS, Ngo H-VV, Claussen JC, Martinetz T. Characterization of K-Complexes and Slow Wave Activity in a Neural Mass Model. PLoS Comput Biol. 2014; 10(11):e1003923. https://doi.org/10.1371/journal.pcbi.1003923 PMID: 25392991

104. Vo T, Bertram R, Wechselberger M. Multiple Geometric Viewpoints of Mixed Mode Dynamics Associated with Pseudo-plateau Bursting. SIAM J Appl Dyn Syst. 2013; 12(2):789-830.

105. Bensaid S, Modolo J, Merlet I, Wendling F, Benquet P. COALIA: A Computational Model of Human EEG for Consciousness Research. Front Syst Neurosci. 2019; 13. https://doi.org/10.3389/fnsys.2019.00059 PMID: 31798421