Reductionist modeling of calcium-dependent dynamics in recurrent neural networks

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.

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:

$\begin{array}{ll}{C}_m\frac{dv}{dt}=-\left(g_l(v-E_L)+I_{AHP}\right)+I_{syn}+I_{app}& (1)\end{array}$

where C_m is the membrane capacitance, which has a value of 1μF, g_l is the leakage conductance with g_l = 0.18μS, E_L is the leak reversal potential with E_L = −60mV, and I_syn denotes the sum of synaptic currents. In this study, I_AHP is Ca²⁺ dependent K⁺ current (Huguet et al., 2016; Sloper and Powell, 1979; Bazhenov et al., 2001a,b; Terman et al., 2002). The applied current I_app is set to a value of 0.2nS.

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

The intrinsic I_AHP current is described by the following:

$\begin{array}{ll}{I}_{AHP}=g_{AHP}(v-E_K)\frac{\left[Ca\right]}{\left[Ca\right]+k_1}& (2)\end{array}$

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

$\begin{array}{ll}{\left[Ca\right]}^{\prime }=-k_{ca}\left[Ca\right]& (3)\end{array}$

with the decay rate parameter $k_{Ca}=0.001m{s}^{-1}$.

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:

$\begin{array}{ll}{C}_m\frac{dv}{dt}=-g_l(v-E_L)+I_{stim}& (4)\end{array}$

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

2.1 Synaptic currents

The synaptic inhibition between the inhibitory cells is given as I_syn = I_I→I+I_E→I. The currents I_I→I and I_E→I are modeled as

$I_{I\to I}=g_i{s}_i(v-E_{GABA})$

$I_{E\to I}=g_e{s}_e(v-E_{AMPA})$

with peak synaptic conductance g_i = 25μS and g_e = 4μS, respectively. The synaptic activation variables s_i and s_e are reset to 1 with every action potential and decay exponentially according to the equations

$s_i^{\prime }=-{\beta }_i{s}_i$

and

$s_e^{\prime }=-{\beta }_e{s}_e$

with decay rates given as ${\beta }_i=0.1m{s}^{-1}$ and ${\beta }_e=2m{s}^{-1}$, 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.

Figure 1

Figure 1. 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.

Figure 2

Figure 2. 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 (I_AHP), which gradually decreases the excitability of the aIC (see Figure 3, second panel).

Figure 3

Figure 3. 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 I_AHP and I_I→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 I_AHP and I_I→I currents (Figure 4). In Figures 4A, B, different initial conditions for Ca²⁺ concentrations and varying values of the Ca²⁺ decay parameter (k_Ca) 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²⁺ levels and other key model parameters.

Figure 4

Figure 4. 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 $C{a}_1^{2+}(0)=0mM$ and $C{a}_2^{2+}(0)=3mM$ vs. $C{a}_1^{2+}(0)=0mM$ and $C{a}_2^{2+}(0)=4mM$, both with a calcium decay rate of $k_{Ca}=0.001{\mathrm{\text{ms}}}^{-1}$. (B) In the second simulation, k_Ca is varied: $k_{Ca}=0.001{\mathrm{\text{ms}}}^{-1}$ vs. $k_{Ca}=0.008{\mathrm{\text{ms}}}^{-1}$ with initial calcium values of $C{a}_1^{2+}(0)=0mM$ and $C{a}_2^{2+}(0)=3mM$ for IC1 and IC2, respectively.

3.2 Approximation of pre-excitation potential values

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

$\frac{d\left[Ca\right]}{dt}\cong o(ϵ)\text{and }\frac{d{s}_i}{dt}\cong o(ϵ).$

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

$-(I_l+I_{AHP}+I_{I\to I})+I_{app}=0$

or

$-(g_l(v-E_L)+g_{AHP}{x}_{Ca}(v-E_K)+g_i{s}_i(v-E_{GABA}))+I_{app}=0$

where $x_{Ca}:=\frac{\left[Ca\right]}{\left[Ca\right]+k_1}$ (see Equation 2). Solving for v, we obtain v as a function of [Ca]

$\begin{array}{ll}{v}_{sub}(\left[Ca\right]):=E_K+\frac{g_l(E_L-E_K)+g_i{s}_i(E_{GABA}-E_K)+I_{app}}{g_l+g_{AHP}{x}_{Ca}+g_i{s}_i}& (5)\end{array}$

If we can estimate [Ca] and s_i 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 s_i decays exponentially from its initial value of 1, with decay rate β_i. Thus, s_i can be approximated explicitly as $s_i=e^{-{\beta }_i{t}_3}$, where t₃ 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, v_tresh. 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:

$\frac{d[Ca]}{dt}=\{\begin{array}{ll}{I}_C,\hfill & \text{if }v=v_{tresh}\hfill \\ -k_{Ca}[Ca]\hfill & \text{otherwise }\hfill \end{array}$

Here, the constant I_C 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 a₀ and s₀ 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, v_tresh. 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 a₀+I_C. Thus, at the conclusion of an inter-spike interval, we obtain the following:

${\left[Ca\right]}_a=(a_0+I_C)e^{-k_{Ca}{t}_3}$

Here, t₃ 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:

${\left[Ca\right]}_s=s_0{e}^{-k_{Ca}{t}_3}$

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 $r=e^{-k_{Ca}{t}_3}$. Thus, at the conclusion of the first inter-spike interval, we have:

${\left[Ca\right]}_a(1)=a_{0}r+I_{C}r$

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

${\left[Ca\right]}_a(2)=(a_{0}r+I_{C}r)r+I_{C}r=a_0{r}^2+I_C(r^2+r)$

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

${\left[Ca\right]}_a(k)=a_0{r}^k+I_C(r^k+r^{k-1}+...+r)=a_0{r}^k+I_{C}r\frac{1-r^k}{1-r}.$

By replacing the expression $\frac{I_C}{1-r}$ with $\overline{A}$, we derive the following:

$\begin{array}{ll}{\left[Ca\right]}_a(k)=r^k{a}_0+\overline{A}(1-r^k)& (6)\end{array}$

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

$\begin{array}{ll}{\left[Ca\right]}_s(1)=s_{0}r,{\left[Ca\right]}_s(2)=s_0{r}^2,...,{\left[Ca\right]}_s(n)=s_0{r}^n.& (7)\end{array}$

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 s_i values alone is not enough to estimate the NSPB; we should, in addition, quantify the effects of [Ca] and s_i 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 s_i 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).

Figure 5

Figure 5. 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 ([Ca²⁺]) 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 I_AHP 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 x₀ and y₀ 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

$\begin{array}{ll}{n}_{xy}(x_0,y_0):=ceil(n_{Ca}(x_0,y_0))& (8)\end{array}$

with

$\begin{array}{ll}{n}_{Ca}(x_0,y_0)=\frac{ln\left(\frac{-m_2+\sqrt{m_2^2-4m_1{m}_3}}{2m_1}\right)}{lnr}& (9)\end{array}$

where for x ∈ ℝ ceil(x): = min{z ∈ ℤ:z > x}. Here m₁, m₂, m₃ are functions of x₀, y₀ and system parameters.

Proof. Let v_a(n): = v_sub([Ca]_a(n)) and v_s(n): = v_sub([Ca]_s(n)). Also define x₀: = [Ca]_a(0) and y₀: = [Ca]_s(0). Assuming the interchange between the active and silent inhibitory cells will occur right after when v_a(n) = v_s(n) and using the Equation 5, we obtain

$\begin{array}{ll}\frac{g_l(E_L-E_K)+I_{app}}{g_l+g_{AHP}{x}_n}=\frac{g_l(E_L-E_K)+g_i{s}_i(E_{GABA}-E_K)+I_{app}}{g_l+g_{AHP}{y}_n+g_i{s}_i}& (10)\end{array}$

with

$\begin{array}{ll}{x}_n:=\frac{{\left[Ca\right]}_a(n)}{{\left[Ca\right]}_a(n)+k_1}& (11)\end{array}$

$\begin{array}{ll}{y}_n:=\frac{{\left[Ca\right]}_s(n)}{{\left[Ca\right]}_s(n)+k_1}& (12)\end{array}$

(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:

$\begin{array}{ll}\frac{g_l+g_{AHP}{x}_n}{g_l(E_L-E_K)+I_{app}}=\frac{g_l+g_{AHP}{y}_n+g_i{s}_i}{g_l(E_L-E_K)+g_i{s}_i(E_{GABA}-E_K)+I_{app}}.& (13)\end{array}$

Let us introduce the following parameters.

$\begin{array}{l}a=\frac{g_l}{g_l(E_L-E_K)+I_{app}}\\ b=\frac{g_{AHP}}{g_l(E_L-E_K)+I_{app}}\\ c=\frac{g_l+g_i{s}_i}{g_l(E_L-E_K)+g_i{s}_i(E_{GABA}-E_K)+I_{app}}\\ d=\frac{g_{AHP}}{g_l(E_L-E_K)+g_i{s}_i(E_{GABA}-E_K)+I_{app}}.\end{array}$

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

$\begin{array}{l}a+b{x}_n=c+d{y}_n\end{array}$

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

$\begin{array}{ll}a+b\frac{{\left[Ca\right]}_a(n)}{{\left[Ca\right]}_a(n)+k_1}=c+d\frac{{\left[Ca\right]}_s(n)}{{\left[Ca\right]}_s(n)+k_1}& (14)\end{array}$

The Equation 14 can be rewritten as follows:

$\begin{array}{l}a+b\left(1-\frac{k_1}{{\left[Ca\right]}_a(n)+k_1}\right)=c+d\left(1-\frac{k_1}{{\left[Ca\right]}_s(n)+k_1}\right)\end{array}$

Rearranging this equation, we obtain the following equation:

$\begin{array}{ll}\frac{-b}{{\left[Ca\right]}_a(n)+k_1}+\frac{d}{{\left[Ca\right]}_s(n)+k_1}=\frac{-a-b+c+d}{k_1}& (15)\end{array}$

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 $m:=\frac{-a-b+c+d}{k_1}$, we obtain the following equation:

$\begin{array}{l}\frac{-b}{r^n{a}_0+\overline{A}(1-r^n)+k_1}+\frac{d}{s_0{r}^n+k_1}=m\end{array}$

This equation can be rewritten as follows:

$\begin{array}{l}\frac{-b}{(a_0-\overline{A})r^n+\overline{A}+k_1}+\frac{d}{s_0{r}^n+k_1}=m\end{array}$

Combining the fractions, we obtain the following equation:

$\begin{array}{ll}\frac{\left[d(a_0-\overline{A})-b{s}_0\right]r^n-b{k}_1+d(\overline{A}+k_1)}{s_0(a_0-\overline{A})r^{2n}+\left[s_0(\overline{A}+k_1)+k_1(a_0-\overline{A})\right]r^n+k_1(\overline{A}+k_1)}=m& (16)\end{array}$

Rearranging the Equation 16, we obtain the following equation:

$\begin{array}{l}m{s}_0(a_0-\overline{A})r^{2n}+\\ \left(m\left[s_0(\overline{A}+k_1)+k_1(a_0-\overline{A})\right]-\left[d(a_0-\overline{A})-b{s}_0\right]\right)r^n\\ +(\overline{A}+k_1)\left[m{k}_1-d\right]+b{k}_1=0& (17)\end{array}$

Setting

$\begin{array}{l}{m}_1=m{s}_0(a_0-\overline{A})\\ m_2=\left(m\left[s_0(\overline{A}+k_1)+k_1(a_0-\overline{A})\right]-\left[d(a_0-\overline{A})-b{s}_0\right]\right)\\ m_3=(\overline{A}+k_1)\left[m{k}_1-d\right]+b{k}_1\end{array}$

the Equation 17 becomes the following

$\begin{array}{l}{m}_1{r}^{2n}+m_2{r}^n+m_3=0\end{array}$

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

$\begin{array}{l}{r}^n=\frac{-m_2+\sqrt{m_2^2-4m_1{m}_3}}{2m_1}\end{array}$

As a result, we obtain the following equation:

$\begin{array}{l}n=\frac{ln\left(\frac{-m_2+\sqrt{m_2^2-4m_1{m}_3}}{2m_1}\right)}{lnr}\end{array}$

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 g_AHP and g_i, for both the continuous model and the discrete formula. NSPB is inversely related to the conductance of the AHP current, g_AHP, because as g_AHP increases, less calcium is needed to hyperpolarize the active IC. On the other hand, NSPB increases with the inhibitory conductance g_i, as higher inhibitory conductance requires more calcium to hyperpolarize the active IC, resulting in more spikes per burst.

Figure 6

Figure 6. 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 (g_AHP) and inhibitory synaptic current conductance (g_i) 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.

Figure 7

Figure 7. 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 n_th 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 C_x(k) and C_y(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

$n_{xy}(k):=n_{xy}(C_x(k),C_y(k))$

as the burst number corresponding to [Ca] values C_x(k) and C_y(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:

$\begin{array}{l}{C}_x(k+1)=r^{n_{xy}(k)}{C}_y(k)\\ C_y(k+1)=r^{n_{xy}(k)}{C}_x(k)+\overline{A}(1-r^{n_{xy}(k)}).\end{array}$

Let F:ℝ² → ℝ² be defined as $F(\vec{w})=r^{n_{xy}(\vec{w})}D+d(\vec{w})$ with, $D:=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$ and $d(\vec{w}):=\left[\begin{array}{c}0\\ \overline{A}(1-r^{n_{xy}(\vec{w})})\end{array}\right]$. Then the discrete map takes the following form:

$\begin{array}{ll}\vec{c}(k+1)=F(\vec{c}(k))\text{with }\vec{c}(k):=\left[\begin{array}{c}{C}_x(k)\\ C_y(k)\end{array}\right].& (18)\end{array}$

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

Lemma. Given n ∈ ℕ⁺, let $F_n:{ℝ}^2\to {ℝ}^2$ be defined as $F_n(\vec{x})=r^{n}D\vec{x}+{\vec{d}}_n$ with ${\vec{d}}_n=\left[\begin{array}{c}0\\ \overline{A}(1-r^n)\end{array}\right]$. Then the map

$\begin{array}{ll}\vec{x}(k+1)=F_n(\vec{x}(k))& (19)\end{array}$

has stable fixed points $D_n:={(I_2-r^{n}D)}^{-1}{\vec{d}}_n$.

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

Proof. Since n_xy(D_k) = k, by lemma, D_k is a fixed point of the map (Equation 18). Again by the continuity of the function n_Ca, there exist a δ neighborhood N_k(δ) of D_k, such that k−1 < n_Ca(N_k(δ)) < k. ∀ E_k ∈ N_k(δ), $||F(D_k)-F(E_k){||}_2=r^k||D(D_k-E_k){||}_2\le r^k||D_k-E_k{||}_2<r^k\delta$ where ||x||₂ denotes Euclidean norm of x ∈ ℝ². This shows that D_k 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. n_Ca(x₀, y₀) given in Equation 9 is a continuous function of the parameter g_AHP. For g_AHP ≫ 1, n_Ca(D_k) < 1 and as g_AHP → 0, n_Ca(D_k) → ∞; that is, neurons either interchange in the initial cycle or never interchange, respectively. Then, by continuity of n_Ca, there exists a g_AHP such that k−1 < n(D_k) < k, and hence, n_xy(D_k) = k. By the theorem, D_k 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 g_AHP and g_i. 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.

Figure 8

Figure 8. Bifurcation diagrams illustrating the dependence of stable bursting solutions (NSPB) on parameters g_AHP (left) and g_i (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 I₁, I₂, ..., I_m+1. In response to the EC excitation, I₁, the IC with the lowest calcium level, spikes and inhibits all other ICs. As calcium builds up in I₁, I₂ 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 C_x(k) and C_yi(k), with i = 1, 2, ..., m denote calcium levels of aIC and sICs right after the k_th interchange, respectively. Set n_xy(k) = n_xy(C_x(k), C_y1(k)), then calcium amounts, and corresponding burst numbers can be followed with the following map:

$\begin{array}{ll}{C}_x(k+1)& =r^{n_{xy}(k+1)}{C}_{y1}(k)\\ C_{y1}(k+1)& =r^{n_{xy}(k+1)}{C}_{y2}(k)\\ .........\\ C_{y(m-1)}(k+1)& =r^{n_{xy}(k+1)}{C}_{ym}(k)\\ C_{ym}(k+1)& =r^{n_{xy}(k+1)}{C}_x(k)+\overline{A}(1-r^{n_{xy}(k+1)}).\end{array}$

Define

${\overrightarrow{Ca}}_m(k):={\left[\begin{array}{l}{C}_x(k)\hfill \\ C_{y1}(k)\hfill \\ \mathrm{...}\hfill \\ C_{ym}(k)\hfill \end{array}\right]}_{(m+1)\times 1},{\overrightarrow{d}}_m(k):={\left[\begin{array}{l}0\hfill \\ \mathrm{...}\hfill \\ 0\hfill \\ \overline{A}(1-r^{n_{xy}(k+1)})\hfill \end{array}\right]}_{(m+1)\times 1}\text{and}$

$\begin{array}{l}{E}_m:={\left[\begin{array}{cccccc}0& 1& 0& 0& ...& 0\\ 0& 0& 1& 0& ...& 0\\ .....\\ 0& 0& 0& ...& 0& 1\\ 1& 0& 0& ...& 0& 0\end{array}\right]}_{(m+1)\times (m+1)}.\end{array}$

Then we obtain following non-linear map,

$\begin{array}{ll}{\overrightarrow{Ca}}_m(k+1)=r^{n_{xy}(k+1)}{E}_m{\overrightarrow{Ca}}_m(k)+{\vec{d}}_m(k).& (20)\end{array}$

3.9 s-periodic Solutions

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

$\begin{array}{l}{\vec{d}}_l^m:={\left[\begin{array}{c}0\\ ..\\ 0\\ \overline{A}(1-r^l)\end{array}\right]}_{(m+1)\times 1}\text{and}\\ F_l^m(\vec{x}):=r^l{E}_m\vec{x}+{\vec{d}}_l^m,\vec{x}\in {ℝ}^{(m+1)}.\end{array}$

For a given ${\vec{n}}_s:={(n_1,n_2,...,n_s)}_{1\times n_s}$, we consider the following function of $\vec{x}$

$\begin{array}{r}{F}_{{\vec{n}}_s}^m(\vec{x}):=F_{n_s}^m(F_{n_s-1}^m(...F_{n_1}^m(\vec{x})))\\ =r^{n_1+n_2+...+n_s}{E}_m^s\vec{x}+r^{n_2+...+n_s}{E}_m^{s-1}{\vec{d}}_{n_1}^m+\\ r^{n3+...+n_s}{E}_m^{s-2}{\vec{d}}_{n_2}^m+...+r^{n_s}{E}_m{\vec{d}}_{n_{s-1}}+{\vec{d}}_{n_s}.\end{array}$

If we let ${\vec{d}}_{n_{s}m}:=r^{n_2+...+n_s}{E}_m^{s-1}{\vec{d}}_{n_1}^m+r^{n3+...+n_s}{E}_m^{s-2}{\vec{d}}_{n_2}^m+...+r^{n_s}{E}_m{\vec{d}}_{n_{s-1}}^m+{\vec{d}}_{n_s}^m$

$F_{{\vec{n}}_s}^m$ takes the following form:

$\begin{array}{ll}{F}_{{\vec{n}}_s}^m(\vec{x})=r^{n_1+n_2+...+n_s}{E}_m^s\vec{x}+{\vec{d}}_{n_{s}m}.& (21)\end{array}$

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

Lemma. Let ${\vec{n}}_s,m\mathrm{\text{and}}s$ be given as above. Then the map ${\vec{x}}_{k+1}=F_{{\vec{n}}_s}^m({\vec{x}}_k)$ has stable fixed point $D_{{\vec{n}}_s}$ given as follows:

$\begin{array}{ll}{D}_{{\vec{n}}_s}:={(I_{m+1}-r^{n_1+..+n_s}{E}_m^s)}^{-1}{\vec{d}}_{n_s}.& (22)\end{array}$

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

Theorem. Let ${\vec{n}}_s=(n_1,..,n_s)\mathrm{\text{and}}{D}_{{\vec{n}}_s}=(v_1,..,v_{m+1})$ be defined as above for m, s ∈ ℤ⁺ with s < m. If $n_i-1<n_{xy}(r^{n_1+...+n_{i-1}}(v_i,v_{i-1}))<n_i\mathrm{\text{for}}i=1,2,...\text{s}$, then $D_{{\vec{n}}_s}$ 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 (g_AHP, g_i, k_ca) 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 (g_AHP, g_i, k_ca, 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.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MZ: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. TD: Data curation, Funding acquisition, Validation, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Gen AI was used in the creation of this manuscript.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Bazhenov, M., Stopfer, M., Rabinovich, M., Abarbanel, H. D., Sejnowski, T. J., and Laurent, G. (2001a). Model of cellular and network mechanisms for odor-evoked temporal patterning in the locust antennal lobe. Neuron 30, 569–581. doi: 10.1016/S0896-6273(01)00286-0

PubMed Abstract | Crossref Full Text | Google Scholar

Bazhenov, M., Stopfer, M., Rabinovich, M., Huerta, R., Abarbanel, H. D., Sejnowski, T. J., et al. (2001b). Model of transient oscillatory synchronization in the locust antennal lobe. Neuron 30, 553–567. doi: 10.1016/S0896-6273(01)00284-7

PubMed Abstract | Crossref Full Text | Google Scholar

Behr, J., Empson, R., Schmitz, D., Gloveli, T., and Heinemann, U. (1996). Electrophysiological properties of rat subicular neurons in vitro. Neurosci. Lett. 220, 41–44. doi: 10.1016/S0304-3940(96)13242-0

PubMed Abstract | Crossref Full Text | Google Scholar

Bhattacharjee, A., and Kaczmarek, L. K. (2005). For k+ channels, na+ is the new ca2+. Trends Neurosci. 28, 422–428. doi: 10.1016/j.tins.2005.06.003

PubMed Abstract | Crossref Full Text | Google Scholar

Channell Jr, P., Cymbalyuk, G., and Shilnikov, A. (2007). Applications of the poincare mapping technique to analysis of neuronal dynamics. Neurocomputing 70, 2107–2111. doi: 10.1016/j.neucom.2006.10.091

Crossref Full Text | Google Scholar

Ermentrout, B., and Terman, D. H. (2010). Foundations of Mathematical Neuroscience. Princeton NJ: CiteSeerX.

Google Scholar

Gold, M. S., Shuster, M. J., and Levine, J. D. (1996). Role of a Ca2+-dependent slow afterhyperpolarization in prostaglandin E2-induced sensitization of cultured rat sensory neurons. Neurosci. Lett. 205, 161–164. doi: 10.1016/0304-3940(96)12401-0

PubMed Abstract | Crossref Full Text | Google Scholar

Guckenheimer, J., Gueron, S., and Harris-Warrick, R. M. (1993). Mapping the dynamics of a bursting neuron. Philos. Trans. R. Soc. Lond. B Biol. Sci. 341, 345–359. doi: 10.1098/rstb.1993.0121

PubMed Abstract | Crossref Full Text | Google Scholar

Huguet, G., Joglekar, A., Messi, L. M., Buckalew, R., Wong, S., Terman, D., et al. (2016). Neuroprotective role of gap junctions in a neuron astrocyte network model. Biophys. J. 111, 452–462. doi: 10.1016/j.bpj.2016.05.051

PubMed Abstract | Crossref Full Text | Google Scholar

Ibarz, B., Casado, J. M., and Sanjuán, M. A. (2011). Map-based models in neuronal dynamics. Phys. Rep. 501, 1–74. doi: 10.1016/j.physrep.2010.12.003

Crossref Full Text | Google Scholar

Lee, K., Duan, W., Sneyd, J., and Herbison, A. E. (2010). Two slow calcium-activated afterhyperpolarization currents control burst firing dynamics in gonadotropin-releasing hormone neurons. J. Neurosci. 30:6214–6224. doi: 10.1523/JNEUROSCI.6156-09.2010

PubMed Abstract | Crossref Full Text | Google Scholar

Liu, Z., Tatar, S., Ulusoy, S., and Zeki, M. (2015). Structural stability for the morris-lecar neuron model. Appl. Math. Comput. 270, 261–268. doi: 10.1016/j.amc.2015.08.031

Crossref Full Text | Google Scholar

Lüthi, A., and McCormick, D. A. (1998). H-current: properties of a neuronal and network pacemaker. Neuron 21:9–12. doi: 10.1016/S0896-6273(00)80509-7

PubMed Abstract | Crossref Full Text | Google Scholar

Matthews, E. A., Linardakis, J. M., and Disterhoft, J. F. (2009). The fast and slow afterhyperpolarizations are differentially modulated in hippocampal neurons by aging and learning. J. Neurosci. 29, 4750–4755. doi: 10.1523/JNEUROSCI.0384-09.2009

PubMed Abstract | Crossref Full Text | Google Scholar

Moustafa, A. A., Cohen, M. X., Sherman, S. J., and Frank, M. J. (2008). A role for dopamine in temporal decision making and reward maximization in parkinsonism. J. Neurosci. 28, 12294–12304. doi: 10.1523/JNEUROSCI.3116-08.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Perez-Orive, J., Bazhenov, M., and Laurent, G. (2004). Intrinsic and circuit properties favor coincidence detection for decoding oscillatory input. J. Neurosci. 24, 6037–6047. doi: 10.1523/JNEUROSCI.1084-04.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Rubin, J., and Terman, D. (2000). Analysis of clustered firing patterns in synaptically coupled networks of oscillators. J. Mathemat. Biol. 41, 513–545. doi: 10.1007/s002850000065

PubMed Abstract | Crossref Full Text | Google Scholar

Rubin, J. E., and Terman, D. (2012). Explicit maps to predict activation order in multiphase rhythms of a coupled cell network. J. Mathem. Neurosci. 2, 1–28. doi: 10.1186/2190-8567-2-4

PubMed Abstract | Crossref Full Text | Google Scholar

Sloper, J., and Powell, T. P. S. (1979). A study of the axon initial segment and proximal axon of neurons in the primate motor and somatic sensory cortices. Philos. Trans. R. Soc. Lond. B Biol. Sci. 285, 173–197. doi: 10.1098/rstb.1979.0004

PubMed Abstract | Crossref Full Text | Google Scholar

Terman, D., Ahn, S., Wang, X., and Just, W. (2008). Reducing neuronal networks to discrete dynamics. Phys. D: Nonlinear Phenomena 237, 324–338. doi: 10.1016/j.physd.2007.09.011

PubMed Abstract | Crossref Full Text | Google Scholar

Terman, D., Kopell, N., and Bose, A. (1998). Dynamics of two mutually coupled slow inhibitory neurons. Phys. D: Nonlinear Phenomena 117, 241–275. doi: 10.1016/S0167-2789(97)00312-6

Crossref Full Text | Google Scholar

Terman, D., Rubin, J. E., Yew, A., and Wilson, C. (2002). Activity patterns in a model for the subthalamopallidal network of the basal ganglia. J. Neurosci. 22, 2963–2976. doi: 10.1523/JNEUROSCI.22-07-02963.2002

PubMed Abstract | Crossref Full Text | Google Scholar

Traub, R. D., Contreras, D., Cunningham, M. O., Murray, H., LeBeau, F. E., Roopun, A., et al. (2005). Single-column thalamocortical network model exhibiting gamma oscillations, sleep spindles, and epileptogenic bursts. J. Neurophysiol. 93, 2194–2232. doi: 10.1152/jn.00983.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Wallén, P., Robertson, B., Cangiano, L., Löw, P., Bhattacharjee, A., Kaczmarek, L. K., et al. (2007). Sodium-dependent potassium channels of a slack-like subtype contribute to the slow afterhyperpolarization in lamprey spinal neurons. J. Physiol. 585, 75–90. doi: 10.1113/jphysiol.2007.138156

PubMed Abstract | Crossref Full Text | Google Scholar

Wang, X.-J., Rinzel, J., and Rogawski, M. A. (1991). A model of the t-type calcium current and the low-threshold spike in thalamic neurons. J. Neurophysiol. 66, 839–850. doi: 10.1152/jn.1991.66.3.839

PubMed Abstract | Crossref Full Text | Google Scholar

Wilson, R. I., and Laurent, G. (2005). Role of gabaergic inhibition in shaping odor-evoked spatiotemporal patterns in the drosophila antennal lobe. J. Neurosci. 25, 9069–9079. doi: 10.1523/JNEUROSCI.2070-05.2005

PubMed Abstract | Crossref Full Text | Google Scholar

Zeki, M., and Balcı, F. (2019). A simplified model of communication between time cells: Accounting for the linearly increasing timing imprecision. Front. Comput. Neurosci. 12:111. doi: 10.3389/fncom.2018.00111

PubMed Abstract | Crossref Full Text | Google Scholar

Zeki, M., and Balcı, F. (2020). A simple three layer excitatory-inhibitory neuronal network for temporal decision-making. Behav. Brain Res. 383:112459. doi: 10.1016/j.bbr.2019.112459

PubMed Abstract | Crossref Full Text | Google Scholar

Zeki, M., and Moustafa, A. A. (2017). Persistent irregular activity is a result of rebound and coincident detection mechanisms: a computational study. Neural Netw. 90, 72–82. doi: 10.1016/j.neunet.2017.03.007

PubMed Abstract | Crossref Full Text | Google Scholar