A unified physiological framework of transitions between seizures, status epilepticus 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 status epilepticus (SE) 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 SE and DB, we show that SE-like events and DB can also occur at the single cell level. Thus, seizures, SE 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 status epilepticus.


Introduction
Seizures are part of the repertoire of built-in activities of neuronal networks as they can be triggered in most brain regions from most species [1]. Several conceptual frameworks have been proposed to explain seizure dynamics [2][3][4][5][6][7][8]. The predominant framework assumes that the majority of seizure onsets and offsets correspond to bifurcations in the electrophysical variables [1,9], although there exist other non-bifurcation types [10]. This framework has been generalized by Saggio and col. [11,12]. 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 [1,12]. A qualitative analysis of the Epileptor reveals that seizures, SE and DB co-exist, and that multiple types of transitions from one type of activity to the other are possible, as verified experimentally [11,13,14]. Since it is phenomenological, the Epileptor model does not provide direct insight regarding underlying biophysical mechanisms. The phenomenological model imposes strong constrains in terms of dynamics. Numerous neuronal network models, including biophysically realistic ones, have been developed to study seizures, SE or DB mechanisms [2,4,[15][16][17][18][19][20][21]. These models contain too many parameters to perform a detailed bifurcation analysis, thus preventing bridging the gap between phenomenological and biophysical approaches. However, with the guidance of phenomenological modeling, design of neuronal spiking network including several biophysical features has been performed [3]. In this work, transitions between states of the neuronal network are ensured by slow variable representing extracellular environmental fluctuation.
Although seizures and DB are generally observed at the neuronal network scale, their dynamical equivalence can be found at the single cell level [22][23][24][25][26][27][28][29], which can be used to study human epilepsy [30][31][32]. Moreover, previous works [3,33,34] show that dynamical features are preserved when going from the network to the single cell level. It therefore seems appropriate to consider a biophysical model at the single cell scale exhibiting dynamic properties identified in the generic model and in which the transitions between the different states are provided by a slow variable describing the variations within the extracellular milieu.
Bursting activity in neurons can be described in terms of bifurcations [35,36], and different single cell biophysical models have been proposed, which can model SLE and DB, but not SElike events [22,25,26,29,[37][38][39][40][41], although, to the best of our knowledge, these activities have not been observed experimentally in isolated neurons. 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 [22,25,26,29,[37][38][39][40], oxygen level [38,42], volume [37,40] and interaction with glial cells [29,40]. 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, including the extracellular slow variations, could be reduced to four variables, while presenting a bursting pattern corresponding to the most common ones encountered at the network level, i.e. the SN/homoclinic pair [1]. Since SLEs, SE and DB pertain to the dynamic repertoire of The Epileptor and of biological neuronal networks [11,13], biophysical models should be able to reproduce all three types of activities. The goal of the present study is to identify candidate mechanisms from physiology robustly leading to the time scale separation and trajectories in SLEs. We will take guidance from the careful dynamic analyses performed in previous works.
In generic models, the dynamics is well understood [1,11] but these models rarely offer direct biophysical insight as they use abstract parameters. Important works have been done to understand the link between phenomenological and biophysical models [3,43]. In order to explore the dynamics repertoire, the high dimensionality of detailed biophysical models must be reduced. A minimal model of interictal and SLEs has been introduced as Epileptor-2, using increase in [K]o to trigger burst discharges and restoration of the K + gradient via the sodiumpotassium pump to stop SLEs [22]. The work of Saggio et al. [11] is a generalization of the dynamics found in Epileptor and Epileptor-2. However, Epileptor-2 does not produce the SN/homoclinic bifurcation consistently found experimentally, and the model does not generate SE or DB. Moreover, when the fast discharges stop, the slow variables of Epileptor-2 continue oscillating, whereas that of Epileptor 1 do not. This can be understood on the basis of the results of Saggio et al. [11]. In the present reduction, we used a Hodgkin-Huxley-like single cell model, and we imposed several constraints: SLE, SE and DB [13], as well as the SN/homoclinic bifurcation must be present. We have thus generated a minimal slow/fast system preserving biophysical representation and satisfying these constraints. A variable acting on a slow time scale is necessary to drive the system through different activities (e.g. from SLE to DB). In a cell, numerous processes can occur on a slow time scale and act as a slow variable, including changes in ion concentration, metabolism, phosphorylation levels, and transcription. Ion homeostasis regulation is critical to maintaining neuron function, and many sub-cellular mechanisms are involved in this process [44]. Fluctuations of ion concentrations in the extracellular space modulate the electrophysiological activity of a single neuron [26,42].
Computational simulations show that potassium could be responsible for local synchronization [55] and is an important parameter in neural dynamics [26,39,41,42]. In addition, in experimental models, the transition to DB correlates with a much larger increase of [K]o as compared to SLEs [13,56]. 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, SE and DB, reproducing patterns found in single neurons recorded experimentally during seizures.

Results
Our goal was to construct a biophysical single neuron model that can reproduce the different firing patterns recorded when extracellular potassium is experimentally increased, while keeping it sufficiently simple to allow a bifurcation analysis. The model is schematized in   [61], the corresponding intracellular activity is not known, however the sustained firing pattern in the model cell resembles the regular field activity recorded during SE-like events in vitro [61]. The sustained DB at the single cell level corresponds to what is observed experimentally during network spreading depolarization when [K]o reaches high levels [62].
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 detailed these dynamical interactions for the different patterns of activity, following the order of appearance when [K]bath increase. Nernst potential of potassium with specific Y axis on the right side. 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.

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 [11], and confirmed by numerical methods (using SymPy [63] and SciPy

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)  (eq.20). This is consistent with the observations described in [45].   Fig. 3) as also reported experimentally [13]. 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. 2e and 2f We conclude that the model behaves as expected from the biological observations, when experimentally increasing of [K]bath. These simulations were obtained when using a "healthy" situation, i.e. as if recording a "control" neuron. In the next section, we model a "pathological" situation for which where the regulatory mechanisms of neuronal homeostasis are affected.

Influence of other parameters
We then aimed to identify relevant parameters that could describe "healthy" and "pathological" states. Experimental data show that impairment in potassium buffering by glial cells leads to pathological behavior [65][66][67][68]. Three model parameters correspond to homeostasis regulation, involving two mechanisms: the ion exchange capacity between compartments (ε and γ parameters in the model), and the maximum capacity of the Na/K pump (ρ in the model). A variation of ε corresponds to a degradation of the interaction with glial cells [67,[69][70][71] which normally ensures the regulation of the extracellular concentration of K + .
Homeostasis of intracellular ions is also critical, and a variation of γ corresponds to an impairment of such mechanisms, not detailed in this model, such as co-transporters and exchangers [72,73]. Changes in both parameters can also be considered. In the model, they thanks to their specific dynamics and resulting shapes [11,74].
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 SE, the ATP concentration augments due to high needs and then decreases [24]. 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). All these observations show that the model presents a behavior consistent with experimental observations. Importantly, the biophysical model is able to reproduce general patterns of activities (i.e. periodic events) as generated by the phenomenological model [13].
Phenomenological models, which present a minimal number of variables and parameters, allow an exhaustive study of the dynamics. The biophysical model used here contains too many parameters for an exhaustive study of the dynamics, but reducing the number of variables will allow a comparison with the generic model. In the next section, we analyze the dynamics of the model.

Dynamical observations
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. This type of phenomenon, corresponding to slow fast systems, has been extensively studied from a theoretical point of view, in particular for neural activity [11,35]. In this subsection, to allow a better correspondence with the theoretical framework, 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 [11,74]: 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 [11], 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. In our model, two equations describe the slow subsystem (Eq. 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.  considering the variables from the slow subsystem as parameters. We used a numerical methods with SymPy [63] and SciPy [64] 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.

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 [1,23,60,61,75] and described by generic models [1,11] . 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 (SLE and SE-like event), via the variation of only one parameter. All parameter values were within biophysical ranges (Table1) [72,73]. The model has two main limitations. The fast system describes only intrinsic excitability and does not include synaptic currents. And, 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 [1]; 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+ [76], Cl - [77,78] and neurotransmitters during seizures [79,80]. Likewise, decreasing extracellular Ca 2+ leads to seizures [50], which are characterized by SN/homoclinic bifurcations [1]. Since it is possible to trigger similar SLEs via totally different biophysical mechanisms [1], we propose that the K +dependent mechanism we describe, is one among many the possible paths leading to the 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 [3]. 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 SE and DB, which corresponding pattern have been found in dynamics of our model. However, 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 [21,26,[38][39][40]42,43,81]. Descriptions of ion concentration dynamics for bursting have been done by Barreto et al. [39], 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 and al. [37]. As different models can lead to similar dynamics [82], this may suggest that different minimalist models are possible to obtain a unified framework. In our work, we proposed a conductancebased model of the neuronal membrane, exhibiting an extended repertoire of behavior and introducing status epilepticus-like events in a unified framework. Another difference with previous work 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 [1]. Our model does not take into account variation of volume or oxygen homeostasis as in [37] but, only variations of ion concentrations, driven by diffusion of potassium from EB. 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 Øyehauget al. [40] 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 [21,[38][39][40], 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. [22] proposed a biophysical model (Epileptor-2) of ictal activities based on the Epileptor [1], 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 saddlenode/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 (see Fig 10. in [22]), 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. Although the Epileptor-2 is not an "intrinsic" Slow/Fast dynamical system, indeed, this model does not describe an independent node as it takes in account the influence from synaptic inputs from neuronal population. In our model, the observed dynamics, is only due to internal interactions between three compartments.
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 [1,11,13], 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, status epilepticus and depolarization block [1,83]. 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 [72,84,85] and the membrane capacitance in [86]. 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).
(1)  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 [87] the activation variable of K+ channels is determined by the equation (eq.12): (12) dn dt = α n (1 − n) − β n n 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 [87,88], 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 has been build thanks to experiments done on the squid giant axon, which present differences from on recording of mammalians 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 [87]. With the same logic, we can consider the percentage of all population of channels opened. But because this is a very fast mechanism [72], it can be considered as an instantaneous function of V [74] (eq.10). Krinskii and Kokoz[89] 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 [29,90], electroneutrality permits the Eq. (12), and so to the Eq. (13). The ratio (Cm γ)/ωi is very small (<10 -5 ) and so, the right-hand side of Eq.(13) 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 [91]. So, in our model, the chloride concentration remains constant.