# A unified physiological framework of transitions between seizures, sustained ictal activity and depolarization block at the single neuron level

Damien Depannemaecker, Anton Ivanov, Davide Lillo, Len Spek, Christophe

Bernard, Viktor Jirsa

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

## To cite this version:

Damien Depannemaecker, Anton Ivanov, Davide Lillo, Len Spek, Christophe Bernard, et al.. A unified physiological framework of transitions between seizures, sustained ictal activity and depolarization block at the single neuron level. Journal of Computational Neuroscience, 2022, 10.1007/s10827-022-00811-1. hal-02980442v2

HAL Id: hal-02980442

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

Submitted on 20 Jan 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A unified physiological framework of transitions between seizures, sustained ictal activity and depolarization block at the single neuron level

Damien Depannemaecker¹,²,³ · Anton Ivanov¹ · Davide Lillo¹ · Len Spek¹,⁴ · Christophe Bernard¹ · Viktor Jirsa¹

Received: 21 February 2021 / Revised: 10 November 2021 / Accepted: 3 January 2022

© The Author(s) 2022

## Abstract

The majority of seizures recorded in humans and experimental animal models can be described by a generic phenomenological mathematical model, the Epileptor. In this model, seizure-like events (SLEs) are driven by a slow variable and occur via saddle node (SN) and homoclinic bifurcations at seizure onset and offset, respectively. Here we investigated SLEs at the single cell level using a biophysically relevant neuron model including a slow/fast system of four equations. The two equations for the slow subsystem describe ion concentration variations and the two equations of the fast subsystem delineate the electrophysiological activities of the neuron. Using extracellular $K^+$ as a slow variable, we report that SLEs with SN/homoclinic bifurcations can readily occur at the single cell level when extracellular $K^+$ reaches a critical value. In patients and experimental models, seizures can also evolve into sustained ictal activity (SIA) and depolarization block (DB), activities which are also parts of the dynamic repertoire of the Epileptor. Increasing extracellular concentration of $K^+$ in the model to values found during experimental status epilepticus and DB, we show that SIA and DB can also occur at the single cell level. Thus, seizures, SIA, and DB, which have been first identified as network events, can exist in a unified framework of a biophysical model at the single neuron level and exhibit similar dynamics as observed in the Epileptor.

**Author Summary:** Epilepsy is a neurological disorder characterized by the occurrence of seizures. Seizures have been characterized in patients in experimental models at both macroscopic and microscopic scales using electrophysiological recordings. Experimental works allowed the establishment of a detailed taxonomy of seizures, which can be described by mathematical models. We can distinguish two main types of models. Phenomenological (generic) models have few parameters and variables and permit detailed dynamical studies often capturing a majority of activities observed in experimental conditions. But they also have abstract parameters, making biological interpretation difficult. Biophysical models, on the other hand, use a large number of variables and parameters due to the complexity of the biological systems they represent. Because of the multiplicity of solutions, it is difficult to extract general dynamical rules. In the present work, we integrate both approaches and reduce a detailed biophysical model to sufficiently low-dimensional equations, and thus maintaining the advantages of a generic model. We propose, at the single cell level, a unified framework of different pathological activities that are seizures, depolarization block, and sustained ictal activity.

**Keywords** Biophysical neuron model · Seizures · Status epilepticus · Depolarization block · Bifurcations · Slow-fast system

## 1 Introduction

Since seizures can be triggered in most brain regions from most species, it has been proposed that simple mathematical rules should be sufficient to describe such a basic form of physiological activity, particularly their dynamics. Several conceptual frameworks have been proposed to explain seizure dynamics (Depannemaecker et al., 2021; Naze, 2015; Naze et al., 2015; Soltesz & Staley, 2008; Staley, 2015; Stefanescu et al., 2012; Y. Wang et al., 2017; Wendling et al.,

Action Editor: Steven J. Schiff.

Christophe Bernard and Viktor Jirsa equally contributing last authors.

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

Damien Depannemaecker
damien.d@cnrs.fr

Extended author information available on the last page of the article
2016). The predominant framework posits that the major-
ity of seizure onsets and offsets correspond to bifurcations
(Jirsa et al., 2014; Saggio et al., 2017), although there exist
other non-bifurcation types (Blenkinsop et al., 2012). This
framework has been generalized by Saggio and colleagues
(Saggio et al., 2017, 2020). A phenomenological mathemati-
cal model, called the Epileptor, describes the dynamics of
a majority of seizures recorded in drug-resistant patients,
and most seizures recorded in experimental models (Jirsa
et al., 2014; Saggio et al., 2020). A qualitative analysis of the
Epileptor reveals that seizures, sustained ictal activity (SIA)
and depolarization block (DB) co-exist, and that multiple
types of transitions from one type of activity to the other
are possible; predictions that were verified experimentally
(El Houssaini et al., 2015; Houssaini et al., 2020; Saggio
et al., 2017). Importantly, the majority of seizures recorded
in patients and experimental models are characterized by a
saddle node (SN) bifurcation at the onset and an homoclinic
bifurcation at the offset (Jirsa et al., 2014). Since it is phe-
nomenological, the Epileptor model does not provide direct
insight regarding underlying biophysical mechanisms.

Detailed biophysical network models have been devel-
oped to investigate seizure mechanisms (Rodrigues et al.,
2015; Santhakumar et al., 2005; Tejada et al., 2014).
These models contain too many parameters and vari-
ables to perform a detailed analysis of their dynamics
repertoire. Although seizures have always been identi-
fied experimentally as network events, their equivalent
in terms of dynamics can be observed at the single cell
level in biophysical models (Bikson et al., 2003; Bragin
et al., 1997; Chizhov et al., 2018; Cressman et al., 2009;
Hübel & Dahlem, 2014; Kager et al., 2000; Lietsche et al.,
2016; McCormick & Contreras, 2001). This suggests that,
in terms of dynamics, core mechanisms already exist at
the single cell level to generate SLE and DB. Obviously,
investigating mechanisms is easier at the single cell level
with a biophysical model than in a network of hundreds of
connected neurons. Building on the proposal that bursting
activity in neurons can be described in terms of bifurca-
tions (E. Izhikevich, 2007; E. M. Izhikevich, 2000), dif-
ferent biophysical single cell models have been proposed
to study SLE and DB, but not SIA (Barreto & Cressman,
2011; Chizhov et al., 2018; Cressman et al., 2009; Hübel
& Dahlem, 2014; Kager et al., 2000; Øyehaug et al.,
2012; Ullah & Schiff, 2010; Wei et al., 2014a, b). They
are slow/fast systems, where a slow subsystem drives the
fast subsystem between different states. In such models,
the studied fast subsystem delineates the neuronal mem-
brane electrophysiological activities. The slow subsystem
can be represented by variations of different slow vari-
ables including ion concentration (Barreto & Cressman,
2011; Chizhov et al., 2018; Cressman et al., 2009; Hübel
& Dahlem, 2014; Kager et al., 2000; Øyehaug et al., 2012;

Wei et al., 2014a, b), oxygen level (Wei et al., 2014a, b), volume (Øyehaug et al., 2012; Wei et al., 2014a, b) and interaction with glial cells (Hübel & Dahlem, 2014; Øyehaug et al., 2012). These models provide mechanistic insights, in particular how the slow variable influences neuronal activity, including the transitions from “healthy” regimes to “pathological” ones like SLEs and DB. However, none of these models show a bursting pattern with a SN bifurcation at the onset and an homoclinic bifurcation at the offset of the event. The goal of the present study was to find the minimal conditions to account for SLE, SIA and DB in a Hodgkin-Huxley-like single cell model, constrained by SN and homoclinic bifurcations at onset and offset, respectively.

A variable acting on a slow time scale is necessary to drive
the system through different activities (e.g. from SLE to DB).
Fluctuations of ion concentrations in the extracellular space
modulate the electrophysiological activity of a single neuron
(Cressman et al., 2009; Wei et al., 2014a, b). The present
work focuses on extracellular potassium concentration ([K]₀)
because it increases during seizures (de Curtis et al., 2018;
Fisher et al., 1976; Fröhlich et al., 2008; Lux et al., 1986; L.
Wang et al., 2016), even in the absence of synaptic activi-
ties (de Almeida et al., 2008; Jefferys & Haas, 1982). In
addition, in experimental models, the transition to DB cor-
relates with a much larger increase of [K]₀ as compared to
SLEs (El Houssaini et al., 2015; Gloveli et al., 1995). Theo-
retical works show that potassium could be responsible for
local synchronization (Durand et al., 2010) and that it is an
important parameter controlling neural dynamics (Barreto &
Cressman, 2011; Cressman et al., 2009; Ullah & Schiff, 2010;
Wei et al., 2014a, b). We here consider the slow modulatory
effects of [K]₀ variations. In our model, the slow sub-system
describes ionic concentration variations. The fast subsystem
characterizes the dynamics of trans-membrane ion flows
through voltage-gated and the sodium-potassium pump, and
so allows tracing the membrane potential. We report that this
single cell model accounts for the SN/homoclinic bifurcation
pair and that it reproduces SLEs, SIA and DB, reproducing
patterns found in single neurons recorded experimentally dur-
ing seizures.

## 2 Results

Our goal was to construct a biophysical model of a sin-
gle neuron that can reproduce the different firing patterns
recorded when extracellular potassium is increased, while
keeping it sufficiently simple to allow a bifurcation analysis.
The model is schematized in Fig. 1 (see methods section for
the equations). It is a simplification of the classical Hodgkin-
Huxley formalism, which also includes close surround-
ing environment with three compartments (external bath,
extracellular space, intracellular space). The corresponding code is available on github at https://github.com/ddepann/TheModel.

Numerous experiments show that seizures and SLEs are associated with an increase in $[K]_o$ (Fisher et al., 1976; Fröhlich et al., 2008) and that increasing of external [K] can trigger SLEs (S. F. Traynelis & Dingledine, 1988; Stephen F. Traynelis & Dingledine, 1989). The model presented here takes into account the regulation of potassium, via the possible diffusion towards the external bath compartment and its associated potassium concentration $[K]_{bath}$. Changing $[K]_{bath}$ parameter will strongly influence the regulation of extracellular potassium by allowing or not the removal of excess potassium from the extracellular compartment. When $[K]_{bath}$ is low, the bath compartment can pump out the extracellular potassium; but it fails to do so when it is saturated by

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Diagram of characteristics and mechanisms described by the model. Three compartments are represented. A passive diffusion of potassium exists between the external bath and the extracellular space. $Na^+$, $K^+$ and $Cl^-$ ions can be exchanged between the extracellular and intracellular compartments via the Na/K-pump and voltage-gated channels. This model can reproduce the typical patterns of the membrane potential $V_m$, shown in the bottom left subplot, including, from top to bottom, spike train, tonic firing, bursting, seizure like events (SLE), sustained ictal activity (SIA) and depolarization block (DB)</figcaption></figure>

potassium. We thus explored the response of the model as the concentration of $[K]_{bath}$ was increased. The gradual increase in potassium led to 7 sequential qualitative firing patterns: Resting State (RS), Spike Train (ST), Tonic Spiking (TS), Bursting, Seizure-like events (SLE), Sustained Ictal Activity (SIA), and Depolarization Block (DB) (Note that what is called here Spike Train also corresponds to another type of burster from a dynamical point of view (E. Izhikevich, 2007; Saggio et al., 2017)). The corresponding changes of membrane potential for all these patterns are shown in Fig. 2.

The number of firing patterns is higher than in the original Hodgkin-Huxley model. This is due to the fact that the model takes into account the variations of concentration, as evidenced by the variation of the Nernst potentials. The changes in Nernst potential for sodium and potassium ion species are shown in Fig. 2. The simulations are initialized with values observed in a "healthy" resting situation. In some cases, the Nernst potentials display a transient change before reaching a sustained low amplitude oscillations following action potentials, as observed during RS, TS, SIA and DB. During periodic events, (ST, Bursting, SLE), larger oscillations are observed in Nernst potentials. These oscillations are directly linked to the observed oscillations in the slow variables of the model (Eq. (3) and Eq. (4)) describing concentration changes. The rate of oscillation of the slow variables thus explains the duration of periodic events, in line with the assumed essential role of ionic homeostatic regulation.

Each of the firing patterns can be associated to a different behavior, observable experimentally at different scales. The correspondence is established on the basis of their shape and their order of appearance as $[K]_{bath}$ is increased. Tonic and bursting patterns are prototypical. We consider the activity shown on Fig. 2d as SLE at the neuronal scale, as it is similar to the activity typically recorded in individual neurons (Haglund & Schwartzkroin, 1990), in particular the transient episode of depolarization block, in different experimental preparations during SLEs at the network scale (e.g. Figure 6 in (Uva et al., 2013); Fig. 1 in (Bikson et al., 2003) or Fig. 8 in (Jirsa et al., 2014)). Although it is possible to generate SIA *in vitro* (Quilichini et al., 2002), neurons have not been recorded in these conditions. However, the sustained firing pattern in our model cell resembles the regular field activity recorded during SIA *in vitro* (Quilichini et al., 2002). The sustained DB at the single cell level corresponds to what is observed experimentally during network spreading depolarization when $[K]_o$ reaches high levels (Somjen, 2001).

Increasing $[K]_{bath}$ leads to different regimes of variation of external potassium (Fig. 3). These different regimes are associated with a specific dynamic (i.e. type of bifurcation) of the excitability of the membrane. It is therefore possible to link the membrane potential to the variations in extracellular potassium, because of exchanges existing between
<figure><img src="image_5.png" /></figure>

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

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

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

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

<figure><img src="image_10.png" /><figcaption>**Fig. 2** Qualitative mode of behavior of the membrane potential and Nernst potentials. In blue: time series of the membrane potential $V_m$ for the following patterns of activity: (a) Spike train, (b) Tonic spiking (TS), (c) Bursting, (d) Seizure-like Event (SLE), (e) Sustained Ictal Activity (SIA), (f) Depolarization Block (DB). In red: Nernst potential of sodium, in green: Nernst potential of potassium with their specific Y axis on the right side of the panels. If the value of</figcaption></figure>

$[K]_{\text{bath}}$ stays below 6 mM, the system remains in resting state around -72 mV. Specific patterns of activities start to appear with a diminution of the Nernst potential of sodium and an increase of the Nernst potential of potassium. When periodic events are occurring (panels c and d), oscillations can also be observed in the Nernst potential of both ions
compartments (i.e. via the slow variable), as shown in Fig. 4. In the next subsection, we detail these dynamical interactions for the different patterns of activity, following the order of appearance when $[K]_{\text{bath}}$ increases.

## 2.1 Resting states, spike train and tonic spiking in low $[K]_{\text{bath}}$

Resting state is found when $[K]_{\text{bath}}$ is around the normal value of $[K]_o$ (called $[K]_{0,o}$, see Methods section). If $[K]_{\text{bath}}$ is smaller than $[K]_{0,o}$ the membrane potential slowly hyperpolarizes, due to a diffusion of potassium in the direction of the external bath. When $[K]_{\text{bath}}$ slightly increases (> 7 mM), ST appears through a SNIC bifurcation. The offset is also a SNIC bifurcation. In this case, the onset and offset bifurcations can be easily identified by their characteristic features (Saggio et al., 2017), and confirmed by numerical methods (using SymPy (Meurer et al., 2017) and SciPy (Millman & Aivazis, 2011) libraries). With higher value of $[K]_{\text{bath}}$ (> 8 mM), TS occurs. In this condition, $[K]_o$ stabilizes (Fig. 3), and the neuron fires at a constant frequency (Fig. 2b). The occurrence of regular spiking due to an increase of $[K]_o$ through diffusion from the bath is consistent with experiments (El Houssaini et al., 2015; Strauss et al., 2008).

## 2.2 Bursting and seizure-like events

Bi-stable behavior occurs when the slow system starts to oscillate when $[K]_{\text{bath}}$ is further increased. The model (with parameters listed in Table 1) displays bursting and SLE, successively. Bursts are square-wave bursts (SN/Homoclinic bifurcations) and SLEs also show SN and Homoclinic

<table><caption>Table 1 Parameters values</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Membrane capacitance</td><td>C<sub>m</sub></td><td>1 nF</td></tr><tr><td>Gating time constant</td><td>τ<sub>n</sub></td><td>0.25 ms</td></tr><tr><td>Chloride conductance</td><td>g<sub>Cl</sub></td><td>7.5 nS</td></tr><tr><td>Maximal potassium conductance</td><td>g<sub>K</sub></td><td>22 nS</td></tr><tr><td>Maximal sodium conductance</td><td>g<sub>Na</sub></td><td>40 nS</td></tr><tr><td>Potassium leak conductance</td><td>g<sub>K,l</sub></td><td>0.12 nS</td></tr><tr><td>Sodium leak conductance</td><td>g<sub>Na,l</sub></td><td>0.02 nS</td></tr><tr><td>Intracellular volume</td><td>ω<sub>i</sub></td><td>2160 µm<sup>3</sup></td></tr><tr><td>Extracellular volume</td><td>ω<sub>o</sub></td><td>720 µm<sup>3</sup></td></tr><tr><td>Intra/extra cellular volume ratio</td><td>β = ω<sub>i</sub>/ω<sub>o</sub></td><td>3</td></tr><tr><td>Conversion factor</td><td>γ</td><td>0.04 mol/C·µm<sup>3</sup></td></tr><tr><td>Diffusion time constant</td><td>ε</td><td>0.001 ms<sup>-1</sup></td></tr><tr><td>Maximal Na/K pump current</td><td>ρ</td><td>250 pA</td></tr></tbody></table>

bifurcations at onset and offset, respectively (See supplementary information: S1, S2, S3, S4). Here, the slow subsystem oscillates in a self-sustained manner (Fig. 4a-f), generating recurrent bursting or SLEs, with important variations of $[K]_o$, due to oscillations in the slow subsystem. The combined effects of oscillations of $\Delta[K]_I$ and $[K]_g$ explain the changes in the Nernst potential of potassium (and sodium, which is linked to potassium in the model), thus changing neural excitability. During spiking activity, voltage-gated potassium channels open increasing potassium current $I_K$. The influence of $I_K$ in the equation of $\Delta[K]_I$ (Eq. 3), explains the decrease of $\Delta[K]_I$, hence the increase of $[K]_o$ through equations (Eq. 16) and (Eq. 20). This is consistent with the observations described in (Fisher et al., 1976). The increase in $[K]_o$ starts with the occurrence of burst and SLEs. Thus,

**Fig. 3** Variation of extracellular potassium concentration as a function of $[K]_{\text{bath}}$. Minimal and maximal external potassium $[K]_o$ and mean (dash line) concentration observed during simulations done for different values of the parameter $[K]_{\text{bath}}$. Due to diffusion from the external bath, increasing $[K]_{\text{bath}}$ leads to variations in $[K]_o$. Different patterns are observed for each range of $[K]_{\text{bath}}$: (a) resting state, (b) spike train, (c) regular spiking, (d) burst, (e) seizure like event, (f) sustained ictal activity, (g) depolarization block. The periodic events (spike train, burst and seizure-like event) correspond to the range of $[K]_{\text{bath}}$ where $[K]_o$ periodically oscillates

<figure><img src="image_11.png" /></figure>
<figure><img src="image_12.png" /></figure>

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

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

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

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

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

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

<figure><img src="image_19.png" /></figure>
►Fig. 4 Time series of membrane potential, $\Delta[K]_i$, $[K]_g$, and $[K]_o$. Numerical integration, with x-axis in millisecond. (a) spike train with $[K]_{bath}=7.5$ mM, (b) spike train with $[K]_{bath}=7.5$ mM and $\gamma=0.04$, $\varepsilon=0.002$, (c) Burst with $[K]_{bath}=12.5$ mM, (d) Burst with $[K]_{bath}=12.5$ mM, and $\gamma=0.06$, $\varepsilon=0.002$, (e) SLE with $[K]_{bath}=16$ mM, (f) SLE with $[K]_{bath}=16$ mM and $\gamma=0.08$, $\varepsilon=0.0008$, (g) SIA with $[K]_{bath}=17.5$ mM, (h) DB with $[K]_{bath}=20$ mM. If not specified, the parameter values used here are the reference parameters described in the method section. Variations of $\Delta[K]_i$ and $[K]_g$ induce different patterns of activity. The combined effects lead to the observed variations in $[K]_o$. The time scale of the slow variables $\gamma$ and $\varepsilon$ influence the shape of $V_m$ allowing the system to exhibit SN or SNIC bifurcation at the onset of the events

it is not the cause of the event but a consequence of homeostasis dysregulation (i.e. augmentation of $[K]_{bath}$).

## 2.3 Steady states, SIA events and DB, in high $[K]_{bath}$ conditions

SIA events (Fig. 2g) appears for $[K]_{bath}$ around 17.5 mM (Figs. 3 and 4), i.e. above the threshold value for SLEs as reported experimentally (El Houssaini et al., 2015). If no other mechanisms act to stop it, these oscillations remain constant (analogous to refractory status epilepticus). Permanent DB occurs for even higher values of $[K]_{bath}$ (> 18.0 mM, Fig. 3) as also reported experimentally (El Houssaini et al., 2015). In these cases, after a peak value (Fig. 4h), $[K]_o$ stabilizes, explaining the short range of variation (Fig. 3). These steady-states start like a SLEs (Fig. 4e-f), then the slow variables stabilize and $[K]_o$ remains constant at a high value (Fig. 4g-h).

We conclude that the model reproduces all the transitions between resting state, spike train, regular spiking, burst, seizure like event, sustained ictal activity, and depolarization block as seen experimentally (El Houssaini et al., 2015) when external potassium increases. In particular, the single cell model reproduces the experimental behavior of neurons recorded in networks generating such activities (El Houssaini et al., 2015).

These simulations were obtained when using a “healthy” situation, corresponding to a neuron recorded in a non-pathological context. In epilepsy, the regulatory mechanisms of neuronal homeostasis are affected (Boison et al., 2013; McDonald et al., 2018; Zilberter & Zilberter, 2017). In the next section, we model a "pathological" situation to obtain insights of what might happen in chronic epilepsy.

## 2.4 Analysis of the model in a pathological context

Glial cells normally ensure the regulation of the extracellular concentration of $K^+$ (Coulter & Steinhäuser, 2015; Kofuji & Newman, 2004; Olsen et al., 2015; Walz, 2000), which is impaired in epilepsy (Coulter & Steinhäuser, 2015; Hubbard

& Binder, 2016; Rangroo Thrane et al., 2013; Scholl et al., 2009). Our model does not include glial cells but it is sensitive to $[K+]_o$, the parameter that glial cells control. To model a pathological context characterized by glial cell dysfunction without making the model more complex, we approximate potassium buffering by its diffusion between the extracellular compartment and the bath, varying the parameter $\varepsilon$. Homeostasis of intracellular ions is also important to consider. The parameter $\gamma$ can be considered biophysically as a conversion factor, but also phenomenologically as the parameter that links the evolution of ion concentration with the activity of the membrane. Thus, the impairment of mechanisms not included in the model, such as co-transporters and exchangers (Hille, 2001; Kandel et al., 1981), which will affect the evolution of ion concentration, can be approximated phenomenologically by changes of $\gamma$. Varying the time constants of the slow subsystem ($\varepsilon$ and $\gamma$), leads to different bi-stable behaviors. Two examples are shown in Fig. 4(b) with $\gamma=0.04$, $\varepsilon=0.002$, (d) $\gamma=0.06$, $\varepsilon=0.002$, and (f) $\gamma=0.08$, $\varepsilon=0.0008$. In these examples, potassium concentration oscillations are affected leading to a change in the duration of the events. For burst and SLE shown in Fig. 4. d and f, the model exhibits a different class of onset bifurcation. For both, a saddle-node on invariant cycle (SNIC) bifurcation at the onset and homoclinic bifurcation at the offset can be identified, based on their specific dynamics and resulting shapes (E. Izhikevich, 2007; Saggio et al., 2017).

The other key parameter to consider is the pump rate $\rho$. The Na/K-ATPase is described by Eq. (8) in the model. In a biological neuron, the pump depends on ATP and during status epilepticus, the ATP concentration increases due to high needs and then decreases (Lietsche et al., 2016). The ATP concentration is not taken into account in the model, but the maximal Na/K-pump rate is modulated by the parameter $\rho$. This parameter also influences the shape of $I_{pump}$ response as a function of $[Na]_i$ and $[K]_o$ (Fig. 5a). For large values of $\rho$, the pump is activated for lower value of $[Na]_i$ and $[K]_o$ (Fig. 5a). We find that burst duration changes with $\rho$ for a fixed $[K]_{bath}$ (Fig. 5b), where a faster activation (higher $\rho$) leads to shorter bursts. The augmentation of $\rho$ does not necessary lead to an increase of $I_{pump}$; it affects the general dynamics of the whole system (Fig. 5c).

Thus, changes related to the Na/K-ATPase affect mainly the duration of the events, while impairment of mechanisms related to the regulation of $K^+$ concentration affects the type of pattern (i.e. onset/offset bifurcation types).

The biophysical model is able to reproduce general patterns of activities (i.e. periodic events) as generated by the phenomenological model (El Houssaini et al., 2015). However, it still contains too many parameters for an exhaustive study of its possible dynamics. With a reduced number of variables, we can make a detailed study of the dynamics for a given set of parameters.
## 2.5 Dynamics of the model

The previous sections describe biophysical aspects in relation to biological observations. In this last section we show how the model can make the link with the theoretical framework. Since the biophysical model contains few differential equations, it is possible to use the tools of dynamical systems theory to directly compare its behavior with the generic model.

The model can be divided into the fast (V, n, respectively Eq. 1 and Eq. 2) and the slow subsystems ($Δ[K]_i$, $[K]_g$, respectively Eq. 3 and Eq. 3). The slow system can oscillate and drive the fast system between different behaviors, in particular switching between resting state and fast oscillations to obtain bursting-like activity. In this subsection, we call burster a system allowing these periodic events. To create the oscillation in the slow subsystem, theoretical works show that two mechanisms are possible (E. Izhikevich, 2007; Saggio et al., 2017): Slow-Wave (SW) burster, where the slow subsystem is made of two equations, independent of the fast system, or Hysteresis-Loop (HL) burster where the slow subsystem is made of only one equation that depends on the fast system. Each has typical onset/offset bifurcation pairs. These specific paths for bursting have been identified in the generic model (Saggio et al., 2017), and are reproduced in Fig. 6a. We first verified if the relations between the equations of the slow and fast systems allow the existence of the mechanisms described previously. Because $I_K$ (Eq. (6))

**Fig. 6** Comparison with the generic model. (a) Paths for bursting activity of the generic model proposed by Saggio et al. adapted with authorization from (Saggio et al., 2017), for hysteresis-loop burster (left) and slow-wave burster (right), the yellow star corresponds to the center of the region captured by our model. (b) Bifurcation diagram of our model, where the white area corresponds to 'resting state only' region, the dark red corresponds to a depolarized region, and the light-red region is the region of bi-stability. The yellow star corresponds to the point also found in the generic model, where the SH, SNIC and SN bifurcations intersect. In the top diagram, the green line corresponds to the path taken by the burster, in the bottom one to the path taken by the SLE. (c) Classes of bursters found in the model, and the corresponding path in the generic model

depends on V and n, the Eq. (3) depends on the fast system. This corresponds to a relation that exists in an HL burster. The second equation of the slow subsystem, Eq. (4), also depends on the Eq. (3), through the Eq. (20). Thus, there exists a relation between the two equations of the slow system, enabling oscillation such as in a SW burster. These relations between the variables of our model allow obtaining the two types of bursters previously described.

We therefore tested for possible correspondences between our model and the generic model. We were able to identify the regions in the generic model capturing the dynamics reproduced by our model in Fig. 6a. The center of the region of interest has been marked with a yellow star in Fig. 6a. for the generic model and its correspondence in the bifurcation diagram of our model in Fig. 6b. In this bifurcation diagram

<figure><img src="image_20.png" /><figcaption>**Fig. 5** Influence of the activity of the Na/K-pump. (a) $I_{\text{pump}}$ function for $\rho=25$ (green), $\rho=1000$ (orange), $\rho=2500$ (blue). The initial slope when the system moves away from the concentrations at rest is affected, explaining the modification of the influence of $I_{\text{pump}}$ in the dynamic of the system. (b) Burst duration as a function of $\rho$ for $[K]_{\text{bath}}=14.0$ mM. Bursts have shorter durations for higher value of $\rho$. (c) Minimal and maximal pump current, $I_{\text{pump}}$, observed during simulation done with $[K]_{\text{bath}}=14.0$ mM. The range of $I_{\text{pump}}$ decreases for higher $\rho$ values</figcaption></figure>

<figure><img src="image_21.png" /></figure>
<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c)</figcaption></figure>

<table><thead><tr><td>Onset/Offset bifurcations</td><td>corresponding path in Saggio et al.</td><td>time series of V<sub>m</sub></td></tr></thead><tbody><tr><td>Saddle-Node/Homoclinic</td><td>c2</td><td></td></tr><tr><td>sup. Hopf/Homoclinic</td><td>c10</td><td></td></tr><tr><td>Saddle-Node/sup. Hopf</td><td>c3</td><td></td></tr><tr><td>SNIC/Homoclinic</td><td>c6</td><td></td></tr><tr><td>SNIC/sup. Hopf</td><td>c7</td><td></td></tr><tr><td>SNIC/SNIC</td><td>c5</td><td></td></tr></tbody></table>

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

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

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

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

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

<figure><img src="image_29.png" /></figure>
<table><caption>**Table 2** Physiological reference values</caption><thead><tr><th></th><th>Ion</th><th>Concentration</th></tr></thead><tbody><tr><td>External bath</td><td>[K]<sub>bath</sub></td><td>[2–30] mM</td></tr><tr><td rowspan="3">Extracellular</td><td>[K]<sub>0,o</sub></td><td>4.8 mM</td></tr><tr><td>[Na]<sub>0,o</sub></td><td>138 mM</td></tr><tr><td>[Cl]<sub>0,o</sub></td><td>112 mM</td></tr><tr><td rowspan="3">Intracellular</td><td>[K]<sub>0,i</sub></td><td>140 mM</td></tr><tr><td>[Na]<sub>0,i</sub></td><td>16 mM</td></tr><tr><td>[Cl]<sub>0,i</sub></td><td>5 mM</td></tr></tbody></table>

we show two possible paths of our model, for burst behavior (Fig. 6b, top) and for SLE (Fig. 6b, bottom). It crosses regions of stable resting state (in white), depolarized (red), and bistable (light red). It is therefore possible to establish a non-exhaustive list of the correspondences between the paths of the two models. The paths for the periodic events have been listed in Fig. 6c. The spike train, Bursting and SLE behaviors correspond to paths, c5, c2 and c10, respectively. The bursting behavior with changes in ε and γ (Fig. 4b) that represents the SNIC/SH bifurcation corresponds to the path c6. The model proposed here, consistent with biophysics, fits into the framework of the generic model.

Since our biophysical model reproduces the bifurcations of the generic model for different types of network activities, it becomes possible to investigate the ionic mechanisms underlying the onset/offset bifurcations. The fast subsystem can be described fixing all parameters (Tables 1 and 2) and considering the two slow variables as parameters. Fixed points can thus be found for different values of Δ[K]ᵢ and [K]ᵣ as shown in Fig. 7. Importantly, some parameter values allow a bi-stable behavior. It is thus possible to understand

the direct relationship between the biophysical variations in potassium concentration and the type of bifurcations by observing the trajectory of the membrane potential in this space for periodic events identified previously. During periodic oscillatory behavior, the neuron is initially in resting state (blue plane). The membrane potential slowly increases due to the rise in extracellular potassium, until it reaches a SN (green plane) and then encounters a limit cycle. The slow subsystem then drives it to a negative value of Δ[K]ᵢ, were the limit cycle meets a SN producing homoclinic bifurcation. These bifurcations are observed at the onset and offset of bursting and SLE behaviors in the model. To have a better understanding of these trajectories, animations with the dynamics of the fast subsystem are available in supplementary material (Figs. S1, S2, S3, S4). We therefore have here a means of bringing together the biophysical aspects, described previously, with the phenomenological vision of dynamical systems approach.

## 3 Discussion

The aim of this work is to develop a minimal biophysical model at single neuron level based on time scale separation, where the system is able to reproduce the dynamics which have been identified in experiments (Bikson et al., 2003; Jirsa et al., 2014; Quilichini et al., 2002; Somjen, 2001; Uva et al., 2013) and described by generic models (Jirsa et al., 2014; Saggio et al., 2017). For this purpose, we developed a three-compartment model: a cell equipped with voltage-gated channels to generate action potentials, and Na⁺/K⁺ pump to maintain stable ion concentration, an extracellular space surrounding the cell and an external

<figure><img src="image_30.png" /><figcaption>**Fig. 7** Fixed points of the fast subsystem. Fixed point of the fast subsystem (V<sub>m</sub>) considering the variables from the slow subsystem as parameters. We used a numerical methods with SymPy (Meurer et al., 2017) and SciPy (Millman & Aivazis, 2011) libraries, to find the roots and the eigenvalues of the Jacobians of the 2D fast subsystem,</figcaption></figure>

and thus the stability considering the existence and the sign of real and imaginary parts of the eigenvalues of the Jacobians. Blue: stable node, green: saddle node, cyan: stable focus, magenta: unstable focus, red: unstable node. Two different angles of view are presented, illustrating the manifold that permits bi-stability
bath that can uptake/release potassium from/to extracellular space. We managed to describe the interaction between these compartments using a system of four differential equations describing two fast and two slow variables. The fast variables delineate excitability while the slow ones, outline potassium changes from the first and third compartments. The sodium concentration changes are not excluded from our model but are linked to potassium through the electroneutrality principle. We have shown that despite its simplicity the model was able to mimic six electrophysiological behaviors classically recorded in neurons and neuronal networks, via the variation of only one parameter. All parameter values were within biophysical ranges (Table1) (Hille, 2001; Kandel et al., 1981). Recent experimental work, using selective ionic pumps to deliver locally $K^+$, confirmed observations from the model, the elevation of $[K^+]$ modified the spiking profile of a bursting neuron in the hippocampal slice (Arbring Sjöström et al., 2021). However, the model has two main limitations. The fast system describes only intrinsic excitability and does not include synaptic currents. The slow system is based (only) on potassium concentration. Introducing synaptic inputs would increase the dimension of the system. We propose that synaptic inputs would act as a noise generator increasing the probability to reach the bifurcation as demonstrated experimentally (Jirsa et al., 2014); including them should not change the general behavior of the model. Furthermore, ion homeostasis is not reduced solely to potassium. Potassium is just one candidate among many others for the slow system. Numerous studies have reported large changes in concentration of $Ca^{2+}$ (Heinemann et al., 1986), $Cl^-$ (Miles et al., 2012; Raimondo et al., 2015) and neurotransmitters during seizures (Chapman et al., 1984; During & Spencer, 1993). Likewise, decreasing extracellular $Ca^{2+}$ leads to seizures (Jefferys & Haas, 1982), which are characterized by SN/homoclinic bifurcations (Jirsa et al., 2014). Since it is possible to trigger similar SLEs via totally different biophysical mechanisms (Jirsa et al., 2014), we propose that the $K^+$-dependent mechanism we describe, is one among many the possible paths leading to the same end point. In our model, changes in potassium constitute the causal factor driving the neuron through different types of activities. Although similar changes in potassium are measured experimentally when networks (and not cells) undergo such transitions, causality has not been demonstrated experimentally, only correlation. Another limitation exists due to the formalism used. If $[K]_{bath}$ tends to zero then membrane potential hyperpolarize until the Nernst potential are is longer defined due to a division by zero. We reach here the limit of the conductance-based model from Hodgkin-Huxley formalism. Due to the expression of the Nernst potential, if the ratio $[K]_o/[K]_i$ approaches zero, then the $I_K$ current increases towards infinity, which is not

physiologically plausible. Another factor to consider is that the dynamics of the single cell is driven by slow changes of extracellular variables, which, in a biological system, is shared with neighboring cells. So, these slow variables can also be responsible for the genesis of network activity (Naze et al., 2015). As these mechanisms exist both at the network and single neuron level, it would be simplistic to conclude that a seizure at the network level is due to the combined expression of seizures at the single cell level. Since a neuronal network can be seen as a complex system of many components, coupled in a non-linear manner, seizures may just be an emergent property, perhaps taking advantage of the fact that they are already encoded at the single cell level. The same consideration applies to other pathological activities such as SIA and DB, which corresponding pattern have been found in dynamics of our model.

We only studied the dynamics for variations of few chosen parameters based on physiological observations identified in previous works. The parameters explored here show that the model can produce different combinations of onset/ offset bifurcations. Numerous studies used ion concentration variations in biophysical models to generate various types of activity (Barreto & Cressman, 2011; Bernard et al., 2014; Cressman et al., 2009; Florence et al., 2009; Krishnan et al., 2015; Øyehaug et al., 2012; Wei et al., 2014a, b). Descriptions of ion concentration dynamics for bursting have been done by Barreto et al. (Barreto & Cressman, 2011), based on a slow/fast system. In this work, the bifurcations for SLEs are SNIC and Hopf. This approach, based on ion concentration dynamics, permits the unification of spike, seizure and spreading depression proposed by Wei et al. (2014a, b). As different models can lead to similar dynamics (Prinz et al., 2004), this suggests that different minimalist models are possible to obtain a unified framework. In comparison to previous works (Barreto & Cressman, 2011; Cressman et al., 2009), our model does not take into account $[Ca]_2^+$, and includes a constant leak current for $Na^+$ and $K^+$. We reduced the fast subsystem to only two equations. Also, we use only differential equation for the evolution $[K]^+$ and not for $[Na]^+$, we consider the evolution through the interdependence between both. Although a number of similar biophysical elements are taken into account, the system of equations obtained is different. We have a different model and so, a different structure of the phase space. This structural difference is important because it explains the different dynamics that the model can reproduce. This explains why we get a different repertoire in terms of types of bifurcations.

Here, we propose a conductance-based model of the neuronal membrane, exhibiting an extended repertoire of behavior and introducing sustained ictal activity in a unified framework. Another difference with previous works is that our model can exhibit bi-stable modes saddle-node/
homoclinic bifurcations, which are the most commonly observed in recordings from patients and experimental animal models (Jirsa et al., 2014). Our model does not take into account variation of volume or oxygen homeostasis as in (Wei et al., 2014a, b) but, only variations of ion concentrations, driven by diffusion of potassium from the external bath. It seems intuitive that other biological variables could be considered as slow variables to drive the fast subsystem in a reduced biophysical model. The work of Øyehaug et al. (2012) presents interesting dynamical features with saddle-node/homoclinic bifurcations for SLEs. However, this model is much more complex as it describes numerous biological features and mechanisms. In comparison to previous works (Barreto & Cressman, 2011; Krishnan et al., 2015; Øyehaug et al., 2012; Wei et al., 2014a, b), our model is reduced to only four equations. We sought to include only a minimal number of mechanisms necessary to reproduce neural dynamics. Chizhov et al. (2018) proposed a biophysical model (Epileptor-2) of ictal activities based on the Epileptor (Jirsa et al., 2014), using different differential equations. In high potassium conditions, Epileptor-2 produces bursts of bursts, described as ictal-like discharges. However, the most common form of seizure belongs to the saddle-node/homoclinic form, which starts with low voltage fast activity, and ends with bursts slowing down in a logarithmic fashion. The latter was reproduced in the present model, including the period during which neurons stop firing (depolarization block) after seizure onset. Another difference lies in K⁺ dynamics. In Epileptor-2, neuronal firing ends when extracellular K⁺ returns to baseline level (Fig. 10 in Chizhov et al., 2018), whereas in the present model, there is a delay, as consistently found experimentally, as a result of glial cell action. This phenomenon in our model can be visualized by observing the evolution of [K]₀ in Fig. 4.

In conclusion, we developed a biophysical model of a single neuron that, despite its simplicity, is able to generate, in a unified framework, many patterns of neuronal network activity found in experimental recording as well as in generic mathematical models. We show that transition from physiological to paroxysmal activity can be obtained by variation of model parameters relating to ion homeostasis while excitability parameters remained constant. Thus, we proposed a simple biophysical model comparable to generic models (El Houssaini et al., 2015; Jirsa et al., 2014; Saggio et al., 2017), offering the possibility of a biological interpretation of observed dynamics. Neuronal networks increase in complexity from flies to humans, but the basic properties of neurons are roughly conserved. The present study shows that acting on an external variable allows single neurons to go through various patterns of activities, which are also found at the network level in the form of seizures, sustained ictal activity and depolarization block (Cunliffe et al., 2015; Jirsa et al., 2014). We propose that they constitute one of the most

primitive forms of activities, appearing as soon as neurons are present.

## 4 Materials and methods

In this project we aim to build a minimal biophysical model that describes different electrophysiological states of a single neuron, the model is schematized in Fig. 1. The model describes three compartments: the intracellular space (ICS), the extracellular (ECS) space and the external bath (EB). Parameters chosen correspond to values observed in whole cell recording. The ion exchange between the ICS and the ECS is carried out by the current flowing through the sodium, potassium, and chloride voltage-gated channels (Eq. (5), (6) and (7)), and by the sodium-potassium pump generated current (Eq. (8)). Parameters values for these currents have identified in (Hamada et al., 2003; Hille, 2001; Läuger, 1991) and the membrane capacitance in (Golowasch et al., 2009). Passive diffusion of potassium exists (Eq. (4)), between EB and ECS. The EB is mimicking the K⁺ buffering of vasculature/astrocytes. In ICS and ECS actualization of potassium and sodium concentrations are done (Eq. (14)-(20)). The γ parameter has the same unit as the inverse of the Faraday constant, and it is a scaling parameter that permit to include all the mechanisms not detailed in this model which affect the concentration variations (such as co-transporter, exchangers). The values of all the parameters used are given in Table 1 and physiological reference and initial values are given in Table 2 and Table 3.

The model is a slow-fast dynamical system based on 4 equations. The fast system describes the membrane potential Eq. (1) and potassium conductance gating variable Eq. (2). The slow system describes intracellular potassium concentration variation Eq. (3) and extracellular potassium buffering by external bath Eq. (4).

$$ \frac{dV}{dt} = -\frac{1}{C_m}(I_{Cl} + I_{Na} + I_K + I_{pump}) \quad (1) $$

$$ \frac{dn}{dt} = \frac{n_{\infty}(V) - n}{\tau_n} \quad (2) $$

$$ \frac{d\Delta[K]_i}{dt} = -\frac{\gamma}{\omega_i}(I_K - 2I_{pump}) \quad (3) $$

$$ \frac{d[K]_g}{dt} = \varepsilon([K]_{bath} - [K]_o) \quad (4) $$

With currents:

$$ I_{Na} = (g_{Na,l} + g_{Na}m_{\infty}(V)h(n))(V - 26.64\log(\frac{[Na]_o}{[Na]_i})) \quad (5) $$
<table><caption>**Table 3** Initial values</caption><thead><tr><th>Variable</th><th>Initial value</th></tr></thead><tbody><tr><td>[K]<sub>o</sub></td><td>[K]<sub>0,o</sub></td></tr><tr><td>[Na]<sub>o</sub></td><td>[Na]<sub>0,o</sub></td></tr><tr><td>[Cl]<sub>o</sub></td><td>[Cl]<sub>0,o</sub></td></tr><tr><td>[K]<sub>i</sub></td><td>[K]<sub>0,i</sub></td></tr><tr><td>[Na]<sub>i</sub></td><td>[Na]<sub>0,i</sub></td></tr><tr><td>[Cl]<sub>i</sub></td><td>[Cl]<sub>0,i</sub></td></tr><tr><td>Δ[K]<sub>i</sub></td><td>0</td></tr><tr><td>[Kg]</td><td>0</td></tr><tr><td>V</td><td>-70 mV</td></tr><tr><td>n</td><td>n<sub>∞</sub>(-70)</td></tr></tbody></table>

$$I_K = (g_{K,l} + g_K n)(V - 26.64 \log(\frac{[K]_o}{[K]_i})) \qquad (6)$$

$$I_{Cl} = g_{Cl}(V + 26.64 \log(\frac{[Cl]_o}{[Cl]_i})) \qquad (7)$$

$$I_{pump} = \rho \frac{1}{1 + \exp(\frac{1}{2}(21 - [Na]_i))} \frac{1}{1 + \exp(5.5 - [K]_o)} \qquad (8)$$

And conductance variables:

$$n_{\infty}(V) = \frac{1}{1 + \exp(\frac{1}{18}(-19 - V))} \qquad (9)$$

$$m_{\infty}(V) = \frac{1}{1 + \exp(\frac{1}{12}(-24 - V))} \qquad (10)$$

$$h(n) = 1.1 - \frac{1}{1 + \exp(-8(n - 0.4))} \qquad (11)$$

<figure><img src="image_31.png" /><figcaption>**Fig. 8** modification in gating variables. (a) $n_{\text{inf}}$ of our model in blue, and $n_{\text{inf}}$ and $1/\tau$ of the Hodgkin-Huxley model respectively in dash blue and red, function of the membrane potential. (b) Response of</figcaption></figure>

The fast subsystem of the model, (Eq. (1)&(2)), is a reduction and simplification of conductance-based models, first described by Hodgkin-Huxley (HH). From the original publication (Hodgkin & Huxley, 1952) the activation variable of K+ channels is determined by the equation (Eq. 12):

$$\frac{dn}{dt} = \alpha_n(1 - n) - \beta_n n \qquad (12)$$

where $\beta(V)$ and $\alpha(V)$ are the voltage-dependent rate constants determining the probability of transitions between, respectively, opened and closed state of the ion channel. To simplify the model, we propose to describe the variable n, through the voltage-dependent parameter $n_{\text{inf}}(V)$ and a constant parameter $\tau_n$. In our model, $n_{\text{inf}}(V)$ is the probability to find a channel at open state at a given membrane potential while $\tau_n$ is the fixed time constant that described the speed for channels to respond to the change of membrane potential. Based on available data in the literature (Bekkers, 2000; Hodgkin & Huxley, 1952), and considering that the mean number of channels opened at a given potential is constant, we could qualitatively estimate this relationship (Eq. 9). In the HH model, the time constant is dependent on the membrane potential due to the formalism used (Eq. 12). The HH model was constructed using experiments performed on the giant squid axon, which differ from mammalian neurons. We compare the $n_{\text{inf}}(V)$ of our model and $1/\tau(V)$, and $n_{\text{inf}}(V)$ of the HH model in Fig. 8(a). The shape has been kept from the HH model but starts to increase for lower values of membrane potential. For the voltage-gated sodium channels, variables for opening, m, and for closing, h, have been described (Hodgkin & Huxley, 1952). With the same logic, we can consider the percentage of all population of channels opened. But because this is a very fast mechanism (Hille, 2001), it can be considered as an instantaneous function of V (E. Izhikevich, 2007) (Eq. 10). Krinskii and Kokoz

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

the fast subsystem of our model to step current stimulation (red) with three different values of $\tau$ (0.1, 0.25, 0.5 ms). The value of $\tau$ influence the frequency rate spike for a same injected current
(Krinskii & Kokoz, 1973) showed that n(t) + h(t) is almost constant, so h can be considered as a function of n. Because of the previous modification, we adapted this fitting to obtain the equation of h(n) (Eq. 11). Due to these simplifications, the interdependence of gating variables makes the spiking rate dependent on τ, as shown in Fig. 8(b).

To be able to take into account concentration variation limiting the number of equations we applied reductions. Inspired by the work of Hübel (Hübel, 2015; Hübel & Dahlem, 2014), electroneutrality permits the Eq. (13), and so to the Eq. (14). The ratio (C_m γ)/ω_i is very small (<10⁻⁵) and so, the right-hand side of Eq. (14) could be considered to be zero. The chloride concentration changes are assumed to be small and regulated by mechanisms which are not described in our model (Doyon et al., 2016). So, in our model, the chloride concentration remains constant.

$$ \frac{dV}{dt} = \frac{\omega_i}{C_m \gamma} \left( \frac{d\Delta[K]_i}{dt} + \frac{d\Delta[Na]_i}{dt} + \frac{d\Delta[Cl]_i}{dt} \right) \quad (13) $$

$$ \Delta[K]_i + \Delta[N\alpha]_i + \Delta[Cl]_i = \frac{C_m \gamma}{W_i} (V - V_0) \quad (14) $$

Thanks to these reductions, concentration variations are calculated as follow:

$$ \Delta[N\alpha]_i = -\Delta[K]_i \quad (15) $$

$$ \Delta[N\alpha]_o = -\beta\Delta[N\alpha]_i \quad (16) $$

$$ \Delta[K]_o = -\beta\Delta[K]_i \quad (17) $$

$$ [K]_i = [K]_{0,i} + \Delta[K]_i \quad (18) $$

$$ [Na]_i = [Na]_{0,i} + \Delta[Na]_i \quad (19) $$

$$ [Na]_o = [Na]_{0,o} + \Delta[Na]_o \quad (20) $$

$$ [K]_o = [K]_{0,o} + \Delta[K]_o + [K]_g \quad (21) $$

All simulations were obtained thanks to numerical methods using odeint function from SciPy library (Millman & Aivazis, 2011).

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s10827-022-00811-1.

**Funding** This research was partially supported by the European Union's Horizon 2020 research and innovation program under grant agreement No. 785907 (SGA2), and No. 945539 (SGA3) Human Brain Project and was also financed in part by the "Coordenação de Apresentação de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001", the National Research Council "Conselho Nacional de

Desenvolvimento Científico e Tecnológico (CNPq)" and the Family Allowance Fund (CAF de Marseille, France).

## Declarations

**Conflict of interest** The authors declare no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## References

Arbring Sjöström, T., Ivanov, A. I., Bernard, C., Tybrandt, K., Poxson, D. J., Simon, D. T., & Berggren, M. (2021). Design and Operation of Hybrid Microfluidic Iontronic Probes for Regulated Drug Delivery. *Advanced Materials Technologies*, 6(2), 2001006. https://doi.org/10.1002/admt.202001006

Barreto, E., & Cressman, J. R. (2011). Ion concentration dynamics as a mechanism for neuronal bursting. *Journal of Biological Physics*, 37(3), 361–373. https://doi.org/10.1007/s10867-010-9212-6

Bekkers, J. M. (2000). Distribution and activation of voltage-gated potassium channels in cell-attached and outside-out patches from large layer 5 cortical pyramidal neurons of the rat. *The Journal of Physiology*, 525(Pt 3), 611–20. https://doi.org/10.1111/j.1469-7793.2000.T01-2-00611.X

Bernard, C., Naze, S., Proix, T., & Jirsa, V. K. (2014). Modern concepts of seizure modeling. *International Review of Neurobiology*, 114, 121–153. https://doi.org/10.1016/B978-0-12-418693-4-00006-6

Bikson, M., Hahn, P. J., Fox, J. E., & Jefferys, J. G. R. (2003). Depolarization Block of Neurons During Maintenance of Electrographic Seizures. *Journal of Neurophysiology*, 90(4), 2402–2408. https://doi.org/10.1152/jn.00467.2003

Blenkinsop, A., Valentin, A., Richardson, M. P., & Terry, J. R. (2012). The dynamic evolution of focal-onset epilepsies - combining theoretical and clinical observations. *European Journal of Neuroscience*, 36(2), 2188–2200. https://doi.org/10.1111/j.1460-9568.2012.08082.x

Boison, D., Sandau, U. S., Ruskin, D. N., Kawamura, M. J., & Masino, S. A. (2013). Homeostatic control of brain function – new approaches to understand epileptogenesis. *Frontiers in Cellular Neuroscience*, 0(JUNE), 109. https://doi.org/10.3389/FNCELI.2013.00109

Bragin, A., Penttonen, M., & Buzsáki, G. (1997). Termination of epileptic afterdischarge in the hippocampus. *The Journal of neuroscience : the official journal of the Society for Neuroscience*, 17(7), 2567–79. https://doi.org/10.1523/JNEUROSCI.17-07-02567.1997

Chapman, A. G., Westerberg, E., Premachandra, M., & Meldrum, B. S. (1984). Changes in Regional Neurotransmitter Amino Acid Levels in Rat Brain During Seizures Induced by l-Allylglycine, Bicuculline, and Kainic Acid. *Journal of Neurochemistry*, 43(1), 62–70. https://doi.org/10.1111/j.1471-4159.1984.tb06679.x

Chizhov, A. V., Zefirov, A. V., Amakhin, D. V., Smirnova, E. Y., & Zaitsev, A. V. (2018). Minimal model of interictal and ictal discharges “Epileptor-2.” *PLOS Computational Biology*, 14(5), e1006186. https://doi.org/10.1371/journal.pcbi.1006186
Coulter, D. A., & Steinhäuser, C. (2015). Role of astrocytes in epilepsy. *Cold Spring Harbor Perspectives in Medicine*, 5(3), a022434.
https://doi.org/10.1101/cshperspect.a022434

Cressman, J. R., Ullah, G., Ziburkus, J., Schiff, S. J., & Barreto, E. (2009). The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. *Journal of Computational Neuroscience*, 26(2), 159–170. https://doi.org/10.1007/s10827-008-0132-4

Cunliffe, V. T., Baines, R. A., Giachello, C. N. G., Lin, W. H., Morgan, A., Reuber, M., et al. (2015). Epilepsy research methods update: Understanding the causes of epileptic seizures and identifying new treatments using non-mammalian model organisms. *Seizure*. W.B. Saunders Ltd. https://doi.org/10.1016/j.seizure.2014.09.018

de Almeida, A.-C.G., Rodrigues, A. M., Scorza, F. A., Cavalheiro, E. A., Teixeira, H. Z., Duarte, M. A., et al. (2008). Mechanistic hypotheses for nonsynaptic epileptiform activity induction and its transition from the interictal to ictal state-Computational simulation. *Epilepsia*, 49(11), 1908–1924. https://doi.org/10.1111/j.1528-1167.2008.01686.x

de Curtis, M., Uva, L., Gnatkovsky, V., & Librizzi, L. (2018). Potassium dynamics and seizures: Why is potassium ictogenic? *Epilepsy Research*, 143, 50–59. https://doi.org/10.1016/j.eplepsyres.2018.04.005

Depannemaecker, D., Destexhe, A., Jursa, V., & Bernard, C. (2021, August 1). Modeling seizures: From single neurons to networks. *Seizure*. W.B. Saunders Ltd. https://doi.org/10.1016/j.seizure.2021.06.015

Doyon, N., Vinay, L., Prescott, S. A., & De Koninck, Y. (2016). Chloride Regulation: A Dynamic Equilibrium Crucial for Synaptic Inhibition. *Neuron*, 89(6), 1157–1172. https://doi.org/10.1016/j.neuron.2016.02.030

Durand, D. M., Park, E.-H., & Jensen, A. L. (2010). Potassium dif- fusive coupling in neural networks. *Philosophical transactions of the Royal Society of London. Series B, Biological sciences*, 365(1551), 2347–62. https://doi.org/10.1098/rstb.2010.0050

During, M. J., & Spencer, D. D. (1993). Extracellular hippocampal glutamate and spontaneous seizure in the conscious human brain. *The Lancet*, 341(8861), 1607–1610. https://doi.org/10.1016/0140-6736(93)90754-5

El Houssaini, K., Ivanov, A. I., Bernard, C., & Jursa, V. K. (2015). Seizures, refractory status epilepticus, and depolarization block as endogenous brain activities. *Physical Review E*, 91(1), 010701. https://doi.org/10.1103/PhysRevE.91.010701

Fisher, R. S., Pedley, T. A., Moody, W. J., & Prince, D. A. (1976). The role of extracellular potassium in hippocampal epilepsy. *Archives of neurology*, 33(2), 76–83. https://doi.org/10.1001/archneur.1976.00500020004002

Florence, G., Dahlem, M. A., Almeida, A.-C.G., Bassani, J. W. M., & Kurths, J. (2009). The role of extracellular potassium dynamics in the different stages of ictal bursting and spreading depression: A computational study. *Journal of Theoretical Biology*, 258(2), 219–228. https://doi.org/10.1016/j.jtbi.2009.01.032

Fröhlich, F., Bazhenov, M., Iragui-Madoz, V., & Sejnowski, T. J. (2008). Potassium dynamics in the epileptic cortex: New insights on an old topic. *The Neuroscientist : A Review Journal Bringing Neurobiology, Neurology and Psychiatry*, 14(5), 422–433. https://doi.org/10.1177/1073858408317955

Gloveli, T., Albrecht, D., & Heinemann, U. (1995). Properties of low Mg2+ induced epileptiform activity in rat hippocampal and entorhinal cortex slices during adolescence. *Brain research. Developmental brain research*, 87(2), 145–52. https://doi.org/10.1016/0165-3806(95)00069-P

Golowasch, J., Thomas, G., Taylor, A. L., Patel, A., Pineda, A., Khalil, C., & Nadim, F. (2009). Membrane capacitance measurements revisited: Dependence of capacitance value on measurement

method in nonisopotential neurons. *Journal of Neurophysiology*, 102(4), 2161–2175. https://doi.org/10.1152/jn.00160.2009

Haglund, M. M., & Schwartzkroin, P. A. (1990). Role of Na-K pump potassium regulation and IPSPs in seizures and spreading depres- sion in immature rabbit hippocampal slices. *Journal of Neuro-physiology*, 63(2), 225–239. https://doi.org/10.1152/jn.1990.63.2.225

Hamada, K., Matsuura, H., Sanada, M., Toyoda, F., Omatsu-Kanbe, M., Kashiwagi, A., & Yasuda, H. (2003). Properties of the Na +/K + pump current in small neurons from adult rat dorsal root ganglia. *British Journal of Pharmacology*, 138(8), 1517–1527. https://doi.org/10.1038/sj.bp.0705170

Heinemann, U., Konnerth, A., Pumain, R., & Wadman, W. J. (1986). Extracellular calcium and potassium concentration changes in chronic epileptic brain tissue. *Advances in neurology*, 44, 641–61. PMID: 3518350.

Hille, B. (2001). *Ion channels of excitable membranes*. Sinauer.

Hodgkin, A. L., & Huxley, F. (1952). A quantitative description of membrane current and its application to conduction and excitation in nerve. *Journal of Physiology*, 117(4), 500–544.

Houssaini, K. El, Bernard, C., & Jursa, V. K. (2020). The epileptor model: A systematic mathematical analysis linked to the dynamics of seizures, refractory status epilepticus, and depolarization block. *NeuroE*, 7(2). https://doi.org/10.1523/ENEURO.0485-18.2019

Hubbard, J. A., & Binder, D. K. (2016). *Astrocytes and epilepsy*. Academic Press.

Hübel, N. (2015). *Local dynamics of ion based neuron models for cortical spreading depression, stroke and seizures*. Technische Universität Berlin.

Hübel, N., & Dahlem, M. A. (2014). Dynamics from Seconds to Hours in Hodgkin-Huxley Model with Time-Dependent Ion Concentrations and Buffer Reservoirs. *PLoS Computational Biology*, 10(12), e1003941. https://doi.org/10.1371/journal.pcbi.1003941

Izhikevich, E. (2007). *Dynamical systems in Neuroscience: The geometry of excitability and bursting*. MIT Press.

Jefferys, J. G., & Haas, H. L. (1982). Synchronized bursting of CA1 hippocampal pyramidal cells in the absence of synaptic transmission. *Nature*, 300(5891), 448–50. http://www.ncbi.nlm.nih.gov/pubmed/6292731. Accessed 13 January 2019

Jursa, S., Quilichini, P. P., Ivanov, A. I., & Bernard, C. (2014). On the nature of seizure dynamics. *Brain : A Journal of Neurology*, 137(Pt 8), 2210–2230. https://doi.org/10.1093/brain/awu133

Kager, H., Wadman, W. J., & Somjen, G. G. (2000). Simulated Seizures and Spreading Depression in a Neuron Model Incorporating Inter- stitial Space and Ion Concentrations. *Journal of Neurophysiology*, 84(1), 495–512. https://doi.org/10.1152/jn.2000.84.1.495

Kandel, E. R., Schwartz, J. H., Jessell, T. M., et al. (1981). *Principles of neural science*.

Kofuji, P., & Newman, E. A. (2004). Potassium buffering in the central nervous system. *Neuroscience*, 129(4), 1045–1056. https://doi.org/10.1016/j.neuroscience.2004.06.008

Krinskii, V. I., & Kokoz, Y. M. (1973). Analysis of equations of excit- able membranes-I. Reduction of the Hodgkin-Huxley equations to a second-order system. *Biophysica, 18*, 533–539. Biofizika. 1973 Sep-Oct;18(5):878-85. Russian. PMID: 4751866.

Krishnan, G. P., Filatov, G., Shilnikov, A., & Bazhenov, M. (2015). Electrogenic properties of the Na⁺/K⁺ ATPase control transi- tions between normal and pathological brain states. *Journal of Neurophysiology*, 113(9), 3356–3374. https://doi.org/10.1152/jn.00460.2014

Läuger, P. (1991). *Electrogenic ion pumps*. Sinauer Associates.

Lietsche, J., Imran, I., & Klein, J. (2016). Extracellular levels of ATP and acetylcholine during lithium-pilocarpine induced status epi- lepticus in rats. *Neuroscience Letters*, 611, 69–73. https://doi.org/10.1016/j.neulet.2015.11.028
Lux, H. D., Heinemann, U., & Dietzel, I. (1986). Ionic changes and alterations in the size of the extracellular space during epileptic activity. *Advances in neurology, 44*, 619–39. PMID: 3518349.

McCormick, D. A., & Contreras, D. (2001). On The Cellular and Network Bases of Epileptic Seizures. *Annual Review of Physiology, 63*(1), 815–846. https://doi.org/10.1146/annurev.physiol.63.1.815

McDonald, T., Puchowicz, M., & Borges, K. (2018). Impairments in oxidative glucose metabolism in epilepsy and metabolic treatments thereof. *Frontiers in Cellular Neuroscience, 12*. https://doi.org/10.3389/FNCEL.2018.00274

Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., et al. (2017). SymPy: Symbolic computing in Python. *PeerJ Computer Science, 3*, e103. https://doi.org/10.7717/peerj-cs.103

Miles, R., Blaesse, P., Huberfeld, G., Wittner, L., & Kaila, K. (2012). Chloride homeostasis and GABA signaling in temporal lobe epilepsy. *Jasper's Basic Mechanisms of the Epilepsies*. 4th ed. Bethesda (MD): National Center for Biotechnology Information (US); 2012. PMID: 22787654.

Millman, K. J., & Aivazis, M. (2011). Python for Scientists and Engineers. *Computing in Science & Engineering, 13*(2), 9–12. https://doi.org/10.1109/MCSE.2011.36

Naze, S. (2015). Multiscale Computational Modeling of Epileptic Seizures: from macro to microscopic dynamics. https://www.theses.fr/2015AIXM4023.

Naze, S., Bernard, C., & Jirsa, V. (2015). Computational Modeling of Seizure Dynamics Using Coupled Neuronal Networks: Factors Shaping Epileptiform Activity. *PLOS Computational Biology, 11*(5), e1004209. https://doi.org/10.1371/journal.pcbi.1004209

Olsen, M. L., Khakh, B. S., Skatchkov, S. N., Zhou, M., Lee, C. J., & Rouach, N. (2015). New Insights on Astrocyte Ion Channels: Critical for Homeostasis and Neuron-Glia Signaling. *The Journal of Neuroscience: The Official Journal of the Society for Neuroscience, 35*(41), 13827–13835. https://doi.org/10.1523/JNEUROSCI.2603-15.2015

Øyehaug, L., Østby, I., Lloyd, C. M., Omholt, S. W., & Einevoll, G. T. (2012). Dependence of spontaneous neuronal firing and depolarisation block on astroglial membrane transport mechanisms. *Journal of Computational Neuroscience, 32*(1), 147–165. https://doi.org/10.1007/s10827-011-0345-9

Prinz, A. A., Bucher, D., & Marder, E. (2004). Similar network activity from disparate circuit parameters. *Nature Neuroscience, 7*(12), 1345–1352. https://doi.org/10.1038/nn1352

Quilichini, P. P., Diabira, D., Chiron, C., Ben-Ari, Y., & Gozlan, H. (2002). Persistent epileptiform activity induced by low Mg2+ in intact immature brain structures. *The European journal of neuroscience, 16*(5), 850–60. https://doi.org/10.1046/j.1460-9568.2002.02143.x

Raimondo, J. V., Burman, R. J., Katz, A. A., & Akerman, C. J. (2015). Ion dynamics during seizures. *Frontiers in Cellular Neuroscience, 9*, 419. https://doi.org/10.3389/FNCEL.2015.00419

Rangroo Thrane, V., Thrane, A. S., Wang, F., Cotrina, M. L., Smith, N. A., Chen, M., et al. (2013). Ammonia triggers neuronal disinhibition and seizures by impairing astrocyte potassium buffering. *Nature Medicine, 19*(12), 1643–1648. https://doi.org/10.1038/nm.3400

Rodrigues, A. M., Santos, L. E. C., Covolan, L., Hamani, C., & Almeida, A.-C.G. (2015). pH during non-synaptic epileptiform activity—computational simulations. *Physical Biology, 12*(5), 056007. https://doi.org/10.1088/1478-3975/12/5/056007

Saggio, M. L., Crisp, D., Scott, J. M., Karoly, P., Kuhlmann, L., Nakatani, M., et al. (2020). A taxonomy of seizure dynamotypes. *eLife, 9*. https://doi.org/10.7554/eLife.55632

Saggio, M. L., Spiegler, A., Bernard, C., & Jirsa, V. K. (2017). Fast-Slow Bursters in the Unfolding of a High Codimension

Singularity and the Ultra-slow Transitions of Classes. *Journal of Mathematical Neuroscience, 7*(1), 7. https://doi.org/10.1186/s13408-017-0050-8

Santhakumar, V., Aradi, I., & Soltesz, I. (2005). Role of mossy fiber sprouting and mossy cell loss in hyperexcitability: A network model of the dentate gyrus incorporating cell types and axonal topography. *Journal of Neurophysiology, 93*(1), 437–453. https://doi.org/10.1152/jn.00777.2004

Scholl, U. I., Choi, M., Liu, T., Ramaekers, V. T., Häusler, M. G., Grimmer, J., et al. (2009). Seizures, sensorineural deafness, ataxia, mental retardation, and electrolyte imbalance (SeSAME syndrome) caused by mutations in KCNJ10. *Proceedings of the National Academy of Sciences of the United States of America, 106*(14), 5842–5847. https://doi.org/10.1073/pnas.0901749106

Soltesz, I., & Staley, K. (2008). *Computational neuroscience in epilepsy*. Elsevier.

Somjen, G. G. (2001). Mechanisms of spreading depression and hypoxic spreading depression-like depolarization. *Physiological Reviews, 81*(3), 1065–1096. https://doi.org/10.1152/physrev.2001.81.3.1065

Staley, K. (2015). Molecular mechanisms of epilepsy. *Nature Neuroscience, 18*(3), 367–372. https://doi.org/10.1038/nn.3947

Stefanescu, R. A., Shivakeshavan, R. G., & Talathi, S. S. (2012, December 1). Computational models of epilepsy. *Seizure*. W.B. Saunders. https://doi.org/10.1016/j.seizure.2012.08.012

Strauss, U., Zhou, F.-W., Henning, J., Battefeld, A., Wree, A., Köhling, R., et al. (2008). Increasing Extracellular Potassium Results in Subthalamic Neuron Activity Resembling That Seen in a 6-Hydroxydopamine Lesion, 99(6), 2902–2915. https://doi.org/10.1152/JN.00402.2007

Tejada, J., Garcia-Cairasco, N., & Roque, A. C. (2014). Combined Role of Seizure-Induced Dendritic Morphology Alterations and Spine Loss in Newborn Granule Cells with Mossy Fiber Sprouting on the Hyperexcitability of a Computer Model of the Dentate Gyrus. *PLoS Computational Biology, 10*(5), 1003601. https://doi.org/10.1371/journal.pcbi.1003601

Traynelis, S. F., & Dingledine, R. (1988). Potassium-induced spontaneous electrographic seizures in the rat hippocampal slice. *Journal of Neurophysiology, 59*(1), 259–276. https://doi.org/10.1152/jn.1988.59.1.259

Traynelis, S. F., & Dingledine, R. (1989). Modification of potassium-induced interictal bursts and electrographic seizures by divalent cations. *Neuroscience Letters, 98*(2), 194–199. https://doi.org/10.1016/0304-3940(89)90509-0

Ullah, G., & Schiff, S. J. (2010). Assimilating seizure dynamics. *PLoS Computational Biology, 6*(5), e1000776. https://doi.org/10.1371/journal.pcbi.1000776

Uva, L., Trombin, F., Carriero, G., Avoli, M., & de Curtis, M. (2013). Seizure-like discharges induced by 4-aminopyridine in the olfactory system of the *in vitro* isolated guinea pig brain. *Epilepsia, 54*(4), 605–615. https://doi.org/10.1111/epi.12133

Walz, W. (2000). Role of astrocytes in the clearance of excess extracellular potassium. *Neurochemistry international, 36*(4–5), 291–300. https://doi.org/10.1016/S0197-0186(99)00137-0

Wang, L., Dufour, S., Valiante, T. A., & Carlen, P. L. (2016). Extracellular Potassium and Seizures: Excitation, Inhibition and the Role of Ih. *International Journal of Neural Systems, 26*(08), 1650044. https://doi.org/10.1142/S0129065716500441

Wang, Y., Trevelyan, A. J., Valentin, A., Alarcon, G., Taylor, P. N., & Kaiser, M. (2017). Mechanisms underlying different onset patterns of focal seizures. *PLOS Computational Biology, 13*(5), e1005475. https://doi.org/10.1371/journal.pcbi.1005475

Wei, Y., Ullah, G., & Schiff, S. J. (2014a). Unification of neuronal spikes, seizures, and spreading depression. *Journal of Neuroscience, 34*(35), 11733–11743. https://doi.org/10.1523/jneurosci.0516-14.2014
Wei, Y., Ullah, G., Ingram, J., & Schiff, S. J. (2014b). Oxygen and seizure dynamics: II. Computational modeling. *Journal of Neurophysiology*, **112**(2), 213–223. https://doi.org/10.1152/jn.00541.2013

Wendling, F., Benquet, P., Bartolomei, F., & Jirsa, V. (2016). Computational models of epileptiform activity. *Journal of Neuroscience Methods*, **260**, 233–251. https://doi.org/10.1016/J.JNEUMETH.2015.03.027

Zilberter, Y., & Zilberter, M. (2017). The vicious circle of hypometabolism in neurodegenerative diseases: Ways and mechanisms of metabolic correction. *Journal of Neuroscience Research*, **95**(11), 2217–2235. https://doi.org/10.1002/JNR.24064

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Authors and Affiliations

Damien Depannemaecker¹,²,³ · Anton Ivanov¹ · Davide Lillo¹ · Len Spek¹,⁴ · Christophe Bernard¹ · Viktor Jirsa¹

¹ Aix Marseille Univ, Inserm, INS, Institut de Neurosciences des Systèmes, Aix Marseille University, Marseille, France

² Disciplina de Neurociência, Departamento de Neurologia e Neurocirurgia, Universidade de Federal de São Paulo, São Paulo, Brazil

³ Department of Integrative and Computational Neuroscience, Institut Neuro-PSI, CNRS, Gif-sur-Yvette, France

⁴ Department of Applied Mathematics, University of Twente, Enschede, The Netherlands