Introduction

The biophysical properties of neurons are well-established as critical determinants of both single-cell computation and network behavior. While complex neuron models allow a detailed analysis at the single-cell level, investigation of network dynamics often requires simpler models that can retain essential features of neuronal computation yet ensure computational efficiency [1–3]. This challenge is further intensified when biophysical single-cell properties change dynamically, such as through alterations in the local environment (e.g., temperature fluctuations, changing ionic concentrations) or neuromodulatory effects. To better understand how such modifications influence network dynamics, it is essential to also incorporate these physiological changes into simplified neuron models. Here, we derive phenomenological neuron models of the integrate-and-fire type that quantitatively capture the dynamics of more complex conductance-based models not only for one specific set of the model’s parameters but also when a physiologically relevant parameter is allowed to vary.

In contrast to models fitted to experimentally measured voltage traces at essentially one specific point in biophysical parameter space, the simplified model captures the two-dimensional bi-furcation structure spanned by the input current as well as an additional physiologically relevant parameter of interest. Moreover, our specific interest with the approach taken is to adequately represent the neuronal response not only in terms of the firing rate but – because of its computational relevance – also the bifurcation that underlies the onset of regular firing at threshold, including switches in bifurcation type when the physiological parameter varies. While we chose to develop phenomenological models capturing the influence of extracellular potassium in this study, the approach is equally amenable to many other biological parameters.

In the processing of information within single cells, the dynamics of action-potential generation play an important role. Although influenced by numerous cellular characteristics, qualitative differences between regular action-potential firing (with all-or-none pulses at threshold) can only arise from three distinct mathematical mechanisms. Their classification coincides but is not entirely identical with the excitability classes 1 and 2 described by Hodgkin (1948) [4] and observed earlier by Arvanitaki (1938) [5], which are based on the frequency of firing at spike onset (with arbitrarily low frequencies for class I and an immediate jump to higher frequencies for class II). These two classes are associated with distinct mathematical bifurcations that underlie the onset of regular spiking: class I firing corresponding to a saddle-node-on-invariant circle (SNIC) bifurcation, while class II firing typically arises from a subcritical Andronov-Hopf bifurcation. A third bifurcation, the saddle-homoclinic orbit (HOM) bifurcation, can also occur and combines properties of class I and class II behaviors, depending on the mode of analysis. The type of bifurcation is a critical determinant of the network state. For instance, neurons with SNIC dynamics tend to produce asynchronous activity when weakly coupled within inhibitory networks. In contrast, the transition to HOM dynamics in the same network configuration results in synchronous activity.

The transition from SNIC to HOM is of particular interest, as it is ubiquitous in the underlying mathematical framework (despite a quantitative dependence of the switching point on other parameters) and has the potential to induce dramatic changes in network behavior. One well-established parameter that induces this transition is the concentration of extracellular potassium. It is directly affected by neuronal spiking. While the outflow of K⁺ ions at each spike is typically compensated for by mechanisms such as Na⁺/K⁺-ATPases and glial buffering, these may not suffice to prevent transient changes during intense activity [6]. Such elevations can alter the neuron’s excitability, forming a feedback loop that can create rich dynamics [7]. Moreover, abnormally elevated levels of extracellular potassium are involved in paroxysmal pathological activities that are specific to epilepsy (seizures) and migraine with aura (spreading depolarization) [8].

To include the effect of extracellular potassium concentration in simpler, network-compatible phenomenological models, we start from a conductance-based model including the potassium concentration dynamics and fit a simpler, integrate-and-fire-type model equipped with a potassium variable to the conductance based model, aiming to capture both firing rates and the two-dimensional bifurcation structure with input current and extracellular potassium concentration as parametric dimensions. Fitting the bifurcation structure of the reduced model near the SNL bifurcation accurately captures the conductance-based model behavior under a range of ionic concentrations, associated with qualitative changes in excitability. This permits the study of transitions between those different regimes as potassium evolves dynamically. Accordingly, the phenomenological model captures a transition from class I (SNIC) to HOM firing with increases in extracellular potassium. The approach contrasts the more widespread practice of deriving phenomenological models via a fit of spike trains or voltage traces from electrophysiological recordings, as they have been used in international competitions, see for example [9, 10].

In this paper, we first demonstrate how to quantitatively derive the integrate-and-fire type version of a given class I conductance-based model. We then compare, at the single-neuron level, the dynamics of the two models, and their phase-response curves. Finally, we investigate the effect of potassium on neuronal synchronization at the network level.

Results

Reducing a class I conductance-based model to a QIF neuron with dynamic

We demonstrate our procedure on a Wang-Buzsáki model embedded in an extracellular space with dynamic potassium concentration (Fig. 1a) [7, 11]. We reduce it to a model of type quadratic integrate-and-fire (QIF), which allows us to capture the normal forms of both SNIC and HOM bifurcations, as well as the transition between them at the so-called saddle-node-loop (SNL) point via the reset voltage. Reset values below the voltage at the saddle-node bifurcation of the QIF correspond to SNIC onsets, while larger values correspond to HOM onsets. We carefully choose parameter values of the QIF system and incorporate potassium dependencies to match the homoclinic, saddle-node and SNIC bifurcations branches of the target model near the SNL transition. In particular, potassium concentration feeds back on the reset voltage, inducing the SNL. An additional reset rule, which increments potassium concentration by a small amount at each spike, mimics the flow of ions through voltage-gated potassium channels (Fig. 1b). Including an electrogenic Na⁺/K⁺-ATPase complements the potassium dynamics. The resulting reduced model is, like the original, a slow-fast system where slow changes in potassium drive the neuron across different excitability regimes.

FIG. 1:

Original model

For convenience, we reproduce in “Methods” the definition of the extended Wang-Buzsáki model [7], that we use as example to illustrate the reduction procedure.

Relying on the difference in timescale between and the other variables, the fast subsystem of System (4) can defined as Eqs. (4a) to (4c), where is a parameter. At physiological potassium concentrations spike onset is mediated by a SNIC bifurcation (e.g. Fig. 2A when ), while a larger concentrations it is mediated by a homoclinic bifurcation (Fig. 2B, ) [7]. Fig. 3A shows the SNL bifurcation organizing this transition.

FIG. 2:

FIG. 3:

When the neuron spikes, to isolate the slow components of the dynamics, in our case overall changes in extracellular potassium concentration, the averaging method can be used to define the averaged slow subsystem (Eq. (6)). The strategy consists in averaging, at given concentrations, the derivative of over one limit action potential cycle of the fast subsystem.

Step 1: QIF version of the Wang-Buzsáki model for a constant

We consider the following QIF model:

Its branch of fixed points, defined by I_app − I_sn,0 − I_pump = a(v − v_SN)², folds at (I_app = I_sn,0 + I_pump, v = v_SN). Stable limit cycles emerge via a SNIC bifurcation when the reset voltage v_r is smaller than v_SN (Fig. 3c), while they emerge via a homoclinic bifurcation at when v_r > v_SN (Fig. 3d), with a SNL bifurcation at v_r = v_SN. Let us first determine, for constant, parameter values of System (1) such that its bifurcation diagram with respect to the applied current closely resembles the one of the target model at this concentration.

I_SN,0 and v_SN are simply the coordinates of the saddle-node bifurcation in the target model, in the absence of pump current. I_pump, which in the target model only depends on potassium (Eq. (5d)), can here be treated as a constant negative contribution to the applied current.

We take a and c to be the quadratic normal form coefficients for the saddle-node bifurcation of the target model. They are defined as (see for example [12]):

where p and q are the normalized left and right eigenvectors of the Jacobian matrix and B the Hessian quadratic form of the vector field at the saddle-node:

Since the bias term c is close to 1 µF cm⁻², we set it to this value for simplicity.

Then, if is beyond the value corresponding to the SNL bifurcation in the target model [K⁺]_o,SNL, we choose v_r such that the homoclinic bifurcation in the QIF model takes place at the same value I_hom of applied current as in the target model:

When we have more flexibility in the choice of reset. v_r should be smaller than v_sn to obtain a SNIC bifurcation, and a reasonable matching of the limit cycles’ minimum voltage is appreciable. We choose to define the reset in the SNIC case in the continuity of the homoclinic case (see Fig. 4d). This satisfies the two above criteria.

FIG. 4:

Finally, the threshold voltage v_th allows us to approximate the limit cycles’ maximum voltage.

Step 2: QIF version of the Wang-Buzsáki model with dynamic

We now compute parameter values of System (1) for a range of , and depending on how strongly they vary with potassium, either set them to a constant or replace them with a function of .

Fig. 3a-d shows I_SN,0, v_SN and a for between 4 and 16 mM and v_r only above [K⁺]_o,SNL. We notice that v_SN and the normal form coefficient a are relatively conserved in this interval. We fix them to their value at [K⁺]_o,SNL. For I_SN and v_r we use linear fits, constrained at [K⁺]_o,SNL, that we extrapolate to lower concentrations in the case of v_r. The increase of the reset voltage with reflects the increase of the potassium reversal potential in the conductance-based model.

For the threshold v_th, we use a sigmoid, that we extrapolate to lower concentrations, to fit the peak voltage of the limit cycles in the bifurcation diagram with respect to , when I_app = I_app,SNL and in the absence of pump (Fig. 3e). This enables us to capture the shrinkage in spike amplitude as potassium accumulates.

We incorporate these functions in System (1) and include potassium dynamics to obtain the following model:

with

The evolution of the extracellular potassium concentration depends on the one hand on uptake by the Na⁺/K⁺-ATPase (Eq. (2b)), and on the other hand on increments at each action potentials (Eq. (2d)). When the target model is Wang-Buzsáki, the parameter values of the corresponding QIF version in the form of (2) are listed in Table 1.

TABLE I:

The fast subsystem of this model consists in Eqs. (2a) and (2c), with treated as a parameter. Its bifurcation diagram with respect to input and potassium, and the frequency of the limit cycles, can easily be calculated analytically (see “Methods”). We can see in Fig. 3 that, with our reduction procedure, we achieve a very close match of the bifurcation branches emanating from the SNL bifurcation in the target model. In Wang-Buzsáki, the transition to depolarization block at larger potassium concentrations is mediated by a fold of limit cycles and a subcritical Hopf bifurcation (panel A). The dimensionality of the QIF fast subsystem does not allow for such bifurcations. Instead, we artificially obtain a depolarization block by silencing the dynamics when the threshold voltage v_th drops to the reset level (see gray and green curves in Fig. 4e). Note that, although the firing frequency in the fast subsystem grows to infinity as approaches this point, this does not effectively occur in the full system since is incremented by a constant amount with each spike.

The averaged slow subsystem of System (2) is simply given by

where T is the period of the limit cycles.

Our phenomenological neuron model closely captures the dynamics of the original conductance-based version

Under identical stimulation, the behavior of the phenomenological model resulting from this bifurcation diagram fitting procedure is very similar to that of the target model. In the same way as the original version, it can produce regular tonic firing (Fig. S2), as well as richer firing patterns. We illustrate this point with three scenarios typically occurring in class I conductance-based model [7, 13, 14], that can notably be obtained in the Wang-Buzsáki model (see overview in Fig. 5). They share the common feature of involving switches in firing regimes, caused by an activity-mediated accumulation of extracellular potassium. Those switches rely on: (i) the different regimes associated with each region of their common bifurcation diagram, (ii) slow (on average) changes in . Although the time traces are not exactly the same, we will see that the two models are quantitatively comparable.

FIG. 5:

Scenario 1: fold/homoclinic bursting

We first consider the hysteresis loop bursting mechanism described in [7], which relies on the interplay between slow changes in extracellular potassium concentration and voltage dynamics. It requires only fast-inactivating sodium channels, delayed rectifier potassium channels, and an indispensable electrogenic pump, without the need for any slow ion channels.

Starting from a physiological potassium concentration ( [15]), we stimulate the Wang-Buzsáki neuron with a constant deterministic applied current (I_app = 0.25 µA cm⁻²). This corresponds to an initial condition in the spiking region (Fig. 6a), i.e., where the only attractor is a stable limit cycle with large voltage amplitude (Fig. S1a). There, potassium release through voltage-gated channels outweighs the effect of the pump on average over the spike cycles, leading to a progressive accumulation of extracellular potassium (first part of Fig. 6b).

FIG. 6:

When the system arrives in the area enlarged in Fig. 6a, where the bistable region is bounded from below by the fold bifurcation (lower ) and from above by the homoclinic bifurcation (higher ), it begins to burst (second part of Fig. 6b). For completeness, note that this configuration is only a necessary condition for the bursting behavior to occur (for more details, see [7]). During each spiking phase, the system follows the branch of stable limit cycles (Fig. 6d), in the course of which potassium increases overall (Fig. 6e) because the increments at each action potential are not fully compensated by the pump (Fig. 6c). When it reaches the homoclinic bifurcation, the system switches to the branch of stable equilibria, initiating a quiescent phase. In the absence of spiking, the pump can restore potassium concentration down to the fold bifurcation, at which point the system jumps back to the stable limit cycle branch. This marks the start of a new bursting cycle.

Starting from the same initial voltage and potassium concentration, and stimulated with the same input, the integrate-and-fire version of the model successfully captures this bursting behavior (Fig. 6f-j). The underlying mechanism is the same as in the original model, with fold and homoclinic dynamical bifurcations governing the onset and termination of bursts. The small discrepancy between the two models in the number of spike per burst and interspike interval is due to the divergence in the frequency of limit cycles (Fig. S1d) and slightly different location of the bifurcations in relation to potassium levels.

We should point out here that the intersection between the x-axis and the time derivative of the averaged slow subsystem, visible in Fig. 6e,j, does not correspond to a stable limit cycle of the full system. Although there exists a potassium concentration for which, in the fast subsystem, the contributions of the pump and voltage-gated channels to potassium dynamics cancel each other out on average, in the full system the hysteresis loop forms and bursting takes place. As explained in [16], this results from the fact that the reduction of the dynamics to the averaged slow subsystem is not valid in the close vicinity of the homoclinic bifurcation.

For this first scenario, in accordance with [7], we adopted a configuration where the potassium concentration rises until the neuron starts bursting. However, note that simulations can also be generated in which the neuron remains in the spiking region, for example by implementing a buffering of the extracellular potassium. This can conveniently be done in both original and reduced versions of the model, by introducing a term in the potassium equation. As a result, when the buffering is sufficiently strong, potassium stabilizes within the spiking region at a concentration corresponding to a fixed point of the average slow subsystem which, in this case, equates to a limit cycle of the full system (Fig. S2).

Scenario 2: stochastic bursting

We now examine a switch in firing regime identified in a version of the Traub-Miles model [17] (a class I conductance-based model) accounting for potassium concentration dynamics. By introducing a noise component in the input, Contreras et al. obtained simulations showing that, after a period of tonic firing, the neuron switches to stochastic bursting, a regime where firing is intermittently interrupted [13]. It is possible to reproduce this behavior with the Wang-Buzsáki model.

We stimulate the neuron with colored noise (see Eq. (8) in “Methods”) of mean input I_app = 0.3 µA cm⁻², chosen sufficiently large to not interfere with the fold/homoclinic bursting behavior (Fig. 7a). Initially, like in the previous scenario, the neuron spikes periodically, causing an accumulation of extracellular potassium (first part of Fig. 7b). When the system reaches the bistable region, i.e., when crosses the fold bifurcation of the fast subsystem, the firing pattern switches to stochastic bursting (second part of Fig. 7b). There, the noisy input triggers jumps between the two attractors: the stable fixed point corresponding to a resting state, and the stable limit cycle corresponding to tonic firing (Fig. 7d).

FIG. 7:

As expected from their common bifurcation structure, when confronted with the same stimulus (same realization of the stochastic process) as the target model, the integrate-and-fire version behaves very similarly (Fig. 7e-h).

Compared to the deterministic fold/homoclinic bursting, for which burst initiation and termination occur only at the boundaries of the bistable region, in this scenario the noise can induce switches between attractors within the bistable region (Fig. 7c-d,g-h).

On a side note, we notice a canard-like behavior in the middle of the burst shown in Fig. 7c-d, where for a short duration the system remains near the saddle point (unstable).

Scenario 3: transition to depolarization block

The last scenario we consider is the transition to depolarization block, which results in transient spike inactivation. This phenomenon is a pathological hallmark: it is the signature of cortical spreading depolarization at the single-neuron level, and might also be involved in epileptiform activity [18]. The induction of a depolarization block, in the case where it is caused by an accumulation of extracellular potassium resulting from neuronal firing, has been extensively studied in conductance-based models with dynamic ion concentrations (e.g. [13, 14, 19–22]). We simulate it here with the Wang-Buzsáki model.

To achieve this, we use the same setup as for the stochastic bursting scenario, but without the noise component in the stimulus. As a result, in the bistable region the system remains on the branch of stable limit cycles, allowing the extracellular potassium concentration to build up without interruption. When exceeds the concentration corresponding to a fold of limit cycles of the fast subsystem, the neuron enters depolarization block (Fig. 8a-d). In comparison, in the stochastic bursting scenario, frequent pauses in firing due to the presence of noise make it easier for the regulatory mechanisms—in our model the Na⁺/K⁺-ATPase—to maintain in the bistable region, preventing the pathological block.

FIG. 8:

When stimulated with the same deterministic input, the integrate-and-fire version of the model produces a voltage trace very similar to the depolarization block observed in the Wang-Buzsáki model (Fig. 8f-g). However, unlike the two other scenarios, this behavior does not entirely rely on the same mechanism as in the original model. The path leading to the block is unchanged, but the transition to the block itself cannot occur via the same dynamical bifurcation, due to the design of the reduced model. Instead, we emulate it by placing the system (both the voltage and ) in refractory period as soon as, at larger values, the reset voltage meets the threshold value (Fig. 8e,h). With this constraint, although the dynamics during the block is only sketchily captured, the desired outcome is for the neuron to ceases spiking. The dynamics preceding the block is, for its part, well captured: the QIF voltage trace features both the reduction in spike amplitude and increase in firing frequency (Fig. 8g). These two effects are in agreement with the experimental recordings in [20]. As shown in Fig. S1f,h, the frequency of the limit cycles of the fast subsystem as a function of increases in two stages before the block. This pattern is a remnant of the situation with two homoclinic bifurcations, at lower values of applied current (Fig. S1b,d).

Reshaping of the phase-response curve

We now compare the phase-response curves (PRCs) of the two models. For a spiking neuron, the PRC describes how the phase of the cycle at which a perturbation is received affects the timing of the next action potential. It provides a link between the dynamics of single neurons and behavior at the network scale.

Fig. 9a shows the infinitesimal PRC (iPRC, i.e., for infinitesimal input amplitude) of the Wang-Buzsáki model near spike onset, for different potassium concentrations. Like when the SNL is reached via other biophysical parameters (e.g. see [24] for the capacitance), we can see a transition from almost symmetric at low [K]⁺ to left-tilted PRCs after the SNL bifurcation. Near spike onset, the alteration of the iPRC occurs suddenly at the SNL, as illustrated with the cases and note that for the latter two, the curves are nearly identical.

FIG. 9:

The PRC reshaping at the SNL is well captured by the QIF model (Fig. 9c). In fact, this simplified model offers an intuitive understanding of how an increase in deforms the iPRC. The effects of larger are mainly expressed in the QIF model by a larger reset voltage, which results in skipping the beginning of the spike cycle (Fig. 9d). Since we are dealing here with a one-dimensional system, the iPRC is simply the inverse of the derivative of the spike voltage trace as a function of the phase. Increasing the reset is therefore equivalent to retaining only a scaled version of the iPRC segment corresponding to late phases (Fig. S3). As pointed out in [24], the SNL bifurcation coincides with a reduction of the period at a given distance from spike onset (i.e., sharper f-I curves), resulting in larger iPRC amplitudes (Fig. 9a,c). This period reduction is captured by the QIF version of the model (Fig. 9e).

The changes associated with the SNL bifurcation directly affect the collective behavior of neurons in a network. Based on the theory of weakly coupled oscillators [25], and given a specific network configuration, the iPRC can be used to predict the neurons’ ability to synchronize. For example, in the case of two neurons coupled via weak delta synapses, the coupling function is a scaled version of the iPRC, and its odd part governs how the phase difference between the two neurons evolves over time [23, 26, 27]. When the uncoupled firing frequency differs between the two neurons, the locking range, i.e., the maximal frequency detuning such that there exists a stable fixed point for the dynamics of the phase difference, is twice the amplitude of the odd part of the coupling function [26]. Fig. 9f shows the locking range as a function of potassium, for a voltage perturbation amplitude of 0.15 µV. For both models, the locking range abruptly increases at the SNL, indicating enhanced synchronization ability. In more complex networks, the analysis must be adapted to the particular network architecture (e.g. network topology, synaptic coupling), but the PRC can still serve to predict aspects of the network state [28].

Extracellular potassium influences neurons synchronization in a network

We now investigate the consequences of the changes introduced by the SNL bifurcation at the network scale. In particular, for tonically spiking neurons, we expect that the left-tilted PRCs at larger potassium concentrations induce in-phase synchronization of inhibitory networks [29]. We test this prediction in a basic network configuration. For this, we consider 100 inhibitory neurons, fully coupled with delta synapses. We then impose on all these neurons a common extracellular potassium trace, simulating for instance an external wave reaching the network (Fig. 10a).

FIG. 10:

With both models, the neurons of the network, which were previously out of synchronization (Fig. 10c,f), become in-phase synchronized upon arrival of this wave (Fig. 10d,g).

Disentangling the respective contributions of PRC asymmetry and frequency fluctuations to the effect of the SNL transition on synchronization requires further investigation and lies beyond the scope of this study. Furthermore, we deliberately chose values such that the neurons remain in the tonic spiking region (Fig. 10b,e). Indeed, in the bistable region, where only a fraction of the neurons may be active, the PRC is not a reliable predictor of network behavior. A study of the complex dynamics at work in this region is beyond the scope of this article (see for example [30]).

Discussion

Neuronal dynamics can influence network behaviour, yet the inclusion of biophysical detail in larger network models, both challenges computational efficiency and reduces model generality due to the amount of additional parameters introduced. Commonly, neuron models are fitted to datasets from a single experimental setting and one best-fit parameter set is determined. These models may fail to generalize, if parameters not varied in the experiment are altered in the model. We here derive models that capture the impact of a biophysical parameter across a larger parameter space by fitting not voltage data but bifurcation diagrams. The resulting model captures the dynamical transitions in parameter space and therefore, generalises at the cost of over-fitting at a particular parameter value.

At the scale of a single neuron, we have shown that the phenomenological model accurately reproduces various activity patterns from the original model, without needing to adjust parameter values to each scenario. At the network level, preliminary simulations highlight the important role that extracellular potassium concentration plays in synchronization.

Specifically, we have here shown how to reduce a class I conductance-based model with dynamic potassium concentration to a quadratic integrate-and-fire model. Depending on the concentration of extracellular potassium and its dynamics, single-cell excitability in this model can vary between the SNIC or the HOM type, capturing the transitions in the full conductance-based neuron models. The effect of potassium on the bifurcation structure in the original, conductance-based model is in the reduced QIF model version implemented via a change in the voltage reset variable after a spike. Note that this approach is not limited to potassium; it could be extended to other parameters known to induce an SNL transition (e.g. temperature [26] or capacitance [24]).

A phenomenological neuron model capturing potassium dynamics

Phenomenological models that capture biophysical aspects are particularly valuable for exploring aspects at the network scale through systematic simulations. For example, phenomenological models of spike-frequency adaptation have been used to study network synchronization [2, 31]. In the same vein, our work formulates a simple model of potassium concentration dependence. It provides a minimal description of a dynamics that would typically require a more detailed, higher-dimensional system.

Extracellular potassium concentration is a particularly pertinent variable to retain in simplified neuron models. Originally, models assumed constant ion concentrations. This is motivated by the fact that the amount of relocated ions per action potential may have a negligible effect on the overall concentrations [32, 33], and that concentrations are regulated by mechanisms involving the Na⁺/K⁺-ATPases, cotransporters and glial cells [34, 35]. Nonetheless, in the last decades numerous studies have focused on modeling their dynamics [7, 8, 13, 14, 22, 36–43]. Indeed, experimental studies suggest activity-dependent fluctuations in concentration, in particular transient accumulations of extracellular potassium reaching several millimolars, in case of intense activity or following sensory stimulation [6, 34, 35, 41, 44, 45]. Such changes are expected to be particularly significant locally or in restricted spaces [34, 35]. Moreover, drastic increases in extracellular potassium are known to occur during pathophysiological events such as spreading depolarizations and epileptic seizures, as well as during hypoxia [34]. Both experimental evidence [13, 45–47] and results based on biophysical neuron models (e.g. [7, 8, 13, 14, 22, 36, 37, 43]) indicate that the weakening of the potassium driving force, resulting from activity-mediated extracellular build-up, in turn qualitatively alters neuronal excitability, at both moderate and large elevations. This motivates the need for simple models capturing potassium dynamics, to investigate the consequences of changes in neurons’ intrinsic excitability at the network scale.

There exist other minimalist models with potassium dynamics, mainly in the context of epilepsy [48–50]. The single neuron model by Depannemaecker et al. [48] reproduces a wide range of firing patterns, including seizure-like events, sustained ictal activity and depolarization block, while being relatively simple. It remains nevertheless at a higher level of complexity compared to the present model, with its two-dimensional fast subsystem of type Hodgkin-Huxley and its two slow variables. The “epileptor-2” population model [49] by Chizhov et al. is also four-dimensional, with the advantage of explicitly modeling sodium dynamics. It can successfully simulate interictal and ictal discharges in high potassium conditions. However, the choice of coupling between potassium and voltage dynamics in this model does not permit onset bifurcation switches in individual neurons. In a recent study, Depannemaecker et al. propose an extension of the adaptive exponential integrate-and-fire model (AdEx) that accounts for impaired ionic gradients [50], by introducing a third, slow, variable.

Modeling the critical transitions in parameter space

The question of “what is a good model” has been debated extensively in the epistemological literature and many arguments in favor of simple leaky integrate-and-fire models have been presented [51], for example, if well fitted, they have high predictive power for the voltage state variable. However, if one includes model predictions not only in observable state variables, but also for changed biophysical parameter sets, then only models that share the fundamental bifurcation structure found in the data yield good predictions. Our work, therefore, presents a strategy to derive models that capture the dependency of the system on physiological quantities, by matching the relevant bifurcation diagram. For a related approach in a broader context, see [52]. This way of constructing phenomenological models is a real alternative to approaches based on the fitting of spike trains. It offers advantages that can be particularly beneficial, depending on what the model is meant to achieve.

This is of additional relevance since nerons are not stationary devices, but subject to constant change in their enviromnetal parameters. Sources of variability include temperature [26, 53], ionic concentrations [13, 14] and neuromodulators [54]. These factors may directly depend on the recent history of activation [13, 55]. In some cases, their changes constitute a homeostatic response to other changes, helping to maintain the neurons’ properties in the face of a varying environment [56]. We argue that, when known, fitting the bifurcation diagram with respect to chosen parameters is a judicious way to obtain models tailored to the major source of instability in a given context, as we have done here with the extracellular potassium concentration. Fitting specific spike trains typically ensures model validity only in the idiosyncratic conditions of their construction. In contrast, fitting a bifurcation diagram not only provides better robustness to changing conditions, but also allows for the study of how these changes themselves slowly drive neuronal activity between different regimes.

There already exist canonical phenomenological models that, thanks to their bifurcation landscape, can recapitulate the main features of neural excitability (see for example [57], or [58], which is based on the unfolding of a higher-dimensional bifurcation). However, they are formulated in terms of abstract variables. In some situations, having a model that contains specific bifurcations is not sufficient; their arrangement in terms of biophysical dimensions also matters. For example, the fold/hom bursting mechanism that we studied in the first scenario cannot occur when the unfolding of the SNL bifurcation is aligned with the axis (see [7] for more details). Starting from a canonical model and adapting it to match desired bifurcation branches, as we did here with the QIF model, has the advantage of preserving important quantitative components.

Potassium dynamics and network models

One of the possible applications of our work is to provide a tool for studying the consequences of potassium-induced excitability switches at the network level. By adding potassium diffusion to the network model future studies could explore spatial effects and the interplay between synaptic and diffusive coupling (see [59, 60]).

Possible model extensions

Sodium dynamics

In this work, we restricted ourselves to potassium dynamics, although sodium concentrations are also directly affected by neuronal activity. This choice was made because potassium is the main player in activity-induced SNL transitions and switches in excitability class, while sodium is more involved in spike amplitude reduction [13]. The engagement of the Na⁺/K⁺-ATPase to restore ionic gradients is accounted for in our model by the potassium dependence of the pump current. In [7], it is shown that sodium dynamics is not crucial for the bursting mechanism we studied in the first scenario.

Taking into account the dynamics of the two ion species, for example by tracking extracellular potassium and intracellular sodium concentrations, remains the most precise option [49, 61]. To avoid adding an extra variable, one could rely on electroneutrality to express one concentration as a function of the other [21, 36, 48]. Applying this strategy to our work is an interesting direction for refinement.

Transition to depolarization block

Finally, the way in which depolarization block is obtained in the reduced model (Scenario 3) is not ideal. Instead of being mediated by a supercritical Hopf bifurcation, or a subcritical Hopf coupled to a fold of limit cycles, as in conductance-based models, it was artificially achieved by silencing the dynamics when the reset and threshold voltages meet. Indeed, the fast subsystem of the reduced model, of type quadratic integrate-and-fire, cannot exhibit Hopf bifurcations. Instead, one would need to start from a two-dimensional fast subsystem, like the adEx model. This would not only allow a more accurate modeling of the transition to the block, but also open the possibility of studying the recovery from this pathological state.

Conclusion

In conclusion, we affirm the importance of phenomenological models that capture key underlying biophysical mechanisms. We contribute a model of potassium concentration dependence, readily suitable for integration into large networks. Moreover, our approach, based on the fitting of a bifurcation structure, provides a generalizable framework for deriving similar models for other relevant modulatory variables.

Methods

Extended Wang-Buzsáki [7] model

Model definition

The model obeys the system of differential equations

with currents

and gating functions

The reversal potential for the potassium evolves with according to

where R = 8314.4 mJ K⁻¹ mol⁻¹ is the ideal gas constant, F = 96 484.6 C mol⁻¹ the Faraday’s constant and T = 310 K the temperature. The variables and parameter values are listed in Tables II and III.

TABLE II:

TABLE III:

Slow dynamics

As in [7], for values such that there is a stable limit cycle in the fast subsystem, we define the averaged slow subsystem of System (4) as

where T is the period of limit cycles.

Fast subsystem of the QIF model (Eqs. (2a) and (2c))

Frequency of the limit cycles

Let (tonic spiking region), the period of the limit cycles is given by:

If v_r > v_th and (bistable region):

with

Infinitesimal phase-response curve (iPRC)

Let . If :

If v_r > v_th and :

Noisy input

To introduce noise in the applied current we replace it with:

where η_t is a Ornstein-Uhlenbeck (OU) process defined by

W_t being the Wiener process. The OU process can be interpreted as a smoothly low pass filtered white noise. In our simulations, we used a correlation time τ_η = 20 ms and a standart deviation σ = 0.05 µA cm⁻².

Computational resources

To derive the parameters I_SN,0, v_SN, a and v_vr of the QIF model (1), we performed symbolic computations with the Python library Sympy [63]. In the case of v_sn, we additionally used the numerical solver nsolve.

For numerical continuation, we used the software AUTO-07P [64].

Simulations of the neuron models were carried out with the python package brian2 [65].

Supplementary figures

FIG. S1:

FIG. S2:

FIG. S3:

Acknowledgements

This project has received funding from the European Research Council (ERC) under the Union’s Horizon 2020 research an innovation program (grant agreement no 864243), from the Einstein Foundation Berlin (to SS, EP-2021-621) and from the German Research Council (CRC/TRR 384). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Additional information

Funding

European Research Council

https://doi.org/10.3030/864243

Einstein Foundation Berlin (EP-2021-621)

German Research Council (CRC/TRR 384)