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 incorporate these physiological changes into simplified neuron models as well. 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 the common practice of fitting phenomenological models to spike trains or voltage traces recorded under a single physiological condition [4, 5], we propose to derive phenomenological models by matching a two-dimensional bifurcation structure spanned by the input current as well as an additional physiologically relevant parameter. Here, we demonstrate this approach using the bifurcation diagram of a biologically more complex model (specifically, a conductance-based model including potassium concentration dynamics) as the fitting target. In principle, the same strategy could be applied using a bifurcation diagram reconstructed from experimental data, provided that sufficiently rich measurements are available to characterize the relevant dynamical behavior. Our specific aim is to capture not only the firing rate, but also the bifurcation underlying the onset of regular firing at threshold, including switches in bifurcation type as the physiological parameter varies. While we focus here on the influence of extracellular potassium, the approach is equally applicable to 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 arise from three distinct mechanisms. Their classification coincides but is not entirely identical with the excitability classes I and II described by Hodgkin (1948) [6] and observed earlier by Arvanitaki (1938) [7], 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: the saddle-node-on-invariant circle (SNIC) bifurcation is the canonical class-I onset mechanism, while class II is typically associated with the subcritical Hopf bifurcation. Another bifurcation, the saddle-homoclinic orbit (HOM) bifurcation, shares the arbitrarily low onset frequencies of class I while differing in other properties (see Izhikevich [8] for details). 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 [9].

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 [10, 11]. One well-established parameter that causes 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 [12]. Such elevations alter the neuron’s excitability, forming a feedback loop that can create rich dynamics [13]. 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) [14].

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 this 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 SNIC to HOM firing with increases in extracellular potassium.

Overall, our contribution is twofold: a minimal model of potassium-dependent excitability and, more generally, a bifurcation-guided procedure for embedding physiologically meaningful parameter dependence into canonical reduced neuron models. In the following, 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 [K⁺]_o

We demonstrate our procedure on a Wang-Buzsáki model embedded in an extracellular space with dynamic potassium concentration (Fig. 1a) [13, 15], hereafter referred to as the target model. 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 (see Fig. S1). 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 bifurcation 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 [13] (System (4)), which we use as an example to illustrate the reduction procedure.

Relying on the difference in timescale between [K⁺]_o and the other variables, the fast subsystem of System (4) can be defined as Eqs. (4a) to (4c), where [K⁺]_o is a parameter. At physiological potassium concentrations, spike onset is mediated by a SNIC bifurcation (e.g., when [K⁺]_o = 5 mM, see Fig. 2a), while at larger concentrations it is mediated by a homoclinic bifurcation (e.g., when [K⁺]_o = 12 mM, see Fig. 2b) [13]. Fig. 3a shows the SNL bifurcation organizing this transition.

FIG. 2:

FIG. 3:

When the full system slowly follows a branch of stable fixed points of the fast subsystem’s bifurcation diagram, i.e., when the neuron is quiescent, the dynamics of [K⁺]_o can be studied using the standard slow subsystem (System 6). When the full system slowly follows a branch of stable limit cycles, i.e., when the neuron spikes, each spike causes a rapid increment of [K⁺]_o (note the presence of v in Eq. (4d)). However, these increments are of small amplitude, and the dynamics of [K⁺]_o remains slow on average (see already Fig. 6c). As is customary in such cases [16–18], to isolate the slow component of the dynamics, the averaging method can be applied to define the averaged slow subsystem (System 7). The strategy consists in averaging the derivative of the slow variable [K⁺]_o over one period (i.e., one spike cycle), for each member of the family of fast subsystem stable limit cycles.

Step 1: QIF version of the Wang-Buzsáki model for a constant [K⁺]_o

We consider the following QIF model:

Its branch of fixed points, defined by c(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. 2c), while they emerge via a homoclinic bifurcation at when v_r > v_SN (Fig. 2d), with a SNL bifurcation at v_r = v_SN. Let us first determine, for a given [K⁺]_o, parameter values of System (1) such that its bifurcation diagram with respect to the applied current closely resembles that of the fast subsystem of the target model.

Note that in the target model (System (4)), I_pump depends only on potassium (Eq. (5d)). We can therefore treat this current as a potassium-dependent negative contribution to the applied current. In other words, when [K⁺]_o is fixed, I_pump only horizontally translates bifurcation diagrams with respect to I_app, such as those shown in Fig. 2. The parameter derivation below is performed in the absence of the pump current. I_pump can be added a posteriori, taken to be the same as in the target model (Eq. (5d)).

I_SN,0 is the value of applied current I_app at which the saddle-node bifurcation of interest occurs in the fast subsystem of the target model, in the absence of pump current. v_SN is the membrane potential at the saddle-node. We take a and c to be the quadratic normal form coefficients for the saddle-node bifurcation. They are defined as (see, for example, [19]):

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 cm² µF⁻¹, we set it to this value for simplicity.

Then, if [K⁺]_o is larger than the value at the SNL bifurcation in the fast subsystem of 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 of applied current:

where I_hom,0 is the value of I_app at which the homoclinic bifurcation occurs in the fast subsystem of the target model, in the absence of pump current. Note that I_hom = I_hom,0 + I_pump is the value of I_app at which the homoclinic bifurcation occurs in the presence of a pump current.

When [K⁺]_o < [K⁺]_o,SNL (i.e., in the SNIC case), we have more flexibility in the choice of the reset voltage. v_r should be smaller than v_SN to obtain a SNIC bifurcation. v_r should also be reasonably close to the limit cycles’ minimum voltage in the target model. We choose to define v_r in the SNIC case by extrapolating the linear fit of the values of v_r assigned using Eq. (2) in the homoclinic case (i.e., when [K⁺]_o > [K⁺]_o,SNL) to lower values of [K⁺]_o (see Fig. 4d). This satisfies the two criteria above.

FIG. 4:

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

Step 2: QIF version of the Wang-Buzsáki model with dynamic [K⁺]_o

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

Fig. 4a-d shows I_SN,0, v_SN and a for [K⁺]_o 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,0 and v_r we use linear fits, constrained at [K⁺]_o,SNL, that we extrapolate to lower concentrations in the case of v_r. In conductance-based models, higher [K⁺]_o weakens the potassium concentration gradient, thereby raising E_K. The sodium and potassium reversal potentials typically bound voltage oscillations during spiking (see, for example, [20]), so the minimum voltage of spikes is expected to be higher when [K⁺]_o is larger. This effect is reflected in the QIF model via the increase of the reset voltage v_r (Fig. 4d).

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 [K⁺]_o, when I_app = I_app,SNL and in the absence of pump (Fig. 4e). This enables us to capture the shrinkage in spike amplitude that precedes the fold of limit cycles bifurcation. Note that, following the SNL bifurcation, the peak voltage increases slightly with [K⁺]_o, by less than 1 mV. We do not include this increase in the QIF model.

We incorporate these functions of [K⁺]_o in System (1), in which [K⁺]_o now appears as a parameter. Its bifurcation structure with respect to [K⁺]_o and the applied current I_app closely resembles that of the fast subsystem of the original model (see already Fig. 3). A summary of the reduction procedure up to this stage is provided in the form of an algorithm in Table I. An illustration of this procedure using temperature instead of potassium concentration is shown in Fig. S2. Returning to the potassium example, we then include potassium dynamics to obtain the following model:

with

TABLE I:

The evolution of the extracellular potassium concentration depends, on the one hand, on uptake by the Na⁺/K⁺-ATPase (Eq. (3b), where I_pump is given in Eq. (5d)), and on the other hand, on increments at each action potential (Eq. (3d)). Note that the reset rule Eq. (3d) captures the dependency of potassium dynamics on v in the original model (Eq. (4d)). When the target model is Wang-Buzsáki, the parameter values of the corresponding QIF version in the form of (3) are listed in Table II.

TABLE II:

The fast subsystem of this model consists of Eqs. (3a) and (3c), with [K⁺]_o 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 fast subsystem of 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 [K⁺]_o approaches this point, this does not effectively occur in the full system since [K⁺]_o is incremented by a constant amount with each spike.

The standard slow subsystem of System (3), used to study the dynamics of [K⁺]_o along stable branches of fixed points of the fast subsystem, is given by

As in the original model, along stable branches of limit cycles of the fast subsystem, the dynamics of [K⁺]_o includes a fast component, causing rapid increments at each spike (here, instantaneous jumps via the reset rule Eq. (3d)). The jump size Δ_K is chosen to match the potassium increment per spike in the original model and is therefore small (see already Fig. 6h). Thus, the averaging method can still be employed (see, for example, [8]). The averaged slow subsystem of System (3) is then 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. S4), as well as richer firing patterns. We illustrate this point with three scenarios typically occurring in class I conductance-based models [13, 16, 21], which 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 and (ii) slow (on average) changes in [K⁺]_o. 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 [13], 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 ([K⁺]_o = 5 mM [22]), 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. S3a). 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 (lower [K⁺]_o) by the fold bifurcation and from above (higher [K⁺]_o) by the homoclinic bifurcation, 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 [13]). During each spiking phase, the system follows the branch of stable limit cycles (Fig. 6d), in the course of which potassium increases overall (in the averaged slow subsystem , see purple curve in 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 (in the standard slow subsystem , see the black curve in Fig. 6e), 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 spikes per burst and interspike interval is due to the divergence in the frequency of limit cycles (Fig. S3d) and slightly different location of the bifurcations in relation to potassium levels.

We should point out here that, as explained in Section 9.2.3 of [8], the reduction of the dynamics to the averaged slow subsystem is not valid in the close vicinity of the homoclinic bifurcation. Therefore, the intersection between the x-axis and the time derivative of the averaged slow subsystem (purple curve), visible in Fig. 6e and j, does not imply the existence of a stable limit cycle corresponding to tonic firing in the full system.

For this first scenario, in accordance with [13], 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 averaged slow subsystem, which in this case equates to a limit cycle of the full system (Fig. S4).

Scenario 2: stochastic bursting

We now examine a switch in firing regime identified in a version of the Traub-Miles model [23] (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 [21]. It is possible to reproduce this behavior with the Wang-Buzsáki model.

We stimulate the neuron with colored noise (see Eq. (9) in “Methods”) of mean input I_app = 0.3 µA cm⁻², chosen sufficiently large not to 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 [K⁺]_o 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 observe an interruption of about 80 ms in the middle of the burst highlighted in Fig. 7c, during which the voltage is more depolarized than during interburst intervals. Instead, the system remains near the repelling slow manifold (the saddle branch in Fig. 7d). The transition from following an attracting manifold during the burst to evolving near a repulsive one is a canard-like effect, here induced by noise. Similar canard trajectories can be obtained in a deterministic way, see for example [24].

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 [25]. 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., [16, 21, 26–29]). 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 [K⁺]_o 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 [K⁺]_o 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 in a refractory state as soon as, at larger [K⁺]_o values, the reset voltage meets the threshold value (Fig. 8e and h): the dynamics of both the voltage and [K⁺]_o are stopped, and their values remain fixed thereafter (see Fig. 8f and g). With this constraint, although the dynamics during the block is only sketchily captured, the desired outcome is for the neuron to cease spiking. The dynamics preceding the block, by contrast, is captured satisfactorily: 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 [27]. As shown in Fig. S3f and h, the frequency of the fast subsystem limit cycles as a function of [K⁺]_o 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. S3b and 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 their collective behavior (e.g., synchronization) when connected in a network.

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 [30] for the capacitance), we can see a transition from almost symmetric at low [K⁺]_o to left-tilted PRCs after the SNL bifurcation. Near spike onset, this alteration of the iPRC occurs suddenly at the SNL: the curves for [K⁺]_o = [K⁺]_o,SNL − 4 mM and [K⁺]_o = [K⁺]_o,SNL − 2 mM are still 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 [K⁺]_o deforms the iPRC. 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. Flatter slopes therefore correspond to larger iPRC values. At low [K⁺]_o, the inflection point of the spike is near phase (Fig. S5b), resulting in an almost symmetric iPRC (Fig. S5a). The effects of larger [K⁺]_o are mainly expressed in the QIF model by a larger reset voltage, which results in skipping the beginning of the spike cycle (see Fig. 9d or compare panels b and d in Fig. S5). Increasing the reset is thus equivalent to retaining only a scaled version of the iPRC segment corresponding to late phases (see Fig. 9c or compare panels a and c in Fig. S5), resulting in more left-tilted iPRCs.

As pointed out in [30], the SNL bifurcation coincides with a reduction of the period at a given distance from spike onset (i.e., sharper f-I curves). This period reduction is captured by the QIF version of the model (Fig. 9e).

Note that, at large [K⁺]_o, iPRCs contain a negative region (corresponding to phase delay) early in phase. This feature is not captured in the QIF model, related to the fact that QIF models do not explicitly describe membrane repolarization. Their iPRCs therefore cannot contain negative regions.

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 [32], 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 weakly pulse-coupled neurons (where each spike induces a small, instantaneous voltage jump in the postsynaptic neuron), 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 [9, 31, 33]. 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 [9]. Fig. 9f shows the locking range as a function of extracellular potassium concentration 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 and synaptic coupling), but the PRC can still serve to predict aspects of the network state [34].

Extracellular potassium influences neuronal 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 extracellular potassium concentrations induce in-phase synchronization of inhibitory networks [11]. We test this prediction in a basic network configuration. For this, we consider a network of 100 all-to-all pulse-coupled inhibitory neurons. We then impose on all these neurons a common extracellular potassium trace, simulating, for instance, an external wave reaching the network (Fig. 10a). With both models, the neurons of the network, which were previously out of synchronization (Fig. 10c and f), become in-phase synchronized upon arrival of this wave (Fig. 10d and g).

FIG. 10:

Disentangling the respective contributions of PRC asymmetry and frequency fluctuations to the effect of the SNL transition on synchronization requires further investigation. Furthermore, we deliberately chose [K⁺]_o values such that the neurons remain in the tonic spiking region (Fig. 10b and 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 detailed study of the network dynamics is beyond the scope of this article (see, for example, [35]).

Discussion

Neuronal dynamics can influence network behavior, yet the inclusion of biophysical detail in larger network models both challenges computational efficiency and reduces model generality due to the number 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 when parameters not varied in the experiment are altered in the model. Here, we derive models that capture the impact of a biophysical parameter across a physiologically relevant range by fitting bifurcation diagrams rather than voltage traces tied to specific recording conditions. The resulting models capture the dynamical transitions across the fitted parameter range, at the cost of reduced precision at individual points in parameter space.

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 and HOM types, capturing the regime transitions observed in the conductance-based neuron model. 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 reset voltage. Note that this approach is not limited to potassium; as illustrated for temperature in Fig. S2, it can be extended to other parameters known to induce an SNL transition.

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, 36]. 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. Typically, 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 [37, 38], and that concentrations are regulated by mechanisms involving the Na⁺/K⁺-ATPase, cotransporters, and glial cells [39, 40]. Nonetheless, in the last decades numerous studies have focused on modeling their dynamics [13, 14, 16, 21, 27, 29, 41–47]. 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 [12, 39, 40, 45, 48, 49]. Such changes are expected to be particularly significant locally or in restricted spaces [39, 40]. 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 [39]. Both experimental evidence [21, 49–51] and results based on biophysical neuron models (e.g., [13, 14, 16, 21, 27, 29, 41, 47]) 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 [K⁺]_o 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 [52–54]. The single neuron model by Depannemaecker et al. [52] 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 [53] 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 [54], 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 [55]. 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, preserving the relevant bifurcation structure can improve robustness.

The idea that neuronal models should be understood through their bifurcation structure is not new in computational neuroscience. Beginning with FitzHugh’s phase-plane analysis and reduced models [56, 57], the emphasis gradually shifted from reproducing every ionic current to capturing the essential dynamical mechanisms underlying neuronal behavior. Building on this perspective, Rinzel and Ermentrout demonstrated that classes of neural excitability correspond directly to bifurcation types [58]. Rinzel also introduced slow–fast analysis to explain bursting dynamics [59]. Bertram et al. extended this viewpoint by classifying bursting according to the unfolding of codimension-2 or higher bifurcations of the fast subsystem [60], a framework that was subsequently broadened by Golubitsky et al. [61]. Izhikevich further organized this landscape by classifying bursters according to the pair of codimension-1 bifurcations that bound their silent and active phases [8, 62]. This framework continues to be developed and extended [30, 63, 64], making the construction of reduced models that match the relevant bifurcation diagram a natural next step.

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 [65]. This way of constructing phenomenological models is a real alternative to approaches based on fitting spike trains, in particular because neurons are not stationary devices, but subject to constant change in their environmental parameters. Sources of variability include temperature [9, 66, 67], ionic concentrations [16, 21, 68] and neuromodulators [69]. These factors may directly depend on the recent history of activation [21, 70]. 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 [71]. 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 may result in models whose validity is largely restricted to the conditions of their construction. In contrast, fitting a bifurcation diagram may not only provide better robustness to changing conditions but also allow for the study of how these changes themselves slowly drive neuronal activity between different regimes.

In this study, we use as a target the bifurcation diagram of a conductance-based neuron model. This target is only indirectly constrained by experimental data, through gating functions fitted to voltage-clamp recordings. Ideally, the target bifurcation diagram would instead be obtained directly from experimental measurements. Estimating bifurcation diagrams from experiments, however, remains challenging. Some techniques exist to approximate branches of stable and unstable steady states, and branches of stable limit cycles can be approximated using slowly ramped current-clamp protocols [72]. In practice, however, many factors vary during recordings, and generating reliable diagrams would require tightly controlled experimental repetitions in various conditions, particularly for reconstructing two-dimensional bifurcation diagrams. Developing such protocols is beyond the scope of the present study but represents an important direction for future work.

There exist canonical phenomenological models that, thanks to their bifurcation landscape, can recapitulate the main features of neural excitability (see, for example, [73], or [63], 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 [K⁺]_o axis (see [13] 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.

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 rather involved in spike amplitude reduction [21]. 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 [13], 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 [53, 74]. To avoid adding an extra variable, one could rely on electroneutrality to express one concentration as a function of the other [27, 28, 52]. Applying this strategy to our work is an interesting direction for refinement.

Potassium diffusion and spatial network models

A natural extension of the present model is its use in spatially structured networks in which extracellular potassium can diffuse between local compartments. The reduced model is particularly suited to this purpose because it retains the potassium-dependent transitions in neuronal excitability while remaining computationally inexpensive enough for large network simulations. At the single-neuron level, Fig. S4 already illustrates that adding a simple potassium diffusion term can stabilize extracellular potassium and thereby qualitatively alter the resulting firing regime. In a spatial network, allowing potassium released by active neurons to spread to neighboring extracellular compartments would introduce an additional, non-synaptic interaction between neurons, through which local activity can modify the excitability of nearby cells. Such models could be used to investigate how synaptic and potassium-mediated diffusive coupling interact to shape synchronization, spatial propagation of activity, and transitions between excitability regimes (see [75, 76]).

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 model [13]

Model definition

The model obeys the system of differential equations

with currents

and gating functions

The reversal potential for the potassium evolves with [K⁺]_o 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 III and IV.

TABLE III:

TABLE IV:

Slow dynamics

The standard slow subsystem of System (4) is given by

As in [13], for [K⁺]_o 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 the limit cycle.

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

Frequency of the limit cycles

Let = c(I_app − I_SN,0 − I_pump). If (tonic spiking region), the period of the limit cycles is given by:

If v_r > v_SN and (bistable region):

with

Infinitesimal phase-response curve (iPRC)

Let If :

If v_r > v_SN and :

Noisy input

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

where η_t is an 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 standard deviation σ = 0.05 µA cm⁻².

Computational resources

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

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

Simulations of the neuron models were carried out with the Python package brian2 [80].

Supplementary figures

FIG. S1:

FIG. S2:

FIG. S3:

FIG. S4:

FIG. S5:

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript.

Acknowledgements

This project has received funding from the European Research Council (ERC) under the Union’s Horizon 2020 research and 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

EC | European Research Council (ERC)

https://doi.org/10.3030/864243

• Susanne Schreiber

Einstein Stiftung Berlin (Einstein_Berlin) (EP-2021-621)

• Susanne Schreiber

Deutsche Forschungsgemeinschaft (DFG) (CRC/TRR 384)

• Susanne Schreiber