Skip to main content
Advertisement
  • Loading metrics

Minimal model of interictal and ictal discharges “Epileptor-2”

  • Anton V. Chizhov ,

    Roles Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing

    anton.chizhov@mail.ioffe.ru

    Affiliations Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia, Computational Physics Laboratory, Ioffe Institute, Saint Petersburg, Russia

  • Artyom V. Zefirov,

    Roles Formal analysis, Investigation, Software

    Affiliations Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia, Computational Physics Laboratory, Ioffe Institute, Saint Petersburg, Russia

  • Dmitry V. Amakhin,

    Roles Data curation, Investigation, Methodology, Software

    Affiliation Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia

  • Elena Yu. Smirnova,

    Roles Data curation

    Affiliations Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia, Computational Physics Laboratory, Ioffe Institute, Saint Petersburg, Russia

  • Aleksey V. Zaitsev

    Roles Conceptualization, Supervision, Writing – review & editing

    Affiliations Laboratory of Molecular Mechanisms of Neural Interactions, Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences, Saint Petersburg, Russia, Institute of Experimental Medicine, Almazov National Medical Research Centre, Saint Petersburg, Russia

Correction

12 Sep 2019: Chizhov AV, Zefirov AV, Amakhin DV, Smirnova EY, Zaitsev AV (2019) Correction: Minimal model of interictal and ictal discharges “Epileptor-2”. PLOS Computational Biology 15(9): e1007359. https://doi.org/10.1371/journal.pcbi.1007359 View correction

Abstract

Seizures occur in a recurrent manner with intermittent states of interictal and ictal discharges (IIDs and IDs). The transitions to and from IDs are determined by a set of processes, including synaptic interaction and ionic dynamics. Although mathematical models of separate types of epileptic discharges have been developed, modeling the transitions between states remains a challenge. A simple generic mathematical model of seizure dynamics (Epileptor) has recently been proposed by Jirsa et al. (2014); however, it is formulated in terms of abstract variables. In this paper, a minimal population-type model of IIDs and IDs is proposed that is as simple to use as the Epileptor, but the suggested model attributes physical meaning to the variables. The model is expressed in ordinary differential equations for extracellular potassium and intracellular sodium concentrations, membrane potential, and short-term synaptic depression variables. A quadratic integrate-and-fire model driven by the population input current is used to reproduce spike trains in a representative neuron. In simulations, potassium accumulation governs the transition from the silent state to the state of an ID. Each ID is composed of clustered IID-like events. The sodium accumulates during discharge and activates the sodium-potassium pump, which terminates the ID by restoring the potassium gradient and thus polarizing the neuronal membranes. The whole-cell and cell-attached recordings of a 4-AP-based in vitro model of epilepsy confirmed the primary model assumptions and predictions. The mathematical analysis revealed that the IID-like events are large-amplitude stochastic oscillations, which in the case of ID generation are controlled by slow oscillations of ionic concentrations. The IDs originate in the conditions of elevated potassium concentrations in a bath solution via a saddle-node-on-invariant-circle-like bifurcation for a non-smooth dynamical system. By providing a minimal biophysical description of ionic dynamics and network interactions, the model may serve as a hierarchical base from a simple to more complex modeling of seizures.

Author summary

In pathological conditions of epilepsy, the functioning of the neural network crucially depends on the ionic concentrations inside and outside neurons. A number of factors that affect neuronal activity is large. That is why the development of a minimal model that reproduces typical seizures could structure further experimental and analytical studies of the pathological mechanisms. Here, on a base of known biophysical models, we present a simple population-type model that includes only four principal variables, the extracellular potassium concentration, the intracellular sodium concentration, the membrane potential and the synaptic resource diminishing due to short-term synaptic depression. A simple modeled neuron is used as an observer of the population activity. We validate the model assumptions with in vitro experiments. Our model reproduces ictal and interictal events, where the latter result in bursts of spikes in single neurons, and the former represent the cluster of spike bursts. Mathematical analysis reveals that the bursts are spontaneous large-amplitude oscillations, which may cluster after a saddle-node on invariant circle bifurcation in the pro-epileptic conditions. Our consideration has significant bearing in understanding pathological neuronal network dynamics.

Introduction

A simple canonical mathematical model of epileptic discharges has been proposed by V. Jirsa et al. [1] based on an analysis of experimental recordings of epileptic discharges in humans and animal models. A comparison of the characteristic features of the canonical bifurcations of simple dynamical systems with the elements of the experimental signals revealed the most typical type of transitions to and from the ictal state. A model called Epileptor, which reproduces this behavior, has been proposed. The model can be expressed in a set of six ordinary differential equations [2]. This model led to a series of predictions and applications in physiology, neuroinformatics, and translational medicine [24]. A disadvantage of the Epileptor is a lack of physical meaning for its governing variables. A biophysical interpretation would help clarify the experimental validation of the model’s predictions. In the present work, an alternative model with a similar mathematical simplicity but that operates with physically meaningful variables is proposed. The model is referred to as the Epileptor-2. The set of variables chosen to describe the epileptic state transitions along with the interpretation of the terms included in the equations constitute the primary model predictions. Namely, the model distinguishes between main variables: two ionic concentrations of the extracellular potassium and intracellular sodium, a membrane potential, and a synaptic resource. The selection of these variables as well as the forms of the governing equations with experimental observations taken from both the experiments and the literature are justified and described in detail.

While formulating the model assumptions, special attention was focused on two specific experimental observations. First, in a number of experiments, IDs were formed from a clustered number of short discharges [58]. These elementary discharges are similar to IIDs. During each IID-like event, a principal neuron typically generated a burst of spikes. Thus, the ID represents a burst of bursts of spikes. Previous biophysical seizure-modeling attempts that considered the ionic dynamics in a single neuron model [9,10] or a network [1113] did not show this feature. The second observation suggested a crucial role of the Na-K pump in the termination of each ID. In fact, the extracellular potassium concentration increased during the initial phase of the ID and began to decrease before the end of the ID [4,8,14], whereas the intracellular sodium concentration increased until the actual end of the ID [15]. Some simulations showed a similar behavior of ionic concentrations [9,10]. Therefore, the ionic dynamics may be explained by the activation of the Na-K pump, and the consequent decrease in the extracellular potassium concentration may lead to ID termination [12]. The two experimental observations described have not been modeled together. The Results section of the paper is organized as follows. First, the proposed model is presented, which is then justified. The simulations of the interictal and ictal discharges are described and compared with experimental data in order to validate the main assumptions of the model. Then, the main component of the Epileptor-2’s dynamics is discussed, which is an interictal-like discharge observed as a short burst of spikes. To study these discharges, the case of constant ionic conductance was considered, the Epileptor-2 model was reduced to a second-order system of stochastic differential equations, and its mathematical analysis was carried out. It is shown that the deterministic system of the reduced model has the only stable equilibrium as an attractor; however, under noise, the model exhibits a rather complex behavior of the stochastic generation of bursting, similar to the behavior of the stochastic Hindmarsh-Rose model [16]. For weak noise, random states concentrate near the equilibrium. With an increase of noise intensity, random trajectories can travel far from the stable equilibrium, and along with small-amplitude oscillations around the equilibrium, bursts are observed. These bursts constitute the interictal and ictal events in the simulations of Epileptor-2. The slow oscillations of the potassium and sodium ionic concentrations were then analyzed, which determine the origination and shape of seizures.

Results

Governing equations of the Epileptor-2

The proposed model consists of three subsystems that describe: (i) the ionic dynamics, (ii) the neuronal excitability, and (iii) a neuron-observer (Fig 1A).

thumbnail
Fig 1. Schematic of Epileptor-2 and the described mechanisms of interictal and ictal discharges.

(A) The model consists of the equations that describes the ionic dynamics, the neuronal population excitation and the activity of a single representative neuron driven by the population. (B) Mechanisms of the discharges, explained in the main text.

https://doi.org/10.1371/journal.pcbi.1006186.g001

The proposed population model consists of the following equations: (1) (2) (3) (4)

All the terms of these equations are explained in the next section. Here, only the main variables are introduced, which are as follows: [K]o and [Na]i represent extracellular potassium and intraneuronal sodium concentrations, respectively; V(t) is the membrane depolarization; xD(t) is the synaptic resource; ν(t) is the firing rate of an excitatory population; and an inhibitory population firing rate is assumed to be proportional to ν(t). The dynamics described by these equations is driven ν(t), which is calculated with a sigmoidal input-output function: (5) where [x]+ is equal to x for the positive argument and 0 otherwise. The input current u(t) includes the potassium depolarizing current, the synaptic drive, and the noise ξ(t), respectively: (6)

The potassium reversal potential is obtained from the ion concentrations via the Nernst equation: (7)

The Na+/K+ pump current is taken from Cressman et al. [9] in the form: (8)

The basic set of parameters is as follows. The time constants for ionic dynamics, membrane polarization, and synaptic depression are τK = 100s, τNa = 20s, τm = C/gL = 10ms, and τD = 2s; the concentration increments at spike are δ[K]o = 0.02mM, δ[Na]i = 0.03mM, δxD = 0.01; the Gaussian white noise ξ(t) has zero mean and unity dispersion, ⟨ξ(t)ξ(t')⟩ = τm δ(tt'); the noise amplitude is σ/gL = 25mV; the maximum pump flux is ρ = 0.2mM/s; the volume ration is γ = 10 (close to 7 given in [17]); the postsynaptic charge is Gsyn/gL = 5mVs; the potassium leak conductance is gK,leak/gL = 0.5; the initial extracellular and elevated bath potassium concentrations are and [K]bath = 8.5mM; ; initial and resting intracellular sodium concentrations are the same, . The parameters of the input-output function are the maximal rate νmax = 100Hz, the threshold potential Vth = 25mV, and the gain kν = 20mV.

A representative neuron was modeled with a quadratic integrate-and-fire neuron [18, 19]. The equations for the membrane potential U(t) are as follows: (9) (10) with the parameters: gU = 0.4nS/mV, CU = 200pF, VT = 25mV, Vreset = −50mV, U1 = −60mV, U2 = −40mV; an initial condition is U = −70mV.

Phenomenological derivation and justification of the model

Eq 1. The model states that the elevation of the extracellular potassium concentration [K]o plays a primary role in self-regenerating IDs. Before each ID, the extracellular potassium concentration accumulates after a series of IIDs. Each IID involves the activation of interneurons [5,7], which results in opening GABAergic channels that are permeable for chloride and bicarbonate [20] and chloride flux into neurons (note that due to bicarbonate participation, the reversal potential of the GABAergic current (VGABA) is more positive than that for chloride; thus, even if the membrane voltage converges to VGABA, the chloride is still transported into the neurons). The increase in intracellular chloride activates the KCC2 cotransporter, which carries the potassium out of the cells [21]. Thus, the extracellular potassium concentration increases at each IID, which was proven experimentally (see Fig 12 from [5] for example).

Different types of IIDs were reported in the literature. Some types are based on pure recurrent GABAergic interactions between interneurons [7,22]. In this case, interneurons excite each other due to the depolarizing GABA effect; however, they evoke only subthreshold voltage fluctuations in the principal cells. Other types of IIDs indicate that the IIDs recruit both GABAergic and glutamatergic synaptic interactions [5]. In some cases, interneurons excite the glutamatergic principal neurons due to the depolarizing GABA effect, and then the principal neurons excite both populations [7]. In other cases, the glutamatergic neurons begin the discharges, and interneurons are involved afterward [23]. It was noted that in a rough consideration of the time scale of one-hundred milliseconds, the spiking activity of interneurons is proportional to that of pyramidal cells. It is denoted as ν(t) in the model equations. Thus, the chloride accumulation and subsequent potassium elevation by the potassium-chloride cotransporters is roughly proportional to ν(t), as expressed by the final term in Eq 1. Another contribution to the final term in Eq 1 is the potassium voltage-gated channels. The total current is also roughly proportional to the spiking activity. Also, glutamatergic currents contribute to the potassium flux [24]. Again, they are roughly proportional to the presynaptic firing rate and thus to ν(t). All these contributions are efficiently taken into account by the final term in Eq 1.

Between IIDs, the potassium concentration [K]o tends to return to the bath concentration due to diffusion and glial buffering. These processes are expressed by the first term in the right-hand side (r.h.s) of Eq 1, and the characteristic time is τK. The potassium is pumped into the neurons by the ATP-dependent Na-K pump. In Eq 1, this outflux affects the extracellular potassium concentration [K]o through the second term in the r.h.s.

Eq 2. The intracellular sodium concentration [Na]i increases due to the firing activity ν(t), as expressed by the final term in Eq 1. The same term implies a sodium flux through glutamatergic synapses [25], whose activation is roughly proportional to the firing rate ν(t). The sodium is pumped from the neurons by the ATP-dependent Na-K pump. In Eq 2, this outflux affects [Na]i through the second term in the r.h.s. The leakage and intracellular diffusion are taken into account by the first term in the r.h.s. of Eq 2. The characteristic time of the joint effect is evaluated by τNa.

Eq 3. Describes a mean membrane depolarization due to the input u(t) using a single-compartment leaky neuron model. Not taking into account a spike generation, the voltage V(t) reflects a nominal, extreme level of membrane polarization. Together with Eq 6, it is a stochastic ordinary differential equation, which corresponds to the Ornstein-Uhlenbeck process that provides colored Gaussian noise with the temporal correlation scale τm [26].

Eq 4. Describes the short-term synaptic depression according to the Tsodyks-Markram model [27] written in the rate-dependent form [28].

Eq 5. Describes the sigmoid-shaped dependence of the firing rate on depolarization V(t) (Fig 2). The dependence is not strictly valid for the transient states. Thus, it does not describe synchronization of neurons within a population and therefore not reproduce sharp peaks of the firing rate that are caused by this effect; however, it satisfactorily approximates activity averaged over a few tens of milliseconds [29]. A synaptic conductance is also missing in the relation. Thus its shunting effect is neglected. Thus, inhibition is taken into account implicitly as a negative term in the expression for u(t) (Eq 6).

Eq 6. For the neuronal input current u(t), Eq 6 takes the depolarizing effect of the potassium concentration increase relative to its initial concentration (term 1), the total synaptic current (term 2), and the noise (term 3) into account. The current u(t) that determines the firing rate by Eq 5 does not contain intrinsic voltage-gated currents but reflects the indirect effects of the ionic concentrations on membrane polarization through the voltage-gated and leak channels. Term 1 is a linear approximation of the dependence of a sum of potassium channels on their reversal potential deflection . The potassium leak current provides the main contribution to this term. Thus, the coefficient of proportionality gK,leak is called the potassium leak conductance. At the same time, the chloride leak current, which is normally smaller than the potassium leak current [30], also contributes to this term. This is in accordance with the previous assumption on the cooperative change of intracellular chloride and extracellular potassium concentrations. The change in the chloride leak current depends on the change in the chloride reversal potential, which is roughly proportional to , thus contributing to term 1.

In comparison with term 1, which evaluates the effects of potassium and chloride ions, the effect of the sodium concentration change on membrane polarization is small because the dependence of the resting potential on the sodium concentration (through the leak sodium current) is at least twice lower [30]. Therefore, the effect of sodium on u(t) was not considered.

The expression for the synaptic current (term 2) implies additional assumptions. It consists of two parts, a positive Gsyn ν(t) xD(t) and negative −0.5 Gsyn ν(t). The positive component of term 2 is an excitatory current; its approximation accounts for the synaptic resource xD and its dependence on the firing rate ν(t), which implies that synaptic kinetics are instantaneous and assumes the driving force to be constant for the sake of simplicity. Hence, the coefficient Gsyn is the maximum postsynaptic charge in response to a single presynaptic spike. The critical role of the short-term depression of excitatory synapses in the termination of IIDs has been revealed for glutamate-based [23] and depolarized GABA-based IIDs [31]. The role of short-term plasticity in inhibitory synapses does not seem to be essential. Thus, the negative component of term 2 is inhibition, which is non-depressing and proportional to the excitation with a coefficient 0.5. The value of the factor is roughly estimated to be less than 1 to provide self-sustaining recurrent excitation at the beginning of discharges when the synaptic resource is full. Below, the proportionality and the estimation of the factor are validated based on the experimental results.

The remainder of the equations, Eqs 7 and 8, are canonical [9]. Following the works by Bazhenov et al. [30] and Krishnan and Bazhenov [32], we assume that the sodium pump approximation elaborated for a single neuron [9] is valid for the description of the population level ion dynamics.

Simulations

After applying the basic set of parameter values and manipulating one of the parameters of the Epileptor-2 model, two regimes of activity with IIDs and IDs were stimulated. As explained in this section, in both regimes, the discharges are composed of spontaneous short pulses that last a fraction of a second and correspond to short bursts (SBs) of spikes in single neurons.

Regime with interictal discharges. In the case of a fast relaxation of the potassium concentration (the parameter τK = 10 s was different from the basic parameter values), the model presents a regime with spontaneous short bursts (SBs), which appear shortly after the onset of pro-epileptic bath solution with [K]bath = 8.5mM (Fig 3). In a regime of non-clustered bursts, these SBs represent IIDs. The decreased τK may correspond, for example, to the case of a more efficient extracellular diffusion of potassium ions. This diffusion equilibrates the level of the potassium concentration between IIDs, thus preventing a cluster formation.

thumbnail
Fig 3. Simulation of an interictal regime with spontaneous short bursts (SBs).

(A) The population firing rate. (B) The nominal depolarization. (C) The synaptic resource. (D) The intracellular sodium (red line) and extracellular potassium (blue) concentrations. (E) The membrane voltage in a representative neuron. (F) The membrane voltage on the time interval containing four SBs and marked by the dashed line in E. (G) The ionic flux through the Na-K-pump (orange, right axis) and the ionic concentrations (left axis) zoomed from D. Population variables in A-D and G were calculated with Eqs 18 with the basic parameter set except one different value, τK = 10s. The neuron membrane potential in E, F was obtained from Eqs 9 and 10. Note a lack of a strong clustering of SBs.

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

The mechanism of the IIDs is shown in Fig 1B. The increase of [K]o from the initial level towards the bath concentration due to diffusion leads to a u(t) increase. The generation of SB begins when the input current u(t) leads the depolarization V(t) close to the threshold Vth gL. According to Eq 5, the firing rate increases and recurrently affects u(t) in accordance to Eq 6. The synaptic resource xD(t) begins to vanish, the current u(t) decreases, and SB terminates. During each SB, a representative neuron generates a few spikes (Fig 3F). The spike generation terminates due to the decrease of u(t).

The ionic concentrations [K]o and [Na]i increase at each SB and then relax during the interburst intervals (Fig 3G). The mean levels of the ion concentrations remain approximately constant. Note that the regime with only IIDs is obtained with relatively fast [K]o relaxation (τK = 10s) compared with the characteristic interburst interval. This fast relaxation prevents any significant potassium accumulation during the train of SBs.

Regime with ictal discharges. The simulations of the model based on Eqs (18) with the basic parameter set show a population activity with recursive spontaneous IDs (Fig 4A), shortly after the onset of pro-epileptic bath solution with [K]bath = 8.5mM. Each ID is characterized by a high rate of activity (black line) for 30 seconds. It consists of SBs resembling IIDs (Fig 4B and Fig 5), and each lasts hundreds of milliseconds. The intervals between IDs are about two minutes. The mechanism of IDs is given schematically in Fig 1B. The concentration [K]o increases towards the level [K]bath. At some critical level of [K]o, the network begins to generate SBs. [K]o, the extracellular potassium concentration, increases during each SB and partially relaxes between SBs. The extracellular potassium accumulates (Fig 4A and 4C, blue line). The SB frequency rapidly grows with increasing [K]o, forming a cluster (Fig 4B and 4C). The intensified firing leads to an accumulation of sodium inside neurons (Fig 4A and 4C, red line). The Na-K-pump function strongly depends on [Na]i; therefore, the pump activity rapidly increases (Fig 4A and 4C, orange line). At some point, the pump dominates the outflux of potassium; thus, [K]o begins to decrease. The decline of [K]o decreases the frequency of SBs. The sodium concentration reaches its maximum and begins to restore. While [Na]i is high, the Na-K-pump is highly active, which maintains the reduction of potassium concentration. At some level of [K]o, the SB generation terminates. During the next cycle, towards the next ID, [K]o begins to increase (Fig 4A) due to the diffusion with the bath solution, which dominates at this stage.

thumbnail
Fig 4. Simulations of IDs, each consisting of separate short bursts.

(A) The nominal depolarization (top plot), the intracellular sodium (red line) and extracellular potassium (blue) concentrations (bottom, left axis), and the ionic flux through the Na-K-pump (orange line, bottom plot, right axis) during six IDs. (B) The population firing rate (black), the nominal depolarization (red), and the synaptic resource (violet) during a single ID consisting of some SBs. (C) Zoomed traces from A, bottom during a single ID. Simulations were done with Eqs 18 with the basic parameter set. Note the decrease of [K]o at the high level of peaks of Ipump and following the termination of the ID.

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

thumbnail
Fig 5. Simulations of a single neuron activity during the IDs shown in Fig 4.

(A) Two IDs as bursts of clustered SBs, seen in the membrane voltage. (B) A single ID containing a number of SBs. Membrane voltage. (C) A single SB.

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

The regime with IDs is obtained with relatively slow [K]o relaxation (τK = 100s, in accordance with the basic parameter set) compared with the characteristic interburst interval. This slow relaxation leads to potassium accumulation at some point of the train of SBs. At some critical level of [K]o, depolarization due to the shift of VK(t) (term 1 in Eq 6) overcomes the decrease of the excitation due to the decrease of the synaptic resource xD(t). It begins an ID via the positive feedback provided by the depolarizing effect of the potassium accumulation.

The results of this subsection lead to the conclusion that Epileptor-2 is able to describe both regimes with interictal and ictal discharges. The mechanisms of the discharge generation correspond to those shown in the schematic in Fig 1B. Simulations illustrate crucial roles of each of the processes included in the schematic (S1 Text, section 1). Modifications of Epileptor-2, given in S1 Text, sections 2–4, extend the model to a case of epileptiform activity in a disinhibited network, a case of prominent depolarization block or consideration of more realistic, adaptive neuron-observer.

Comparison of the model with the experiments

The IDs simulated using the proposed model are similar to the those registered in slices of a rat entorhinal cortex (Fig 6). The ictal events occurred with a period of about three minutes (compare Fig 6A with Fig 5A). Each ID consisted of SBs (Fig 6B with Fig 5B). Each SB consisted of several spikes (Fig 6C with Fig 5C). The shift of the baseline of the potential during an ID was observed during the experiment (Fig 6A) and in the average potential V(t) in the model (Fig 4A); it is not present in U(t) due to the approximate nature of the single neuron model used.

thumbnail
Fig 6. Ictal discharges recorded in a pyramidal neuron from a rat entorhinal cortex using an in vitro 4-AP model of epileptic activity.

(A) The full record. (B, C) Zoomed fragments. Note each ID is a burst of SBs of spikes.

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

Ionic dynamics. The transmembrane concentration gradient of the six major ions (K+, Na+, Cl-, Ca2+, H+, and HCO3-) is altered during an epileptic seizure [15]. The most well-studied are the levels of extracellular K+ and Ca2+ and intracellular Cl-, Na+, and H+ concentrations [15]. Changes in the concentration of extracellular K+ and intracellular Cl- and Na+ play the most critical roles in seizure induction and maintenance [15]. The central assumptions of the Epileptor-2 assert that only extracellular K+ and intracellular Na+ concentrations govern the activity, thus excluding intracellular Cl- from consideration. This assumption implies a proportionality between the changes of [K]o and [Cl]in. Fig 7 was reproduced from the study by Raimondo [15] and supports this assumption. The profiles of [K]o and [Cl]in during ID are similar. Presumably, this proportionality is explained by the activity of K-Cl-cotransporters, which evoke a potassium outflux in response to chloride accumulation within the neurons.

thumbnail
Fig 7. Ion concentration changes during seizures.

Modified with permission from [15].

https://doi.org/10.1371/journal.pcbi.1006186.g007

The extracellular potassium concentration increases during a single SB up to 1–2 mM according to data from [33] (see their Fig 7 or Fig 12 in [5]) and [8,14]. This potassium flux originates from both pyramidal cells and interneurons [8]. The potassium voltage-gated channels, glutamatergic receptors, and K-Cl cotransporters contribute to the flux. The increase of [K]o during a single SB results in moderate depolarization on about 1 mV according to Eq 7. In the Epileptor-2, the value of parameter δ[K]o provides the same level of potassium increase after each SB; the parameter τK defines the speed of potassium level relaxation, and the parameter gK,leak controls the level of polarization dependent on [K]o.

According to the experimental observations [5,8,14,15], the extracellular potassium concentration has a peak in the middle of an ID, while the intracellular sodium concentration still gradually increases until the end of an ID [15] (Fig 7). The model reproduces the ionic dynamics (Fig 4A), employing the nonlinear activation of the Na-K-pump. The pump is activated by a high sodium concentration (Fig 4C) according to the approximation (Eq 8). The influx of potassium through the pump overcomes the outflux, and [K]o begins to decrease. Note that the peak of [K]o precedes the decrease of [Na]i in the simulation, the mentioned experiment (Fig 7), and the previous modeling [9]. The decrease of [K]o results in the termination of the ID, as shown by the experiments. The model explains this ID termination by the decrease of the potassium-dependent current (term 1 in the r.h.s. of Eq 6). This suggests a crucial role of the Na-K pump in ID termination.

Indirect support of the hypothesis of the crucial role of the ATP-dependent Na-K-pump in ID termination is provided by the evidence of a drop in the oxygen concentration during an ID, which reflects the activity of energy metabolism [1].

Synaptic currents. The second term in Eq 6 implies a proportionality between the excitatory and inhibitory currents and the firing rate. The experiments support this approximation. Assuming that two neighboring pyramidal neurons receive roughly identical synaptic input during epileptiform events [7], the pairing recordings were made from the neighboring neurons during IDs and IIDs. The excitatory and inhibitory currents from one neuron in a voltage-clamp mode and detected spikes in another neuron in a cell-attached mode were recorded (Fig 8A and 8B). The holding voltages in the voltage clamp corresponded to the reversal potentials of GABAergic and glutamatergic currents, -50 mV and 0 mV, respectively. The current profiles for several SBs were superimposed (Fig 8C). The mean currents were rescaled to correspond to the voltage level -40 mV, which is close to the spike threshold. These mean inhibitory and inverted excitatory currents are shown in Fig 8D. The instantaneous firing rate was calculated from cell-attached recordings as an average number of spikes per interval divided by the length of the interval (12.5 ms). It was found that (i) the shapes of the currents and the instantaneous firing rate are roughly similar; (ii) at spike threshold potential, the excitatory current is bigger than the inhibitory one; (iii) at the initial stage, the peak of the excitatory current is twice larger than that of the inhibitory current; and (iv) the excitatory current decreases faster than the inhibitory one. These findings support the assumptions of the proportionality of the currents to the firing rate, the effect of short-term synaptic depression, and the ratio of excitatory versus inhibitory currents before the action of depression. Therefore, these findings (observations) validate term 2 of Eq 6.

thumbnail
Fig 8. Proportionality between the excitatory and inhibitory currents and the firing rate in the experiment.

(A) Upper trace: A representative recording of glutamatergic synaptic activity in the voltage-clamp mode at the reversal potential of GABA-mediated currents (Vhold = VGABA = −51mV). Lower trace: A corresponding firing activity of the neighboring neuron recorded in the cell-attached voltage-clamp mode. (B) The recordings from the same pair of neurons as in A. The upper trace was recorded at the reversal potential of glutamate-mediated currents (Vhold = VGlut = 0) and represents the GABAergic synaptic activity. The yellow shadow indicates the IDs. (C) Several superimposed IID-like events from A (left) and B (right) that correspond to SBs and constitute IDs. (D) Pooled data from four IDs, each consisting of 50–80 SBs as in C. A histogram represents the firing rate of the neuron recorded in the cell-attached VC mode. Solid lines are the plots of GABA- and glutamate-mediated currents scaled to represent the synaptic input at the membrane potential of -40 mV.

https://doi.org/10.1371/journal.pcbi.1006186.g008

Mathematical analysis of SBs

Based on the simulations of the system of Eqs 18 (Fig 4), the variables xD(t) and V(t) oscillate faster than the dynamics of [K]O and [Na]i, producing the SBs. Thus, the subsystem of Eqs 35 was studied under the conditions of fixed ionic concentrations, i.e., with . The system of Eqs 35 is planar and stochastic with non-smooth r.h.s., as follows: (3’) (4’) (5’)

Nullclines of the correspondent deterministic system (σ = 0) are shown in the phase space (Fig 9A). The deterministic system has a unique stable equilibrium point with V = 0 and xD = 1. There are no other equilibrium points or limit circles. Thus, all possible trajectories converge to this stable steady state. Consider a trajectory with an initial position inside or to the left of the U-shaped nullcline. The trajectory first approaches the U-shaped nullcline, then reaches the bottom point, separates towards the straight line nullcline, and proceeds towards the steady state. Because a voltage trace that corresponds to this trajectory would be shaped as an impulse, the zones inside and to the left of the U-shaped nullcline could be referred to as zones of excitation.

thumbnail
Fig 9. Fast subsystem producing SBs.

(A) The phase space of the system of Eqs 3’5’. Nullclines (solid red lines) and vector-fields (arrows) are calculated for ; dashed nullclines are for . Purple curve is a trajectory. (B) the traces of V(t) and xD(t). The parameter values correspond to the basic set of parameters, except the noise amplitude σ/gL = 40 mV.

https://doi.org/10.1371/journal.pcbi.1006186.g009

Next, the effect of random disturbances on the system (3’-5’) that has the only attractor of the deterministic system was examined. A small amplitude noise that leads to small voltage fluctuations near 0 cannot excite the system. Indeed, small deviations from the equilibrium result in the monotonic convergence to the equilibrium; however, deviations larger than the threshold result in a large excursion before returning to the resting state. This implies a stochastic excitability. As the noise intensity increases, the time between two successive activations decreases. A random trajectory beginning from the stable equilibrium and declined by the noise to the zone of excitation is plotted in the phase plane (Fig 9A) along with corresponding time series (Fig 9B). When the noise intensity is higher than some threshold value, random trajectories can move far from the stable equilibrium, and large-amplitude oscillations occur along with small-amplitude oscillations. The intermittency of the oscillations forms the stochastic bursting.

Finally, the effect of the extracellular potassium concentration on the frequency and magnitude of SBs was examined. The level of [K]O determines an additive current in Eq 3 according to Eqs 6 and 7. An increase in [K]O leads to an inward current. The inward current shifts the U-shaped nullcline left and down and shifts the straight line nullcline to the right (dashed lines in Fig 9A). Consequently, the probability of the excursion in the phase-space to cross the nullcline from left to right is higher. The amplitude of oscillations decreases, and the frequency increases.

The noise brings a new quality to the system, which cannot be observed in the deterministic case. This phenomenon, known as coherence resonance, is a generic feature of noise-driven excitable systems and is mostly independent of the particular details of the model [16,34,35]. For epileptic discharges, the mechanism of stochastic large-amplitude oscillations underlies the SBs observed as interictal events and as primary components of IDs.

Mathematical analysis of IDs as slow oscillations

To analyze an appearance of IDs, the full model was reduced to a slow subsystem. In the simulations (Fig 4), the dynamics of [K]o and [Na]i consisted of slow and low-amplitude fluctuations. The former are much slower than the dynamics of xD, V, and ν. Also, it was noted that Eqs 1 and 2 have slower time scales (τK and τNa) than Eqs 3 and 4 with the time scales τm and τD. These observations allowed for reducing the full model to a slow subsystem that is based on Eqs 1 and 2. The fast subsystem based on Eqs 3 and 4 and discussed in the previous subsection is controlled by the only slow variable [K]o. On the other hand, Eqs 1 and 2 are affected by the fast subsystem only through ν(t). In the simulations, the slow high-amplitude oscillations of [K]o and [Na]i did not follow the fast fluctuations of ν but were controlled by a slow component of ν. Therefore, to extract a slow subsystem, ν(t) was averaged, and the dependence of this average variable on [K]o was evaluated. This was done with the integration of Eqs 37 for the slow ramp of [K]bath (from 3 to 22 mM during 200s) with the other parameters taken from the basic parameter set. The obtained dependence was approximated as follows: (5”) which is valid for [K]o < 20 mM.

After the substitution of ν by in Eqs 1, 2 and 8, the slow subsystem is obtained: (1’) (2’) (8’) where [K]bath = 8.5 mM, , τK = 100 s, τNa = 20 s, δ[K]o = 0.02 mM, δ[Na]i = 0.03 mM, γ = 10, ρ = 0.2 mM/s, as in the basic parameter set.

The reduced model based on Eqs 1’, 2’, 5” and 8’ behaves similarly to the full Epileptor-2, as shown in Fig 10.

thumbnail
Fig 10. Slow subsystem producing IDs.

(A) IDs in the full Epileptor-2 model. (B) Oscillations in the reduced model based on Eqs 1’, 2’, 5” and 8’. The plots repeat those in Fig 4. (C, D, and E) Phase portraits of the reduced model obtained for [K]bath = 3, 6.42017, and 8.5 mM, correspondingly. Note the stable node in C and the limit cycle around the unsteady focus in D.

https://doi.org/10.1371/journal.pcbi.1006186.g010

The slow dynamical system was then analyzed based on Eqs 1’, 2’, 5” and 8’. Under “normal” conditions of a low potassium concentration in the bath solution, [K]bath = 3 mM, the system converges to a fixed point, which is a stable node. The other two equilibrium points are an unstable focus and a saddle. At some critical value of , the stable node couples with the saddle at the boundary that is given by a kink of the function , Eq 5”. This bifurcation is a codimension-two bifurcation in a piecewise-smooth system [36]. To the best of the authors’ knowledge, a full classification of the codimension-two bifurcations of non-smooth systems is still challenging, even in the case of planar systems. Nevertheless, the bifurcation of this system is similar to SNIC (saddle-node on invariant cycle) bifurcation, which can be observed in smooth systems. In this bifurcation, a limit cycle originates around the unsteady focus when the node and saddle disappear. The limit cycle has a zero initial frequency and a finite amplitude. This bifurcation corresponds to the transition to ID generation in the Epileptor-2 model. In the pathological condition of , such as in the case of [K]bath = 8.5 mM, as shown in Fig 4, a stable oscillation of [K]o is observed, which controls the generation of SBs during IDs. Overall, the analysis explains the dynamics of the regime with ID generation due to oscillations of potassium and sodium concentrations, with each ID composed of clustered SBs.

Discussion

Model assumptions

In the present work, a mathematical model of IDs and IIDs is proposed based on a consideration of only a few dominating processes that underlie the generation of the discharges. These processes have been described in the form of a simple mathematical model consisting of four ordinary differential equations written in terms of physically meaningful variables, referred to as Epileptor-2, as an alternative to the well-known but more abstract model Epileptor [1]. Along with a modeled neuron-observer, simulations of Epileptor-2 provide voltage, firing rate, and ionic concentration traces that are close to experimental recordings in many aspects. The model has been derived under the strong assumptions described in the Results section. The most important assumptions are: (i) ion concentrations that affect the discharges are the extracellular potassium and the intracellular sodium concentrations (Fig 1A); (ii) both excitatory and inhibitory currents are mainly controlled by the firing rate of an integrated neuronal population; (iii) only potassium concentration affects the intrinsic neuronal excitability; (iv) the excitatory and inhibitory currents are controlled by the same presynaptic firing rate; and (v) the excitatory current undergoes short-term depression.

The first assumption has been partially validated by experimental observations from the literature [15]. It is also consistent with that from modeling studies [10,11,32]. For instance, it was concluded in the detailed modeling study that it is the coupled dynamics of sodium, potassium, and chloride ions play a critical role in the development and termination of seizures [32]. The main ground for this assumption was relatively high concentrations of the complementary intracellular potassium, and extracellular sodium and chloride ions that are less affected by the ion dynamics and less affect the reversal potentials. Extracellular calcium concentration is shown to be progressively reduced during seizure [36], which decreases the reliability of synaptic transmission and can affect long-range synchronization properties of cell firing during ID [37]. Nevertheless, the mechanism of ID termination explored in this study depends on the specific cellular properties, such as activation of the Na-K pump. This pump is driven by ATPase and represents an intrinsic factor that is not directly affected by the extracellular calcium [32]. The intracellular calcium is connected to the neuronal excitation and short-term synaptic depression. Thus it implicitly affects the input-output function, Eq 5, and the parameters of Eq 4 for the synaptic resource. Based on these arguments the calcium ions are excluded from explicit consideration. In addition, our study associates the intracellular chloride concentration changes with the extracellular potassium, thus keeping only the extracellular potassium and the intracellular sodium concentrations.

The second assumption seems to be too strong if taking the sometimes crucial role of GABAergic interneurons in interictal discharges into account [6,7,22,31,38,39]; however, the types of the discharges studied more often include strong glutamatergic excitation [5,40] or are only based on glutamatergic interactions [23]. Meanwhile, even GABAergic channel-based IIDs are in agreement with Epileptor-2 if they are considered contributors to excitatory but not inhibitory currents, which are described by the second term in the r.h.s. of Eq 6. Indeed, GABA-based excitation implies a depolarized GABA receptor reversal potential. In this case, quantitatively, these receptors contribute to the positive term. The synapses that maintain a healthy level of chloride contribute to the negative term. Interneurons that are presynaptic to both subpopulations of synapses are active simultaneously. Thus, both synaptic currents roughly correlate with the firing rate of the excitatory population, as supposed by the fifth assumption. These arguments extend the applicability of the Epileptor-2 to the case of interneuron-driven discharges.

The third assumption postulates an essential role of the extracellular potassium concentration. It corresponds to the Fertziger and Ranck hypothesis [41]; however, the intracellular chloride concentration is also a potent regulator of neuronal excitability. These two variables are tightly connected, as the accumulation of intracellular chloride due to the activation of GABAergic receptors leads to an increase in KCC2-activity and a consecutive increase in extracellular potassium levels, which in turn directly affects the excitability and results in a less-efficient evacuation of chloride ions from the cell. The latter might at some point result in a failure of inhibition, causing more dramatic effects on neuronal excitability. In some experimental models of epilepsy, the activity of GABAergic interneurons is the primary cause of extracellular potassium build up and the progressive increase of neuronal excitability [5,6]; however, in this model, no particular cell types were introduced, and the firing rate was assumed to be the same for all cell types. A feasible approximation in this case is an assumption of the proportionality between the elevation of extracellular potassium and intracellular chloride concentrations, which allows for accounting for only one type of ion. Thus, in this model, several steps that link firing activity to potassium accumulation were skipped, and a simplified approach where cell firing directly affects the potassium concentration was implemented. Moreover, the change in the level of potassium affects neuronal depolarization via the potassium reversal potential [42]. Therefore, for simplicity, only the potassium concentration was included in the model, and the proportionality between intracellular chloride and extracellular potassium concentrations was assumed. Based on this assumption, the effects of their temporary imbalance [43] was neglected. Consequently, this approach is entirely different from that by Proix et al. [44], who introduced a permittivity variable that captures effects evolving on slow timescales and including extracellular ionic concentrations and energy metabolism. The permittivity variable controls the remainder of the dynamical system of the Epileptor. The advantage of the new approach is the explicit form of the dependence on the potassium concentration.

The fourth assumption has been validated based on the experimental data. More precise matching would require explicit consideration of excitatory and inhibitory populations synaptic components to fit experimental estimates of either cumulative excitatory and inhibitory synaptic conductances [45,46] or AMPA, NMDA and GABA-mediated conductances, as done in [31] for interictal discharges. However, such detailed modeling in the case of ID generation would involve a complex description of the ionic dynamics. Thus such modeling is on the opposite extreme of detail in comparison to the proposed Epileptor-2 model. Another critical point to the fourth assumption is its validity in the case of the depolarization block, which occurs in interneurons during IDs [1,47]. Therefore, the approximation might overestimate the inhibitory currents as well as the potassium accumulation mediated by the accumulation of chloride and the activity of K-Cl-cotransporters during some period within an ID. These effects are expected to influence the duration of IDs. This view is supported by a modified Epileptor-2 model that takes into account the depolarization block, as described in S1 Text, section 3.

The fifth assumption is based on experimental evidences of a short-term depression of glutamatergic synapses [23,48,49]. An enhanced depression of GABAergic synaptic transmission was also found to be a crucial factor for reproducing synchronized synaptic activity in the previous model of interictal discharges. Experimental observations partly confirmed this assumption, as GABA-mediated responses to extracellular stimuli seemed to be occluded [31]. The synaptic depression of inhibitory synapses has not been included in the model for simplicity because its contribution could only reshape the discharges.

In summary, the assumptions of the Epileptor-2 model are in line with a wide range of experimental evidence obtained in different models of epilepsy.

Stochastic oscillations as a mechanism of epileptic discharges

Based on the simulations, both ictal and interictal discharges consist of elementary events, SBs. To the best of the authors’ knowledge, the IDs have not been modeled before as bursts of spike bursts. The SBs have been mathematically analyzed. They were found to be large amplitude stochastic oscillations. The second-order system of differential equations is probably the simplest model of these kernel elements of pathological discharges. This statement is an important prediction, which opposes the inherently stochastic mechanism of the discharge generation to an oscillatory mechanism. The latter explains a stochastic sequence of the discharges due to noise in a deterministic periodic process. In contrast, the former implies that the discharges are impossible without noisy fluctuations. Similar behaviors have been studied for a neuronal model [16].

Bifurcations leading to epileptic discharges

The original Epileptor model explains epileptic discharges in terms of bifurcations [1]. In this mode, the types of bifurcations that correspond to the onset and offset of IDs are saddle-node and homoclinic bifurcations, respectively [1] with the control parameter of “permittivity,” which is the slowest variable. The saddle-node bifurcation at the onset of IDs in the Epileptor was chosen based on experimentally observed features, such as fixed frequency and fixed amplitude of abruptly started oscillations, and a shift of baseline field potential. The homoclinic bifurcation at the offset of IDs was chosen for the shift of the baseline, fixed amplitude, and decaying frequency of oscillations.

In contrast to the permittivity parameter in the Epileptor, this model has [K]O as the control variable of fast oscillations, SBs. The variable [K]O shows slow oscillations with pulses (Fig 10A and 10B). Due to these oscillations, the frequency of SBs oscillates as well, forming IDs; however, no bifurcation occurs in the fast, stochastic system. Instead, the frequency of SBs changes along with the probability of spontaneous SB generation, depending on the level of [K]O, as explained with the mathematical analysis. The three characteristic features of the onset of ID are similar to the ones in the original Epileptor: (i) the amplitude of the fast oscillations is large (Fig 9); (ii) their frequency grows not abruptly but rapidly (Fig 4B) due to the rapid onset of the pulses of [K]O; and the baseline shift of voltage (Fig 4A) occurs due to the depolarizing effect of the high level of [K]O. The offset of ID corresponds to the decaying phase of a slow pulse of [K]O with a rate of frequency decrease similar to that at the onset of ID. Thus, to a large extent, the Epileptor-2 reproduces the main features of IDs that were captured by the original Epileptor.

The analysis explains the scenario of an origination of a regime of ID generation during a change in an external factor, an increase of [K]bath, similar to previous studies conducted with a single neuron model [9,10]. Using a fast-slow analysis, it was found that the scenario corresponds to the SNIC-like bifurcation. This implies that: (i) IDs are seldom at the levels of [K]bath just above the critical value and (ii) the amplitude of [K]O oscillations is always finite, which indicates a similar shape of IDs at different [K]bath. Both predictions are more or less consistent with the experimental observations. The whole system of the Epileptor is of fold/homoclinic type [50].

The Epileptor-2 can not be classified in the same way. Instead, fast and slow dynamics are to be considered separately. In contrast to the Epileptor, the Epileptor-2 is a non-smooth, stochastic dynamical system. The non-smoothness is due to the kink in Eq 5”. This kink reflects the threshold nature of neuronal excitation. This is consistent with the observations of the rapid development of spontaneous spiking in a slice that undergoes a change in potassium concentration from low to high values in a bath. The presence of noise in the Epileptor-2 is also quite natural and realistic. The voltage fluctuation in normal and pathological states and an irregularity of IIDs indicate the intrinsic stochasticity of the network, which validates the predictions of the Epileptor-2 model.

Several mechanisms of seizure recruitment have been tested in computational models with varying levels of mathematical abstraction [32,52,53]. Typically, the most complex, spatially distributed, stochastic models are studied by means of their reduction to the system of ordinary differential equations. For example, the model proposed by Steyn-Ross et al [54] has been modified and studied by Kramer et al. [51,55]. Their main variable is compared to the experimentally registered local field potential. The oscillating regime for this variable is regarded as an ictal state. The transitions to and from the ictal state have been studied in a plane of two parameters of excitation. The obtained bifurcations help to understand the behavior of the full spatially distributed, stochastic system. In comparison with the Epileptor-2, the ionic dynamics is not considered in the Kramer’s model, however the excitation variables might be related to the ionic concentrations. In the regime of ID generation in Epileptor-2, the concentrations oscillate. Therefore, the transitions in the Kramer’s model are to be compared to the gradual change of SB generation in Epileptor-2, from the rare bursting between IDs to the fast bursting during each of the IDs. In Epileptor-2, there are no bifurcations. Instead, the probability of the SB generation gradually changes, with SBs being described by the fast subsystem and always representing stochastic large-amplitude oscillations. In this sense, our model is qualitatively different from the previous models [32,43,52,53]. Meanwhile, in the Epileptor-2, the slow dynamics of the concentrations undergoes bifurcations of the transitions to and from the oscillations governed by one of the external parameters of excitation, like the bath potassium concentration that controls the SNIC-like bifurcation.

Conclusion

In this paper, a simple model of epileptic discharges is proposed, and the fundamental mechanisms described by this model are illustrated. The model is qualitatively and in a sense quantitatively compared to a wide range of experimental recordings. Similar to the known Epileptor model, this model can be applied to investigate the pathological states in a network of coupled oscillators [44]. It can also be used as a navigator in the parameter space for biophysically detailed modeling.

Methods

Ethics statement

All animal procedures followed the guidelines of the European Community Council Directive 86/609/EEC and were approved by the Animal Care and Use Committee of the Sechenov Institute of Evolutionary Physiology and Biochemistry of the Russian Academy of Sciences.

Slice preparation and electrophysiological recordings

Experimental procedures were described previously in details in our recent paper [7]. Shortly, 3-week-old Wistar rats were decapitated and their brains removed rapidly. A vibrating microtome (Microm HM 650 V; Microm; Germany) was used to cut 300-μm thick horizontal slices that contained entorhinal cortex and hippocampus. Artificial cerebrospinal fluid (ACSF) was used through all steps, which had the following composition (in mM): 126 NaCl, 2.5 KCl, 1.25 NaH2PO4, 1 MgSO4, 2 CaCl2, 24 NaHCO3, and 10 dextrose. The ACSF was aerated with carbogen (95% O2/5% CO2). Recordings were made at 30˚C. Pyramidal neurons in deep layers of the entorhinal cortex were visualized using a Zeiss Axioscop 2 microscope (Zeiss; Oberkochen, Germany) equipped with a video camera (PointGrey Grasshopper3 GS3-U3-23S6M-C, FLIR Integrated Imaging Solutions Inc., USA) and differential interference contrast optics. Patch electrodes (3–5 MΩ) were pulled from borosilicate filamented glass capillaries (WPI; UK) on a P-1000 Micropipette Puller (Sutter Instrument; Novato, CA, USA). For current-clamp recordings, a potassium-gluconate-based filling solution was used with the following composition (in mM): 135 K-gluconate, 10 NaCl, 5 EGTA, 10 HEPES, 4 ATP-Mg, and 0.3 GTP (with pH adjusted to 7.25 with KOH). For voltage-clamp recordings, a solution based on cesium-methane-sulfonate (CsMeS) was used with the following composition (in mM): 127 CsMeS, 10 NaCl, 5 EGTA, 10 HEPES, 6 QX314, 4 ATP-Mg, and 0.3 GTP (with pH adjusted to 7.25 with CsOH). For cell-attached voltage-clamp recordings, a sodium-chloride-based filling solution was used. It had the following composition (in mM): 138.5 NaCl, 8.5 KCl, 10 HEPES, 5 EGTA (with pH adjusted to 7.25 with NaOH). Cell-attached voltage-clamp recording of neuron firing activity was performed as described [56].

Recordings were performed with two Model 2400 patch-clamp amplifiers (AM-Systems; Sequim, WA, USA), and an NI USB-6343A/D converter (National Instruments; Austin, TX, USA) using WinWCP5 software (SIPBS; Glasgow, UK). The data were filtered at 10 kHz and sampled at 20 kHz. After formation of the whole-cell configuration, access resistance was less than 15 MΩ and remained stable (< 30% increase) during the experiments in all cells included.

Epileptiform activity was induced with the pro-epileptic solution, containing the following (in mM): 120 NaCl, 8.5 KCl, 1.25 NaH2PO4, 0.25 MgSO4, 2 CaCl2, 24 NaHCO3, 10 dextrose, and 0.05 4-AP. The flow rate in the perfusion chamber was 5–6 ml/min. The liquid junction potentials were measured as described [57], and the holding potential was compensated offline for whole-cell voltage-clamp recordings by subtracting 7 mV.

Modeling and analysis

The new model is based on the previous biophysically detailed considerations of ionic dynamics during the pathological states of brain activity [9,17,30,42]. The ionic dynamics were implemented in a rate-based model for recurrently connected excitatory and inhibitory neuronal populations, where the inhibitory population has been accounted for implicitly, and the firing rate was assumed to be proportional to that of the excitatory population. The firing rate has been described as a rectified sigmoid function of a membrane potential (Fig 2). The membrane potential has been described by Kirchoff’s current conservation law, which was written for a one-compartment neuron. The expressions for the excitatory and inhibitory synaptic currents, the input-output-function, the rate-based equations for the ionic dynamics, etc., have been justified in the Results section. The short-term synaptic depression is described according to the Tsodyks-Markram model [27]. A quadratic integrate-and-fire model was used as a model for a representative neuron [19].

The simulations were performed in the Delphi-7 environment. The mathematical analysis of the stochastic oscillations was performed using Wolfram Mathematica 10 (Champaign, IL, USA). The Euler-Maruyama explicit numerical scheme was applied for the integration of the stochastic ordinary differential equations. The typical value of a time step was 0.5 ms. The results were dependent on the numerical parameter in a similar extent as for different realizations of noise. The numerical realizations of the model are available from the websites: the code in Wolfram Mathematica is at https://yadi.sk/d/927UjbS-3QQhMW; the code and executive file in Delphi-Pascal are at https://drive.google.com/file/d/1AJhAFKLOjvgauBF_6SQzvol8zzmfMaDV/view?usp=sharing and https://drive.google.com/open?id=10ij-Nt780jROcMv9qniUm4rAJNr8WD2j, correspondingly.

Supporting information

S1 Text. Illustrations to the mechanisms of IDs and IIDs.

Disinhibition as a model of epilepsy. Depolarization block. Alternative model of a neuron-observer.

https://doi.org/10.1371/journal.pcbi.1006186.s001

(PDF)

Acknowledgments

The authors are grateful to Elena B. Soboleva for her help in the experiments.

References

  1. 1. Jirsa VK, Stacey WC, Quilichini PP, Ivanov AI, Bernard C. On the nature of seizure dynamics. Brain. 2014; 137(8): 2210–2230.
  2. 2. Nagaraj V, Lamperski A, Netoff TI. Seizure Control in a Computational Model Using a Reinforcement Learning Stimulation Paradigm. Int J Neural Syst. 2017; 27(7):1750012. pmid:28030999
  3. 3. Naze S, Bernard C, Jirsa V. Computational modeling of seizure dynamics using coupled neuronal networks: factors shaping epileptiform activity. PLoS Comput Biol. 2015; 11(5):e1004209. pmid:25970348
  4. 4. Melozzi F, Woodman MM, Jirsa VK, Bernard C. The Virtual Mouse Brain: A Computational Neuroinformatics Platform to Study Whole Mouse Brain Dynamics. eNeuro. 2017; 4(3): pii: ENEURO.0111-17.2017.
  5. 5. Avoli M, D'Antuono M, Louvel J, Kohling R, Biagini G, Pumain R, et al. Network and pharmacological mechanisms leading to epileptiform synchronization in the limbic system in vitro. Prog Neurobiol. 2002; 68(3): 167–207. pmid:12450487
  6. 6. Yekhlef L, Breschi GL, Lagostena L, Russo G, Taverna S. Selective activation of parvalbumin- or somatostatin-expressing interneurons triggers epileptic seizurelike activity in mouse medial entorhinal cortex. J Neurophysiol. 2015; 113(5): 1616–1630. pmid:25505119
  7. 7. Amakhin DV, Ergina JL, Chizhov AV, Zaitsev AV. Synaptic Conductances during Interictal Discharges in Pyramidal Neurons of Rat Entorhinal Cortex. Front Cell Neurosci. 2016; 10: 233. pmid:27790093
  8. 8. Librizzi L, Losi G, Marcon I, Sessolo M, Scalmani P, Carmignoto G, et al. Interneuronal network activity at the onset of seizure-like events in entorhinal cortex slices. J Neuroscience. 2017; 37(43): 10398–10407.
  9. 9. Cressman JR, Ullah G, Ziburkus J, Schiff SJ, Barreto E. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. J Comput Neurosci. 2009; 26(2): 159–170. pmid:19169801
  10. 10. Owen JA, Barreto E, Cressman JR. Controlling seizure-like events by perturbing ion concentration dynamics with periodic stimulation. PLoS One. 2013; 8(9): e73820. pmid:24066075
  11. 11. Ullah G, Cressman JR Jr, Barreto E, Schiff SJ. The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states. II. Network and glial dynamics. J Comput Neurosci. 2009; 26(2): 171–183. pmid:19083088
  12. 12. Krishnan GP, Filatov G, Shilnikov A, Bazhenov M. Electrogenic properties of the Na/K ATPase control transitions between normal and pathological brain states. J Neurophysiol. 2015; 113: 3356–3374. pmid:25589588
  13. 13. Ho ECY, Truccolo W. Interaction between synaptic inhibition and glial-potassium dynamics leads to diverse seizure transition modes in biophysical models of human focal seizures. J Comput Neurosci. 2016; 41(2): 225–244. pmid:27488433
  14. 14. Antonio LL, Anderson ML, Angamo EA, Gabriel S, Klaft ZJ, Liotta A, et al. In vitro seizure like events and changes in ionic concentration. J Neurosci Methods. 2016; 260: 33–44. pmid:26300181
  15. 15. Raimondo JV, Burman RJ, Katz AA, Akerman CJ. Ion dynamics during seizures. Front Cell Neurosci. 2015; 9: 419. pmid:26539081
  16. 16. Bashkirtseva I, Fedotov S, Ryashko L, Slepukhina E. Stochastic Bifurcations and Noise-Induced Chaos in 3D Neuron Model. Int J Biffurcation and Chaos. 2016; 26(12): 1630032–1630021.
  17. 17. Kager H, Wadman WJ, Somjen GG. Simulated seizures and spreading depression in a neuron model incorporating interstitial space and ion concentrations. J Neurophysiol. 2000; 84(1): 495–512. pmid:10899222
  18. 18. Latham PE, Richmond BJ, Nelson PG, Nirenberg S. Intrinsic dynamics in neuronal networks. I. Theory. J Neurophysiol. 2000 Feb;83(2):808–27. pmid:10669496
  19. 19. Izhikevich EM. Simple Model of Spiking Neurons IEEE Transactions on Neural Networks. 2003; 14: 1569–1572. pmid:18244602
  20. 20. Kaila K, Voipio J, Paalasmaa P, Pasternack M, Deisz RA. The role of bicarbonate in GABAA receptor-mediated IPSPs of rat neocortical neurones. J Physiol. 1993; 464: 273–289. pmid:8229801
  21. 21. Payne JA, Rivera C, Voipio J, Kaila K. Cation-chloride cotransporters in neuronal communication, development and trauma. Trends Neurosci. 2003; 26: 199–206. pmid:12689771
  22. 22. Pallud J, Le Van Quyen M, Bielle F, Pellegrino C, Varlet P, Labussiere M et al. Cortical GABAergic excitation contributes to epileptic activities around human glioma. Science Translational Medicine. 2014; 6: 244ra89–244ra89. pmid:25009229
  23. 23. Staley KJ, Longacher M, Bains JS, Yee A. Presynaptic modulation of CA3 network activity. Nature Neuroscience. 1998; 1: 201–209. pmid:10195144
  24. 24. Mayer ML, Westbrook GL. Permeation and block of n-methyl-d-aspartic acid receptor channels by divalent cations in mouse cultured central neurones. J.Physiology 1987; 394: 501–527.
  25. 25. Rose CR, Konnerth A. NMDA Receptor-Mediated Na1 Signals in Spines and Dendrites. J.Neuroscience 2001; 21(12): 4207–4214.
  26. 26. Gerstner W., Kistler W. M., Naud R., and Paninski L. Introduction: neurons and mathematics. In: Neuronal Dynamics. From Single Neurons to Networks and Models of Cognition. Cambridge: University Press; 2014. pp. 395–416.
  27. 27. Tsodyks M, Pawelzik K, Markram H. Neural Networks with Dynamic Synapses. Neural Computation. 1998; 10(4): 821–835. pmid:9573407
  28. 28. Loebel A, Tsodyks M. Computation by ensemble synchronization in recurrent networks with synaptic depression. J Comp Neuroscience. 2002; 13: 111–124.
  29. 29. Chizhov AV. Conductance-Based Refractory Density Approach: Comparison with Experimental Data and Generalization to Lognormal Distribution of Input Current. Biol. Cybernetics. 2017; 111(5–6): 353–364. pmid:28819690
  30. 30. Bazhenov M. Potassium Model for Slow (2–3 Hz) In Vivo Neocortical Paroxysmal Oscillations. J Neurophysiology. 2004; 92: 1116–1132.
  31. 31. Chizhov A, Amakhin D, Zaitsev A. Computational model of interictal discharges triggered by interneurons. PLoS ONE. 2017; 12(10): e0185752. pmid:28977038
  32. 32. Krishnan GP, Bazhenov M. Ionic dynamics mediate spontaneous termination of seizures and postictal depression state. J Neuroscience. 2011; 31(24): 8870–8882.
  33. 33. Avoli M, Barbarosie M, Lucke A, Nagao T, Lopantsev V, Kohlig R. Synchronous GABA-mediated potentials and epileptiform discharges in the rat limbic system in vitro. J Neurosci. 1996; 16: 3912–3924. pmid:8656285
  34. 34. Steuer R. Effects of stochasticity in models of the cell cycle: from quantized cycle times to noise-induced oscillations. J Theor Biology. 2004; 228: 293–301.
  35. 35. Lindner B. Coherence and stochastic resonance in nonlinear dynamical systems. Ph.D. Thesis, Humboldt-Universitaet zu Berlin. 2001.
  36. 36. Bernardo M, Budd C, Champneys AR, Kowalczyk P. Piecewise-smooth Dynamical Systems: Theory and Applications (Vol. 163). London: Springer; 2008.
  37. 37. Heinemann U, Lux HD, Gutnick MJ. Extracellular free calcium and potassium during paroxsmal activity in the cerebral cortex of the cat. Exp Brain Res. 1977; 27: 237–243. pmid:880984
  38. 38. Boucetta S, Chauvette S, Bazhenov M, Timofeev I. Focal generation of paroxysmal fast runs during electrographic seizures. Epilepsia. 2008; 49: 1925–1940. pmid:18616553
  39. 39. Stafstrom CE. Distinct mechanisms mediate interictal and pre-ictal discharges in human temporal lobe epilepsy. Epilepsy Curr. 2011; 11(6): 200–202. pmid:22130144
  40. 40. Huberfeld G, de la Prida LM, Pallud J, Cohen I, Le Van Quyen M, Adam C. et al. Glutamatergic pre-ictal discharges emerge at the transition to seizure in human epilepsy. Nature Neuroscience. 2011; 14: 627–634. pmid:21460834
  41. 41. Fertziger AP, Ranck JB Jr. Potassium accumulation in interstitial space during epileptiform seizures. Exp Neurol. 1970; 26(3): 571–585. pmid:5435740
  42. 42. Frohlich F, Bazhenov M, Iragui-Madoz V, Sejnowski TJ. Potassium Dynamics in the Epileptic Cortex: New Insights on an Old Topic. Neuroscientist. 2007; 14: 422–433.
  43. 43. Buchin A, Chizhov A, Huberfeld G, Miles R, Gutkin BS. Reduced Efficacy of the KCC2 Cotransporter Promotes Epileptic Oscillations in a Subiculum Network Model. J Neurosci. 2016; 36(46): 11619–11633. pmid:27852771
  44. 44. Proix T, Bartolomei F, Chauvel P, Bernard C, Jirsa VK. Permittivity coupling across brain regions determines seizure recruitment in partial epilepsy. J Neurosci. 2014; 34(45): 15009–15021. pmid:25378166
  45. 45. Ziburkus J, Cressman JR, Barreto E, Schiff SJ. Interneuron and pyramidal cell interplay during in vitro seizure-like events. J Neurophysiol. 2006;95(6): 3948–3954. pmid:16554499
  46. 46. Ziburkus J, Cressman JR, Schiff SJ. Seizures as imbalanced up states: excitatory and inhibitory conductances during seizure-like events. J Neurophysiol. 2013; 109(5): 1296–1306. pmid:23221405
  47. 47. Kim CM, Nykamp DQ. The influence of depolarization block on seizure-like activity in networks of excitatory and inhibitory neurons. J Comput Neurosci. 2017; 43(1): 65–79. pmid:28528529
  48. 48. Farisello P, Boido D, Nieus T, Medrihan L, Cesca F, Valtorta F, et al. Synaptic and extrasynaptic origin of the excitation/inhibition imbalance in the hippocampus of synapsin I/II/III knockout mice. Cereb Cortex. 2013; 23(3): 581–593. pmid:22368083
  49. 49. Zaitsev AV, Lewis DA. Functional properties and short term dynamics of unidirectional and reciprocal synaptic connections between layer 2/3 pyramidal cells and fast spiking interneurons in juvenile rat prefrontal cortex.European J Neuroscience. 2013; 38(7): 2988–2998.
  50. 50. Izhikevich EM. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting. 2007, The MIT Press Cambridge.
  51. 51. Kramer M. Pathological pattern formation and cortical propagation of epileptic seizures. J. R. Soc. Interface. 2005; 2: 113–127. pmid:16849171
  52. 52. Meijer HGE, Eissa TL, Kiewiet B, Neuman JF, Schevon CA, Emerson RG et al. Modeling Focal Epileptic Activity in the Wilson–Cowan Model with Depolarization Block. J Math Neuroscience. 2015; 5:7.
  53. 53. Wendling F, Benquet P, Bartolomei F, Jirsa V. Computational models of epileptiform activity. J Neurosci Methods. 2016; 260: 233–251. pmid:25843066
  54. 54. Steyn-Ross M., Steyn-Ross D., Sleigh J., Whiting D. 2003; Theoretical predictions for spatial covariance of the electroencephalographic signal during the anesthetic-induced phase transition: increased correlation length and emergence of spatial self-organization. Phys. Rev. E 68, 021902
  55. 55. Kramer MA, Lopour BA, Kirsch HE, Szeri AJ. Bifurcation control of a seizing human cortex. Phys Rev E. 2006; 73(1): 041928.
  56. 56. Perkins KL. Cell-attached voltage-clamp and current-clamp recording and stimulation techniques in brain slices. J Neurosci Methods. 2006; 154(1–2): 1–18. pmid:16554092
  57. 57. Neher E. Correction for liquid junction potentials in patch clamp experiments. Methods Enzymol. 1992; 207: 123–131. pmid:1528115