Reductionist modeling of calcium-dependent dynamics in recurrent neural networks

TL;DR

This paper introduces a discrete map model to study calcium-dependent bursting behavior in inhibitory neural networks, showing it can predict key dynamics observed in real systems.

Contribution

The novel discrete map approach captures calcium-dependent bursting dynamics in inhibitory networks by exploiting time-scale differences in neural currents.

Findings

The discrete map accurately predicts the number of spikes per burst based on initial calcium levels.

The model explains how individual currents influence system behavior and can determine fixed points and stability.

The approach successfully replicates canonical bursting signatures observed in continuous systems.

Abstract

Mathematical analysis of biological neural networks, specifically inhibitory networks with all-to-all connections, is challenging due to their complexity and non-linearity. In examining the dynamics of individual neurons, many fast currents are involved solely in spike generation, while slower currents play a significant role in shaping a neuron's behavior. We propose a discrete map approach to analyze the behavior of inhibitory neurons that exhibit bursting modulated by slow calcium currents, leveraging the time-scale differences among neural currents. This discrete map tracks the number of spikes per burst for individual neurons. We compared the map's predictions for the number of spikes per burst and the long-term system behavior to data obtained from the continuous system. Our findings demonstrate that the discrete map can accurately predict the canonical behavioral signatures of bursting performance observed in the continuous system. Specifically, we show that the proposed map a) accounts for the dependence of the number of spikes per burst on initial calcium levels, b) explains the roles of individual currents in shaping the system's behavior, and c) can be explicitly analyzed to determine fixed points and assess their stability.

Figures

Network geometries. Simulations of networks consisting of 2 inhibitory cells (IC) and one excitatory cell (EC), as well as 3 ICs and 1 EC, are presented. In both configurations, the systems exhibit all-to-all coupling. Regular arrowheads indicate excitatory connections, while oval arrowheads represent inhibitory connections.

Oscillations in inhibitory cells. The response of the inhibitory cell (IC) to square pulses results in spike frequency adaptation. With each spike, calcium current activation increases, thereby activating the calcium-dependent potassium current. This potassium current has a hyperpolarizing effect, causing the membrane potential to become more negative.

IC competition. In response to synaptic input from the excitatory cell (EC), inhibitory cells (ICs) compete to spike (upper panel). The key factor determining which IC fires is the activation of the calcium-dependent potassium current (AHP), which increases with rising calcium levels (second panel). Prolonged activation of the AHP current makes the membrane potential more negative, decreasing the likelihood of further spiking. This process determines which IC remains active and which becomes silent.

Example network behaviors. Competition among inhibitory cells (ICs) in a 2IC-1EC network results in varying numbers of spikes per burst (NSPB) depending on the initial calcium levels of the ICs and other key network parameters. (A) Upper: Comparison of two cases with initial calcium concentrations of Ca12+(0)=0 mM and Ca22+(0)=3mM vs. Ca12+(0)=0 mM and Ca22+(0)=4 mM, both with a calcium decay rate of kCa=0.001 ms-1. (B) In the second simulation, kCa is varied: kCa=0.001 ms-1 vs. kCa=0.008 ms-1 with initial calcium values of Ca12+(0)=0 mM and Ca22+(0)=3 mM for IC1 and IC2, respectively.

Comparison of the explicit formula with the full system. Heat maps illustrate the stable number of spikes per burst (NSPB) as a function of initial calcium concentrations ([Ca2+]) for the continuous system (left) and the discrete map using the explicit formula (right). The color gradient indicates the number of spikes per burst, with warmer colors (red, orange) corresponding to higher NSPB values and cooler colors (blue, cyan) representing lower NSPB values. The axes represent the initial calcium concentrations for inhibitory neurons. These visualizations demonstrate the accuracy of the discrete map in capturing the complex dynamics observed in the continuous model.

Comparison of the explicit formula with the full system. Heat maps illustrate the stable number of spikes per burst (NSPB) as functions of afterhyperpolarization conductance (gAHP) and inhibitory synaptic current conductance (gi) for the continuous system (left) and the discrete map model (right). The color gradient represents the NSPB values, with lighter colors indicating higher NSPB values and darker colors indicating lower values. The figures show how the parameters influence bursting behavior and highlight the strong qualitative correspondence between the continuous model and the discrete map.

Schematic diagram illustrating the step-by-step derivation of the discrete calcium map from the continuous model. The workflow begins with voltage and calcium dynamics in the continuous system, applies timescale separation, extracts the number of spikes per burst (NSPB), and updates calcium levels analytically.

Bifurcation diagrams illustrating the dependence of stable bursting solutions (NSPB) on parameters gAHP (left) and gi (right), comparing results from the continuous system (black curves) and the discrete map model (red curves). Each point represents the long-term stable number of spikes per burst (NSPB). Both systems exhibit discrete saddle-node bifurcations, where stable NSPB solutions appear or disappear as parameter values vary, demonstrating that the discrete map effectively captures the qualitative dynamical behavior of the continuous model.

Equations

Full text

1 Introduction

The primary goal of this study is to characterize the role of calcium-dependent adaptation in bursting neurons observed in many brain regions (e.g., see Lee et al., 2010; Wallén et al., 2007; Matthews et al., 2009; Gold et al., 1996; Behr et al., 1996; Bhattacharjee and Kaczmarek, 2005; Guckenheimer et al., 1993; Bazhenov et al., 2001a). While these models arise in many neuronal systems (e.g., see Traub et al., 2005; Bazhenov et al., 2001a; Moustafa et al., 2008; Zeki and Moustafa, 2017), there has been almost no mathematical analysis of biologically realistic models. This difficulty arises because the models are highly nonlinear, include a large number of parameters, and typically exhibit a complex structure of oscillatory behavior (e.g., see Wilson and Laurent, 2005; Perez-Orive et al., 2004; Bazhenov et al., 2001b; Zeki and Balcı, 2019, 2020).

In our study, we initially examine a network consisting of one excitatory neuron and two inhibitory neurons, where the inhibitory neurons are connected through lateral inhibition and receive synaptic input from the excitatory neuron. As a slow modulatory current, inhibitory cells are assumed to include an intrinsic slow afterhyperpolarization (sAHP) current (e.g., see Lee et al., 2010; Wallén et al., 2007; Matthews et al., 2009; Gold et al., 1996; Behr et al., 1996; Bhattacharjee and Kaczmarek, 2005; Guckenheimer et al., 1993; Bazhenov et al., 2001a). For example, in the case of two inhibitory neurons, one neuron fires and inhibits the other. With ongoing activity, calcium in the active cell increases, which activates the sAHP current. The sAHP current decreases the excitability of the active inhibitory cell. At some point, the second inhibitory cell manages to take over the activity. In this network, we aim to determine the number of spikes each active cell makes per burst [number of spikes per burst (NSPB)] and how it depends on important model parameters.

Our approach for studying a general class of networks is to first reduce the full system of differential equations to a discrete-time dynamical system using a Poincaré map-like approach. This is done systematically so that every parameter in the full model corresponds to some parameter, or combination of parameters, in the discrete model. The discrete model has many advantages over the full continuous model; it is considerably easier to solve numerically, allowing for a more systematic study of how the model's behavior depends on parameters. Moreover, in many cases, we can mathematically analyze the discrete model, which is typically impossible for the full system of equations, except in very special circumstances (see Liu et al., 2015; Terman et al., 1998; Rubin and Terman, 2000 for examples of mathematical analysis at the single neuron level and on networks with relaxation oscillator neurons, respectively).

We reduce the full model to a discrete one by constructing a map, which keeps track of the active and silent (inhibited) cells. That is, if we know which inhibitory cell is active and for how long it has been active, then the discrete map determines which cell fires during the next bursting period. However, it turns out that it is not enough to simply keep track of active and silent inhibitory cells. We must also know what the calcium levels of these cells are. In some sense, the calcium levels can be thought of as a slow variable in the sense of geometric singular perturbation theory (e.g., see Ermentrout and Terman, 2010). If we know which cell is active during a bursting episode and what the initial calcium levels of all the cells, and then the map determines the current number of spikes per burst and which cell is going to be active during the next bursting period.

Transient synchronization phenomena, observed prominently in biological neuronal networks such as the olfactory systems of insects and mammals, play a crucial role in sensory information processing and neural coding. Bazhenov et al. (2001a) and Bazhenov et al. (2001b), through extensive experimental recordings and computational modeling of the locust antennal lobe, explicitly showed that transient synchronization emerges from inhibitory neuronal competition mediated predominantly by calcium-dependent potassium currents. Despite the biological clarity provided by such experimental and computational studies, rigorous analytical frameworks explicitly capturing and predicting these transient synchronization behaviors remain underdeveloped. Our study explicitly addresses this gap by deriving analytically tractable discrete maps of calcium-dependent potassium currents, providing explicit formulas for the number of spikes per burst (NSPB). Such analytical maps allow systematic analysis of stability, bifurcations, and parameter sensitivity, significantly deepening theoretical insights into transient synchronization and directly complementing the foundational work by Bazhenov et al. (2001a) and Bazhenov et al. (2001b).

We begin this study by considering small excitatory-inhibitory networks with special network architectures. By considering small networks, we are able to more easily demonstrate how the discrete map is constructed. Moreover, we perform a detailed mathematical analysis of the types of solutions that these networks exhibit and how they depend on parameters. The analysis leads to concrete formulas for the number of spikes each cell exhibits during each episode. In particular, we obtained an explicit formula for the number of spikes per burst, depending on initial calcium values and system parameters (see Lemma 3.4). Using this formula, we constructed an explicit map and analyzed the existence and stability of its fixed points for various networks (see Theorem 8). The discrete map leads to a clear understanding of the roles each component of the model, including the ionic currents, plays in generating the network behavior.

2 Materials and methods

Inhibitory cells (IC) were modeled as a single compartment with currents that are governed by integrate-and-fire kinetics as follows:

$$[eqn]$$

where Cm is the membrane capacitance, which has a value of 1μF, gl is the leakage conductance with gl = 0.18μS, EL is the leak reversal potential with EL = −60mV, and Isyn denotes the sum of synaptic currents. In this study, IAHP is Ca^2+^ dependent K^+^ current (Huguet et al., 2016; Sloper and Powell, 1979; Bazhenov et al., 2001a,b; Terman et al., 2002). The applied current Iapp is set to a value of 0.2nS.

Cells are assumed to fire an action potential when the membrane potential reaches a threshold level vT = −50mV. Once the threshold level is crossed, the membrane potential is reset to the level vR = −75mV.

The intrinsic IAHP current is described by the following:

$$[eqn]$$

with maximal conductance gAHP = 50μS. IAHP rate constant k1 = 10. The reversal potential for IAHP is EK = −90mV. The calcium concentration [Ca] is assumed to increase by 1 μM with every action potential and decays exponentially according to the equation

$$[eqn]$$

with the decay rate parameter $[eqn]$ .

The calcium dynamics used in this model follow an exponential decay form that reflects the experimentally observed behavior of calcium-dependent potassium currents in hippocampal neurons, as reported by Sloper and Powell (1979). This modeling approach has also been employed in several biologically detailed neuronal models (Traub et al., 2005; Bazhenov et al., 2001a,b). While more complex calcium handling mechanisms exist, including buffering and multi-compartmental effects, the effective decay rate used here serves as a biologically grounded approximation that enables analytical tractability without sacrificing key dynamical features relevant to the bursting behavior and inhibitory competition under study.

An excitatory cell (EC) is modeled as follows:

$$[eqn]$$

The parameters are identical to those of the inhibitory cell. The applied current for the EC is given as Istim = 2nS.

2.1 Synaptic currents

The synaptic inhibition between the inhibitory cells is given as Isyn = II→I+IE→I. The currents II→I and IE→I are modeled as

$$[eqn]$$ $$[eqn]$$

with peak synaptic conductance gi = 25μS and ge = 4μS, respectively. The synaptic activation variables si and se are reset to 1 with every action potential and decay exponentially according to the equations

$$[eqn]$$

and

$$[eqn]$$

with decay rates given as $[eqn]$ and $[eqn]$ , respectively.

2.2 Network geometry

We analyze networks consisting of both excitatory and inhibitory neurons, focusing on configurations with one excitatory cell and two or more inhibitory cells (see Figure 1). In these networks, all neurons are connected through all-to-all coupling, allowing for direct interactions between every neuron. This architecture enables competition between inhibitory cells in response to excitatory input, a key aspect we explore in this paper, though the precise mechanism of this competition will be detailed later.

Network geometries. Simulations of networks consisting of 2 inhibitory cells (IC) and one excitatory cell (EC), as well as 3 ICs and 1 EC, are presented. In both configurations, the systems exhibit all-to-all coupling. Regular arrowheads indicate excitatory connections, while oval arrowheads represent inhibitory connections.

3 Results

3.1 Simple networks

3.1.1 Network behavior

In Figure 2, the behavior of a single IC in response to a step current is displayed together with the [Ca] dynamics. In response to the input current, inhibitory cells fire, but inter-spike intervals gradually increase with time. This is due to the AHP current that activates with the build-up of calcium (Figure 2). AHP current is a hyperpolarizing current (making the membrane potential more negative). With continued spiking, interspike intervals reach a point at which the increase in calcium is compensated by its slow decay. Hence, inter-spike intervals eventually stabilize to a constant level.

Oscillations in inhibitory cells. The response of the inhibitory cell (IC) to square pulses results in spike frequency adaptation. With each spike, calcium current activation increases, thereby activating the calcium-dependent potassium current. This potassium current has a hyperpolarizing effect, causing the membrane potential to become more negative.

In a 2IC-1EC network, each inhibitory cell (IC) receives synaptic inhibition from the other (see network architecture, Figure 1). Inhibitory cells compete to spike in response to excitatory input from the excitatory cell (EC). The competition between ICs is primarily driven by the calcium levels within the ICs and the strength and decay rate of the slow inhibitory currents between them.

The network functions as follows: When excitatory input is received, one inhibitory cell, called the active IC (aIC), fires an action potential and inhibits the other, referred to as the silent IC (sIC), via inhibitory-to-inhibitory (I-to-I) slow inhibition (see Figure 3). The aIC may continue spiking for several cycles, maintaining inhibition of the sIC. However, with each spike, calcium builds up in the aIC, activating the calcium-dependent potassium current (IAHP), which gradually decreases the excitability of the aIC (see Figure 3, second panel).

IC competition. In response to synaptic input from the excitatory cell (EC), inhibitory cells (ICs) compete to spike (upper panel). The key factor determining which IC fires is the activation of the calcium-dependent potassium current (AHP), which increases with rising calcium levels (second panel). Prolonged activation of the AHP current makes the membrane potential more negative, decreasing the likelihood of further spiking. This process determines which IC remains active and which becomes silent.

The sIC eventually becomes active after a few cycles, depending on the kinetics and strength of the IAHP and II→I currents, and takes over the spiking activity. This allows the roles of the ICs to switch. During the subsequent burst, calcium levels in the former aIC decay while they increase in the newly active sIC (see Figure 3). Through this process, ICs alternate in firing bursts in response to the excitatory input.

The number of cycles required for the interchange between the active IC (aIC) and the silent IC (sIC), which corresponds to the number of spikes per burst (NSPB), varies depending on the kinetic and conductance parameters, as well as the initial values of the activation variables for the IAHP and II→I currents (Figure 4). In Figures 4A, B, different initial conditions for Ca^2+^ concentrations and varying values of the Ca^2+^ decay parameter (kCa) lead to different NSPB outcomes. The NSPB does not stabilize immediately and shows variability from burst to burst. One of the primary objectives of this study is to characterize the evolution of NSPB and explore its variation based on initial Ca^2+^ levels and other key model parameters.

Example network behaviors. Competition among inhibitory cells (ICs) in a 2IC-1EC network results in varying numbers of spikes per burst (NSPB) depending on the initial calcium levels of the ICs and other key network parameters. (A) Upper: Comparison of two cases with initial calcium concentrations of Ca12+(0)=0 mM and Ca22+(0)=3mM vs. Ca12+(0)=0 mM and Ca22+(0)=4 mM, both with a calcium decay rate of kCa=0.001 ms-1. (B) In the second simulation, kCa is varied: kCa=0.001 ms-1 vs. kCa=0.008 ms-1 with initial calcium values of Ca12+(0)=0 mM and Ca22+(0)=3 mM for IC1 and IC2, respectively.

3.2 Approximation of pre-excitation potential values

At the subthreshold levels, the only active currents are IAHP, IL, and the synaptic inhibition II→I. The changes in the IAHP and II→I are very slow within the inter-spike interval. That is, in the ICs model

$$[eqn]$$

If we let ϵ → 0, then the membrane potential value can be approximated as a fixed point of the potential equation; that is,

$$[eqn]$$

or

$$[eqn]$$

where $[eqn]$ (see Equation 2). Solving for v, we obtain v as a function of [Ca]

$$[eqn]$$

If we can estimate [Ca] and si amounts of the inhibitory cells at the end of inter-spike intervals, we can approximate preburst subthreshold potential values using the above formula.

In Equation 5, the synaptic inhibition variable si decays exponentially from its initial value of 1, with decay rate β_i. Thus, si_ can be approximated explicitly as $[eqn]$ , where t3 is the duration of the inter-spike interval. This demonstrates explicitly how the duration of inhibition influences NSPB calculations.

3.3 Calcium dynamics

The concentration of calcium, denoted as [Ca], is assumed to undergo an immediate surge when the membrane potential attains a predetermined threshold level, vtresh. Subsequently, it experiences exponential decay throughout the inter-spike intervals (see Equation 3).

Specifically, for an active inhibitory cell, the calcium dynamics can be represented by the following differential equation:

$$[eqn]$$

Here, the constant IC represents the instantaneous activation of the calcium current. Solving this differential equation explicitly allows us to determine the calcium value at the conclusion of the inter-spike interval based on the initial calcium value. Let a0 and s0 represent the calcium values of the spiking (active) and inhibited (silent) inhibitory cells right at the onset of the inter-spike interval, i.e., just before the membrane potential reaches the firing threshold, vtresh. Consequently, the calcium level of the active cell, denoted as [Ca]a, undergoes exponential decay according to the equation derived from the initial value of a0+IC. Thus, at the conclusion of an inter-spike interval, we obtain the following:

$$[eqn]$$

Here, t3 represents the duration of the inter-spike interval.

Similarly, the calcium level of the inhibited (silent) cell, denoted as [Ca]s, experiences exponential decay according to the same equation. Hence, at the conclusion of the inter-spike interval, we have the following:

$$[eqn]$$

.

Assuming that the spiking inhibitory cell maintains its activity for k consecutive spikes, we can compute the calcium levels of both the spiking ([Ca]a) and silent ([Ca]s) inhibitory cells at the conclusion of the k^th^ inter-spike interval as follows: Let $[eqn]$ . Thus, at the conclusion of the first inter-spike interval, we have:

$$[eqn]$$

. This value can serve as the initial value for the second inter-spike interval. Consequently, we obtain the following:

$$[eqn]$$

At the conclusion of the k^th^ inter-spike interval, we obtain the following:

$$[eqn]$$

By replacing the expression $[eqn]$ with $[eqn]$ , we derive the following:

$$[eqn]$$

Similarly, we determine the calcium level of the inhibited cell ([Ca]s) at the conclusion of the s0 and the parameter r as follows:

$$[eqn]$$

3.4 Explicit formula for number of spikes per burst (NSPB) depending on pre-burst calcium values and on the system constants

Knowing subthreshold [Ca] and si values alone is not enough to estimate the NSPB; we should, in addition, quantify the effects of [Ca] and si variables on the excitability of aIC and sIC. In order to reach this goal, we consider the Equation 4 of the potential given in the model. The idea to calculate the NSPB is that the interchange between inhibitory cells will occur right after the subthreshold potential values of the neurons are equal at the end of an inter-spike interval (ISI).

Now, combining the ISI-end [Ca] and si calculation and the quantification of their effects on the ISI-end potential value, we can estimate how many cycles it takes for the excitability of the sIC to be more than that of the aIC, that is, NSPB. If we observe the excitability in terms of the ISI-end potential values, we should calculate how many cycles calcium should accumulate in the aIC so that the ISI-end potential of the aIC would be equal to or less than that of the sIC, as in Figure 3. By combining nth-calcium iteration and potential approximations, we obtained an explicit formula (see Equation 8) for the NSPB of the aIC depending on the initial [Ca] values and on the other system parameters and compared NSPBs of the continuous model and explicit NSPB formula, corresponding to different initial [Ca] values (see Figure 5C).

Comparison of the explicit formula with the full system. Heat maps illustrate the stable number of spikes per burst (NSPB) as a function of initial calcium concentrations ([Ca2+]) for the continuous system (left) and the discrete map using the explicit formula (right). The color gradient indicates the number of spikes per burst, with warmer colors (red, orange) corresponding to higher NSPB values and cooler colors (blue, cyan) representing lower NSPB values. The axes represent the initial calcium concentrations for inhibitory neurons. These visualizations demonstrate the accuracy of the discrete map in capturing the complex dynamics observed in the continuous model.

We define one burst as a combination of cycles until the interchange between aIC and sIC occurs. To understand when exactly this interchange occurs, we look at the subthreshold potential behavior of the ICs. With the accumulation of [Ca], the potential level of aIC is reduced via the IAHP current. The interchange will occur right after the potential level of the aIC is equal to that of the sIC (Figure 3). Since we have explicit formulas for the potential level depending on pre-excitation [Ca] values of the aIC and sIC after n consecutive cycles, we can combine these to obtain an equation for NSPB in terms of other system parameters and ISI-end [Ca] values.

Lemma. Let x0 and y0 denote the burst-initial [Ca] amounts of aIC and sIC, respectively. Then, the number of spikes per burst, that is, NSPB, can be given by

$$[eqn]$$

with

$$[eqn]$$

where for x ∈ ℝ ceil(x): = min{z ∈ ℤ:z > x}. Here m1, m2, m3 are functions of x0, y0 and system parameters.

Proof. Let va(n): = vsub([Ca]a(n)) and vs(n): = vsub([Ca]s(n)). Also define x0: = [Ca]a(0) and y0: = [Ca]s(0). Assuming the interchange between the active and silent inhibitory cells will occur right after when va(n) = vs(n) and using the Equation 5, we obtain

$$[eqn]$$

with

$$[eqn]$$ $$[eqn]$$

(see Equation 2). This equation can be solved for n explicitly to obtain the desired formula.

Inverting both sides of the Equation 10, we obtain the following:

$$[eqn]$$

Let us introduce the following parameters.

$$[eqn]$$

Substituting the above variables in the Equation 13 we obtain the following simple form.

$$[eqn]$$

After substituting xn and yn from the Equations 11, 12 into this equation, we obtain the following equation:

$$[eqn]$$

The Equation 14 can be rewritten as follows:

$$[eqn]$$

Rearranging this equation, we obtain the following equation:

$$[eqn]$$

We previously calculated the formulas for the calcium levels of the spiking and inhibited neurons after n spikes (see Equations 6, 7, respectively). Substituting the formulas for [Ca]a(n) and [Ca]s(n) into the Equation 15 and defining $[eqn]$ , we obtain the following equation:

$$[eqn]$$

This equation can be rewritten as follows:

$$[eqn]$$

Combining the fractions, we obtain the following equation:

$$[eqn]$$

Rearranging the Equation 16, we obtain the following equation:

$$[eqn]$$

Setting

$$[eqn]$$

the Equation 17 becomes the following

$$[eqn]$$

This equation is a quadratic equation for the variable r^n^ and can be solved to have

$$[eqn]$$

As a result, we obtain the following equation:

$$[eqn]$$

3.5 Comparison of explicit “number of spikes per burst (NSPB)” formula with continuous system simulations

In Figure 5, we illustrate how NSPB depends on the initial calcium levels using both the continuous model (left panel) and the discrete formula (right panel) derived in the previous Lemma. Notably, the discrete formula qualitatively captures the relationship between NSPB and the initial [Ca] values of the inhibitory cells.

When the difference in initial calcium levels between the inhibitory cells is small, both the discrete formula and the continuous model yield a low NSPB, as expected. In such cases, only a few spikes are needed to reduce the excitability of the active inhibitory cell compared to the silent one. Conversely, when the gap between initial calcium levels increases, more spikes are required to hyperpolarize the active IC (see Figure 5).

In Figure 6, we present heat maps showing the dependence of NSPB on key model parameters, specifically gAHP and gi, for both the continuous model and the discrete formula. NSPB is inversely related to the conductance of the AHP current, gAHP, because as gAHP increases, less calcium is needed to hyperpolarize the active IC. On the other hand, NSPB increases with the inhibitory conductance gi, as higher inhibitory conductance requires more calcium to hyperpolarize the active IC, resulting in more spikes per burst.

Comparison of the explicit formula with the full system. Heat maps illustrate the stable number of spikes per burst (NSPB) as functions of afterhyperpolarization conductance (gAHP) and inhibitory synaptic current conductance (gi) for the continuous system (left) and the discrete map model (right). The color gradient represents the NSPB values, with lighter colors indicating higher NSPB values and darker colors indicating lower values. The figures show how the parameters influence bursting behavior and highlight the strong qualitative correspondence between the continuous model and the discrete map.

3.6 Discrete map

To clarify the analytical approach used to derive the discrete map from the continuous calcium dynamics model, we provide a detailed schematic diagram (see Figure 7). The diagram illustrates each step involved in constructing the discrete map: starting from the continuous model defined by differential equations for voltage and calcium dynamics, identifying timescale separation to treat calcium as a slow variable, estimating the subthreshold voltage, calculating the number of spikes per burst (NSPB) explicitly, and updating calcium levels accordingly. Finally, this procedure yields an analytically tractable discrete map.

Schematic diagram illustrating the step-by-step derivation of the discrete calcium map from the continuous model. The workflow begins with voltage and calcium dynamics in the continuous system, applies timescale separation, extracts the number of spikes per burst (NSPB), and updates calcium levels analytically.

Understanding the NSPB formula allows for the calculation of the initial [Ca] values for the next burst. This is achieved by using formulas that estimate the calcium concentrations of inhibitory cells following n-spikes (see Equations 6, 7, respectively). In this manner, a map that calculates the initial calcium values of the (n+1)th burst based on the nth burst's initial calcium values and other system parameters is constructed (see Equation 18).

The constructed map is non-linear. However, using certain linearization techniques, the fixed points of the map can be calculated explicitly, and their existence and stability can be established.

To understand the long-term behavior of the network, we examine the following non-linear map by combining the iterated [Ca] Equations 7, 6 with the equation for NSPB (Equation 8). Let Cx(k) and Cy(k) denote the calcium levels of the aIC and sIC, respectively, at the beginning of the kth interchange of the neurons with k = 1, 2, .... Let us also define

$$[eqn]$$

as the burst number corresponding to [Ca] values Cx(k) and Cy(k) (see Equation 8). Therefore, the discrete map to find k+1_th_ pre-burst calcium values and corresponding NSPB can be set up as follows:

$$[eqn]$$

Let F:ℝ^2^ → ℝ^2^ be defined as $[eqn]$ with, $[eqn]$ and $[eqn]$ . Then the discrete map takes the following form:

$$[eqn]$$

Since r < 1 with the given parameter set, the following lemma follows.

Lemma. Given n ∈ ℕ^+^, let $[eqn]$ be defined as $[eqn]$ with $[eqn]$ . Then the map

$$[eqn]$$

has stable fixed points $[eqn]$ .

Theorem. Let Dk be defined in the above lemma for k ∈ ℕ^+^. If k−1 < nCa(Dk) < k, then Dk is a stable fixed point of the map (Equation 18).

Proof. Since nxy(Dk) = k, by lemma, Dk is a fixed point of the map (Equation 18). Again by the continuity of the function nCa, there exist a δ neighborhood Nk(δ) of Dk, such that k−1 < nCa(Nk(δ)) < k. ∀ Ek ∈ Nk(δ), $[eqn]$ where ||x||2 denotes Euclidean norm of x ∈ ℝ^2^. This shows that Dk is stable as r < 1.

Corollary. Given k ∈ ℕ^+^, there exists a g_AHP_ value such that D_k_ is a fixed point of the map given in Equation 18.

Proof. nCa(x0, y0) given in Equation 9 is a continuous function of the parameter gAHP. For gAHP ≫ 1, nCa(Dk) < 1 and as gAHP → 0, nCa(Dk) → ∞; that is, neurons either interchange in the initial cycle or never interchange, respectively. Then, by continuity of nCa, there exists a gAHP such that k−1 < n(Dk) < k, and hence, nxy(Dk) = k. By the theorem, Dk is a stable fixed point of the map (Equation 18).

3.7 Bifurcation analysis

To explicitly characterize how parameter changes affect the long-term bursting behavior, we conducted a bifurcation analysis, illustrated clearly in Figure 8. These bifurcation diagrams explicitly show the dependence of the stable number of spikes per burst (NSPB) on critical model parameters, specifically gAHP and gi. As these parameters vary, we observe discrete saddle-node (fold) bifurcations, where stable NSPB solutions either appear or disappear, clearly reflecting transitions in the system's qualitative dynamical behavior. These bifurcations are accurately captured by the discrete map, further validating its effectiveness in modeling the system's dynamics.

Bifurcation diagrams illustrating the dependence of stable bursting solutions (NSPB) on parameters gAHP (left) and gi (right), comparing results from the continuous system (black curves) and the discrete map model (red curves). Each point represents the long-term stable number of spikes per burst (NSPB). Both systems exhibit discrete saddle-node bifurcations, where stable NSPB solutions appear or disappear as parameter values vary, demonstrating that the discrete map effectively captures the qualitative dynamical behavior of the continuous model.

3.8 Generalization to all-to-all couplings

The analysis of small networks is generalized to larger ones. An illustrative example is the nIC-1EC network with all-to-all coupling. Techniques used to analyze the behavior of the 2IC-1EC network are generalized to understand the dynamics of the nIC-1EC network.

The dynamics of the larger network are very similar to 2EC-1EC network: In response to EC excitation, one of the ICs fires and inhibits the other ICs. The firing IC competes with the IC that has the lowest [Ca] among the inhibited ICs. The duration of this competition is determined similarly to the 2IC-1EC network.

A discrete map is constructed for this system by ordering ICs with respect to their initial [Ca] levels (see Equation 20).

Similar linearization techniques used in the analysis of the 2IC-1EC network are applied to this system to understand the long-term behavior of the system. Here, more interesting solutions, such as alternating NSPBs, were obtained. Through this analysis, explicit formulas for fixed NSPBs and k-periodic (alternating NSPBs) solutions are derived.

We now consider a network of m+1 ICs and one EC. All ICs are coupled to all other ICs via inhibition. We order ICs according to their initial calcium amounts as I1, I2, ..., Im+1. In response to the EC excitation, I1, the IC with the lowest calcium level, spikes and inhibits all other ICs. As calcium builds up in I1, I2 will be the second IC to spike, as it has the highest level of excitability among all other ICs since it had the second lowest calcium level initially. Here we assume that inhibition is strong enough and calcium builds up fast enough so that the calcium level of the aIC will be more than that of any other sIC after completing its burst. To analyze the network behavior in the long run, we consider the following map.

Let Cx(k) and Cyi(k), with i = 1, 2, ..., m denote calcium levels of aIC and sICs right after the kth interchange, respectively. Set nxy(k) = nxy(Cx(k), Cy1(k)), then calcium amounts, and corresponding burst numbers can be followed with the following map:

$$[eqn]$$

Define

$$[eqn]$$ $$[eqn]$$

Then we obtain following non-linear map,

$$[eqn]$$

3.9 s-periodic Solutions

Given s ∈ ℤ^+^ with s<m+1 and l ∈ ℕ set

$$[eqn]$$

For a given $[eqn]$ , we consider the following function of $[eqn]$

$$[eqn]$$

If we let $[eqn]$

$[eqn]$ takes the following form:

$$[eqn]$$

Again, since r < 1, we obtain the following lemma.

Lemma. Let $[eqn]$ be given as above. Then the map $[eqn]$ has stable fixed point $[eqn]$ given as follows:

$$[eqn]$$

The following theorem can be proved by the same approach taken to prove the 2IC-1EC case.

Theorem. Let $[eqn]$ be defined as above for m, s ∈ ℤ^+^ with s < m. If $[eqn]$ , then $[eqn]$ is an s-periodic stable point of the map (Equation 20).

4 Discussion

In this paper, we developed a new, simple discrete model of calcium-modulated bursting inhibitory neurons. The discrete model explains key observations regarding the bursting activity patterns of the cells and is closely related to the singular perturbation approach used to analyze biological neural networks. The model consists of two main parts: 1) a nonlinear function that calculates the NSPB of the active (firing) cell depending on the initial calcium levels, and 2) functions tracking the calcium levels of inhibitory cells (active and silent) depending on the calculated NSPB. The observed bursting activity converged to a constant NSPB, which might change depending on the model parameters. Importantly, the resulting fixed points of the discrete map are stable. Given the nonlinear nature of the resulting map, this is a significant outcome of the study.

The present study focuses on the development and mathematical analysis of a calcium-dependent bursting model derived from biologically grounded principles. While direct validation against calcium imaging or electrophysiological data is beyond the scope of this work, our findings complement previous modeling and experimental studies on transient synchronization in olfactory networks (Bazhenov et al., 2001a,b), where inhibitory competition and calcium-dependent potassium currents play a key role.

Although rebound-like activity is observed in our model, it is important to note that the underlying mechanism does not rely on classical post-inhibitory rebound (PIR) currents, such as T-type calcium (Wang et al., 1991) or h-currents (Lüthi and McCormick, 1998). Rather, the rhythmic switching behavior between inhibitory neurons arises from calcium accumulation during spiking, which activates a calcium-dependent potassium current. This hyperpolarizes the neuron and temporarily suppresses its excitability, thereby allowing its counterpart to become active. Thus, the rebound-like behavior in our system is driven by a calcium-mediated inhibitory competition mechanism, consistent with findings in biologically grounded models of olfactory networks (Bazhenov et al., 2001a,b).

The proposed model exhibited several key properties related to the spiking behavior underlying the bursting dynamics of inhibitory cells. First, the discrete model contained many of the biophysical parameters of the continuous model. Second, the model could replicate the NSPB using the initial calcium values and other important model parameters. Finally, the model enabled explicit mathematical analysis. We demonstrated that the proposed model could explain the alternation between inhibitory cells and the associated NSPB in all-to-all coupled nIC-1EC networks.

Several previous studies have explored calcium-dependent bursting using discrete maps or simulations (e.g., see Channell Jr et al., 2007; Rubin and Terman, 2012; Ibarz et al., 2011). Unlike heuristic or purely numerical methods common in prior work, our key contribution is an explicit closed-form formula (Equation 8) that analytically predicts the number of spikes per burst (NSPB). This explicit formulation enables direct, analytical insights into how calcium levels and key parameters (gAHP, gi, kca) affect burst lengths, facilitating stability analyses and natural extensions to larger networks.

Although we primarily illustrate our findings around a single stable equilibrium, the discrete map derived here indeed supports multiple stable solutions depending on parameters and initial conditions (gAHP, gi, kca, initial calcium concentrations). Our method explicitly allows the exploration of multistability, which we discuss further in Section 4.7 for larger network structures.

We compared the behavior of the model to the data obtained from the continuous system. However, we assumed that the number of excitatory neurons remained constant (one) in the simulated systems. In Terman et al. (2008), the authors developed a discrete map that can function with different numbers of excitatory cells in the network. However, no explicit mathematical analysis for the long-term behavior of the network activity was shown, possibly due to the complications caused by the varying synaptic excitation to the inhibitory cells at each cycle.

In this study, we assumed that the switch between the active inhibitory cell and the silent inhibitory cell occurs when their subthreshold potential levels are equalized. Generalization to random connections can be done by including synaptic excitation to the compared subthreshold potentials. This type of generalization would make sense in the case of instantaneous activation of excitatory synaptic currents (AMPA), as opposed to NMDA-type excitatory synaptic currents that operate on much slower time scales. We obtained promising preliminary results for random network architectures with varying excitatory connections to the inhibitory cells, but we did not include them here to keep the focus of the current study on the bursting duration of the inhibitory cells.

The proposed model not only generated NSPB similar to the continuous system data but also explicitly predicted the long-term NSPB. In summary, our model predicted various aspects of the bursting behavior in an nIC-1EC network with all-to-all couplings in a way that made explicit mathematical analysis of the long-term behavior possible. Our results revealed a high degree of overlap between the discrete map and continuous network results.