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., 2016). The predominant framework posits that the majority 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 mathematical 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 phenomenological, the Epileptor model does not provide direct insight regarding underlying biophysical mechanisms.

Detailed biophysical network models have been developed to investigate seizure mechanisms (Rodrigues et al., 2015; Santhakumar et al., 2005; Tejada et al., 2014). These models contain too many parameters and variables to perform a detailed analysis of their dynamics repertoire. Although seizures have always been identified 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 bifurcations (E. Izhikevich, 2007; E. M. Izhikevich, 2000), different 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 membrane electrophysiological activities. The slow subsystem can be represented by variations of different slow variables 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]o) 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 activities (de Almeida et al., 2008; Jefferys & Haas, 1982). In addition, in experimental models, the transition to DB correlates with a much larger increase of [K]o as compared to SLEs (El Houssaini et al., 2015; Gloveli et al., 1995). Theoretical 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]o 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 during seizures.

2 Results

Our goal was to construct a biophysical model of a single 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 surrounding environment with three compartments (external bath, extracellular space, intracellular space). The corresponding code is available on github at https://github.com/ddepann/TheModel.

Fig. 1
figure 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 Vm, 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)

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 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 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]bath increases.

Fig. 2
figure 2

Qualitative mode of behavior of the membrane potential and Nernst potentials. In blue: time series of the membrane potential Vm 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 [K]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

Fig. 3
figure 3

Variation of extracellular potassium concentration as a function of [K]bath. Minimal and maximal external potassium [K]o and mean (dash line) concentration observed during simulations done for different values of the parameter [K]bath. Due to diffusion from the external bath, increasing [K]bath leads to variations in [K]o. Different patterns are observed for each range of [K]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]bath where [K]o periodically oscillates

2.1 Resting states, spike train and tonic spiking in low [K]bath

Resting state is found when [K]bath is around the normal value of [K]o (called [K]0,o, see Methods section). If [K]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]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]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]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 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 ∆[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 IK. The influence of IK in the equation of ∆[K]I (Eq. 3), explains the decrease of ∆[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, it is not the cause of the event but a consequence of homeostasis dysregulation (i.e. augmentation of [K]bath).

Table 1 Parameters values

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 ε. Homeostasis of intracellular ions is also important to consider. The parameter γ 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 γ. Varying the time constants of the slow subsystem (ε and γ), leads to different bi-stable behaviors. Two examples are shown in Fig. 4(b) with γ = 0.04, ε = 0.002, (d) γ = 0.06, ε = 0.002, and (f) γ = 0.08, ε = 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 ρ. 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 ρ. This parameter also influences the shape of Ipump response as a function of [Na]i and [K]o (Fig. 5a). For large values of ρ, the pump is activated for lower value of [Na]i and [K]o (Fig. 5a). We find that burst duration changes with ρ for a fixed [K]bath (Fig. 5b), where a faster activation (higher ρ) leads to shorter bursts. The augmentation of ρ does not necessary lead to an increase of Ipump; it affects the general dynamics of the whole system (Fig. 5c).

Fig. 4
figure 4

Time series of membrane potential, ∆[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 γ = 0.04, ε = 0.002, (c) Burst with [K]bath = 12.5 mM, (d) Burst with [K]bath = 12.5 mM, and γ = 0.06, ε = 0.002, (e) SLE with [K]bath = 16 mM, (f) SLE with [K]bath = 16 mM and γ = 0.08, ε = 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 ∆[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 γ and ε influence the shape of Vm allowing the system to exhibit SN or SNIC bifurcation at the onset of the events

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 IK (Eq. (6)) 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.

Fig. 5
figure 5

Influence of the activity of the Na/K-pump. (a) Ipump function for ρ = 25 (green), ρ = 1000 (orange), ρ = 2500 (blue). The initial slope when the system moves away from the concentrations at rest is affected, explaining the modification of the influence of Ipump in the dynamic of the system. (b) Burst duration as a function of ρ for [K]bath = 14.0 mM. Bursts have shorter durations for higher value of ρ. (c) Minimal and maximal pump current, Ipump, observed during simulation done with [K]bath = 14.0 mM. The range of Ipump decreases for higher ρ values

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 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]i and [K]g 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]i, 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.

Table 2 Physiological reference values
Fig. 6
figure 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

Fig. 7
figure 7

Fixed points of the fast subsystem. Fixed point of the fast subsystem (Vm) 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, 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

Fig. 8
figure 8

modification in gating variables. (a) ninf of our model in blue, and ninf and 1/τ of the Hodgkin-Huxley model respectively in dash blue and red, function of the membrane potential. (b) Response of the fast subsystem of our model to step current stimulation (red) with three different values of τ (0.1, 0.25, 0.5 ms). The value of τ influence the frequency rate spike for a same injected current

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 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 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 Ca2+ (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 Ca2+ 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 IK 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]o 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.

Table 3 Initial values

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}}}\left({I}_{Cl}+ {I}_{Na}+ {I}_{K}+ {I}_{pump}\right)$$
(1)
$$\frac{dn}{dt}=\frac{n_\infty\left(V\right)-n}{\tau_n}$$
(2)
$$\frac{d\Delta{\left[K\right]}_i}{dt}=-\frac\gamma{\omega_i}(I_K-2I_{pump})$$
(3)
$$\frac{d[{K]}_{g}}{dt}=\varepsilon ({\left[{K}\right]}_{{bath}}-{\left[{K}\right]}_{{o}})$$
(4)

With currents:

$${I}_{Na}={(g}_{Na,l}+{g}_{Na}{m}_{\infty }\left(V\right)h\left(n\right))(V-26.64{log}(\frac{{\left[Na\right]}_{o}}{{\left[Na\right]}_{i}}))$$
(5)
$${I}_{K}={(g}_{K,l}+{g}_{K}n)(V-26.64{log}(\frac{{\left[K\right]}_{o}}{{\left[K\right]}_{i}}))$$
(6)
$${I}_{Cl}={g}_{Cl}(V+26.64{log}(\frac{{\left[Cl\right]}_{o}}{{\left[Cl\right]}_{i}}))$$
(7)
$${I}_{pump}=\rho \frac{1}{1+{exp}(\frac{1}{2}(21-{\left[Na\right]}_{i}))}\frac{1}{1+{exp}(5.5-{\left[K\right]}_{o})}$$
(8)

And conductance variables:

$${n}_{\infty }\left({V}\right)=\frac{1}{1+{exp}(\frac{1}{18}(-19-{V}))}$$
(9)
$${m}_{\infty }\left({V}\right)=\frac{1}{1+{exp}(\frac{1}{12}(-24-{V}))}$$
(10)
$$h\left({n}\right)=1.1-\frac{1}{1+{exp}(-8({n}-0.4))}$$
(11)

The fast subsystem of the model, (Eq. (1)&(2)), is a reduction and simplification of conductance-based models, first describe 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}}\left(1-{n}\right)-{\beta }_{{n}}{n}$$
(12)

where β(V) and α(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 ninf(V) and a constant parameter τn. In our model, ninf(V) is the probability to find a channel at open state at a given membrane potential while τ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 ninf(V) of our model and 1/τ(V), and ninf(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 (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 (Cm γ)/ωi is very small (< 10–5) 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}(\frac{d\Delta{\left[K\right]}_i}{dt}+\frac{d\Delta{\left[Na\right]}_i}{dt}+\frac{d\Delta{\left[Cl\right]}_i}{dt})$$
(13)
$$\Delta{\left[K\right]}_i+\Delta{\left[N\alpha\right]}_i+\Delta{\left[Cl\right]}_i=\frac{C_{m}\gamma}{W_i}\left({V}-{V_0}\right)$$
(14)

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

$$\Delta{\left[N\alpha\right]}_i=-\Delta{\left[K\right]}_i$$
(15)
$$\Delta{\left[N\alpha\right]}_o=-\beta\Delta{\left[N\alpha\right]}_i$$
(16)
$$\Delta{\left[K\right]}_o=-\beta\Delta{\left[K\right]}_i$$
(17)
$${\left[{K}\right]}_{{i}}= {\left[{K}\right]}_{0,{i}} + \Delta {\left[{K}\right]}_{{i}}$$
(18)
$${\left[{N a}\right]}_{{i}}= {\left[{N a}\right]}_{0,{i}} + \Delta {\left[{N a}\right]}_{{i}}$$
(19)
$${\left[{N a}\right]}_{{o}}= {\left[{N a}\right]}_{0,{o}} + \Delta {\left[{N a}\right]}_{{o}}$$
(20)
$${\left[{K}\right]}_{{o}} = {\left[{K}\right]}_{0,{o}} + \Delta {\left[{K}\right]}_{{o}} + {[{K}]}_{{g}}$$
(21)

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