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

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. Supplementary Information The online version contains supplementary material available at 10.1007/s10827-022-00811-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 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(Saggio et al., , 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 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 . 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.

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/ ddepa nn/ TheMo del.
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 ). 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 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 voltagegated 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) 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 [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 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.

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

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)  . This is consistent with the observations described in (Fisher et al., 1976 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 nonpathological 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.

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 I pump 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 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.

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

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 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 Fig. 7 Fixed points of the fast subsystem. Fixed point of the fast subsystem (V m ) 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 subsys-tem, 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 ; 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+ , 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 . Since it is possible to trigger similar SLEs via totally different biophysical mechanisms , 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 . 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 . 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 saddlenode/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 , 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.

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).
With currents: And conductance variables: 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): 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 n inf (V) and a constant parameter τ n . In our model, n inf (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 n inf (V) of our model and 1/τ(V), and n 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 (12) dn dt = n (1 − n) − n n  . The value of τ 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 -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.
Thanks to these reductions, concentration variations are calculated as follow: All simulations were obtained thanks to numerical methods using odeint function from SciPy library (Millman & Aivazis, 2011

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:// creat iveco mmons. org/ licen ses/ by/4. 0/.