Multistability of bursting rhythms in a half-center oscillator and the protective effects of synaptic inhibition

For animals to meet environmental challenges, the activity patterns of specialized oscillatory neural circuits, central pattern generators (CPGs), controlling rhythmic movements like breathing and locomotion, are adjusted by neuromodulation. As a representative example, the leech heartbeat is controlled by a CPG driven by two pairs of mutually inhibitory interneurons, heart interneuron (HN) half-center oscillators (HCO). Experiments and modeling indicate that neuromodulation of HCO navigates this CPG between dysfunctional regimes by employing a co-regulating inverted relation; reducing Na⁺/K⁺ pump current and increasing hyperpolarization-activated (h-) current. Simply reducing pump activity or increasing h-current leads to either seizure-like bursting or an asymmetric bursting dysfunctional regime, respectively. Here, we demonstrate through modeling that, alongside this coregulation path, a new bursting regime emerges. Both regimes fulfill the criteria for functional bursting activity. Although the cycle periods and burst durations of these patterns are roughly the same, the new one exhibits an intra-burst spike frequency that is twice as high as the other. This finding suggests that neuromodulation could introduce additional functional regimes with higher spike frequency, and thus more effective synaptic transmission to motor neurons. We found that this new regime co-exists with the original bursting. The HCO can be switched between them by a short pulse of excitatory or inhibitory conductance. In this domain of coexisting functional patterns, an isolated cell model exhibits only one regime, a severely dysfunctional plateau-containing, seizure-like activity. This aligns with widely reported notion that deficiency of inhibition can cause seizures and other dysfunctional neural activities. We show that along the coregulation path of neuromodulation, the high excitability of the single HNs induced by myomodulin is harnessed by mutually inhibitory synaptic interactions of the HCO into the functional bursting pattern.

1 Introduction

Rhythmic circuits, called central pattern generators (CPGs), that control oscillatory motor behaviors like walking, swimming, and breathing, demonstrate a remarkable ability to produce patterns that are robust and, at the same time, can flexibly adjust to changes of the environment and behavioral goals (Marder and Calabrese, 1996; Calabrese, 1998; Doi and Ramirez, 2008; Marder, 2012; Alonso and Marder, 2020; Stadele and Stein, 2022). The neuromodulation making these adjustments does so by producing changes in membrane and synaptic currents often with seemingly opposing effects in a complex non-intuitive manner (Brezina et al., 2003a,b; Brezina et al., 2005; Harris-Warrick and Johnson, 2010; Ellingson et al., 2021). How the bursting neurons of these rhythmic circuits, which inherently display a wide range of biophysical parameters of ionic currents (Goldman et al., 2001; Prinz et al., 2003, 2004; Ellingson et al., 2021) respond to neuromodulation and deliver functional activity patterns with a wide range of the cycle period and spike frequency is not well understood.

The complexity of this imposition of neuromodulation on an inherently variable background, stems from its non-linear dynamical nature. Neurons and neuronal networks can exhibit different activity regimes, e.g., silence, spiking, or bursting, and their co-existence (Guttman et al., 1980; Hounsgaard et al., 1984; Hounsgaard and Kiehn, 1989; Lechner et al., 1996; Perrier and Hounsgaard, 2000; Williams et al., 2002; Loewenstein et al., 2005; Paydarfar et al., 2006; Newman and Butera, 2010; Forrest et al., 2012) only one of which may be functionally appropriate in a given context. Transitions between these regimes set limits to the working range of possible control by neuromodulators. Moreover, neuromodulation can be associated with multistability of activity regime. Multistability refers to the co-existence of more than one regimes in the phase space of the neuronal system. Neurons and their circuits can show a switch between regimes in response to short transient signals or transition between patterns in a state-dependent fashion, demonstrating hysteresis, which is the hallmark of multistability (Canavier et al., 1994; Lechner et al., 1996; Butera, 1998; Nadim et al., 2008; Newman and Butera, 2010; Malashchenko et al., 2011a,b; Marin et al., 2013; Parker et al., 2018, 2021). As a product of neuromodulation, multistability can be functional and desirable or dysfunctional and dangerous. Multistable neurons may play a role as toggle switches and memory units in information processing and working memory (Turrigiano et al., 1996; Fall et al., 2005; Loewenstein et al., 2005; Nadim et al., 2008; Forrest et al., 2012) and could be useful elements of multifunctional central pattern generators (Parker et al., 2018, 2021), but multistability can also be a pathological feature. For example, normal regimes can coexist with seizure regimes (Hahn and Durand, 2001; Frohlich and Bazhenov, 2006; Ziburkus et al., 2006; Cressman et al., 2009; Ullah et al., 2009).

In our previous work, we reported on how myomodulin in a half-center oscillator (HCO) of the leech heartbeat central pattern generator comodulates the h-current and the Na⁺/K⁺ pump current, coordinately increasing the former and decreasing the latter, to control the period of activity in a wide range while avoiding dysfunctional regimes (Ellingson et al., 2021). As the most dangerous dysfunctional regime, we considered the plateau-containing regimes observed at the low values of the maximal pump activity. Similar regimes have been described in leech neurons; widespread seizure-like plateau oscillations are induced by substituting external Ca²⁺ with blocking ions and these plateaus are terminated by pump activity (Angstadt and Friesen, 1991; Opdyke and Calabrese, 1994). The reduction of Na⁺/K⁺ pump activity has been also implicated in the transition to seizure-like plateau-containing activity in principal neurons of mouse hippocampal slices (Krishnan et al., 2015). During seizures the Na⁺/K⁺ pump appears to play a crucial role in termination and in postictal depression (Frohlich and Bazhenov, 2006, Ziburkus et al., 2006, Cressman et al., 2009; Ullah et al., 2009, Krishnan et al., 2015).

We hypothesized that given the prevalence of plateaus at low values of maximal pump current, I_PumpMax, in our leech heart interneuron (HN) HCO model, the half-center structure may protect the system from succumbing to underlying pathological single cell dynamics. We compare maps of activity regimes of HCO model with a single HN model. Recently, we discovered that a single HN with augmented dynamics by upregulating persistent sodium and pump currents with dynamic clamp can produce a new type of bursting with high spike frequency and high amplitude of the underlying voltage envelope (Erazo-Toscano et al., 2023). Reviewing the waveforms of the HCO bursting activity along the myomodulin coregulation path (Ellingson et al., 2021), we noticed a change in the waveform and spike frequency suggesting emergence of a new bursting regime. Here, we present the discovery of a unique, functional multistable domain in an HCO model under modulated conditions and investigate the parameter space of modulation in a single cell HN model. We found that, in HCO, a high-spike-frequency high-voltage-amplitude bursting coexists with a low-spike-frequency low-voltage-amplitude bursting. The cycle period of these two regimes was the same. We found that the multistability of the two regimes appears in the HCO model and is not found in the decoupled HN model. We compare contributions of specific currents to generation of the two bursting rhythms. We describe several regimes and demonstrate how the HCO structure prevents the breakdown of functional activity under modulated conditions.

2 Materials and methods

2.1 Model and parameter sweeps

We have incorporated the Na⁺/K⁺ pump dynamics into the canonical model of the leech heart interneuron (Kueh et al., 2016; Ellingson et al., 2021; Toscano et al., 2021; Erazo-Toscano et al., 2023). The model consists of either a single neuron (single HN model) or two reciprocally inhibitory identical neurons (HCO model). It replicates the electrical activity of the HNs of the leech heartbeat CPG under a variety of experimental conditions (Opdyke and Calabrese, 1994; Olsen and Calabrese, 1996; Cymbalyuk and Calabrese, 2000; Cymbalyuk et al., 2002; Sorensen et al., 2004; Olypher et al., 2006; Tobin and Calabrese, 2006; Weaver et al., 2010; Kueh et al., 2016; Ellingson et al., 2021; Toscano et al., 2021). The individual HN was modeled as a single isopotential compartment with Hodgkin-Huxley type intrinsic membrane currents. It has 8 voltage-gated currents: fast Na⁺ current, I_NaF; persistent Na⁺ current, I_P; fast and slow low-threshold Ca²⁺ currents, I_CaF and I_CaS; h-current, I_h; delayed rectifier-like K⁺ current, I_K1; a persistent K⁺ current, I_K2; and a fast-transient K⁺ current, I_KA. The HCO model includes inhibitory spike-mediated synaptic current. The HCO model is utilized as presented in Ellingson et al. (2021), while the single HN model uses the same set of equations without synapses. Parameter sweeps for the single cell model are conducted in the same manner as described in Ellingson et al. (2021).

2.2 Modeling the Na⁺/K⁺ pump

The activity of the pump depends on the extracellular K⁺ concentration ( ${\left[K\right]}_o$ ) and the intracellular Na⁺ concentration ( ${\left[Na\right]}_i$ ). We assume that the external and internal K⁺, and external Na⁺ concentrations do not change significantly during normal operation of heart interneurons and that the pump rate is not voltage dependent. We model the pump current following Angstadt and Friesen (1991) and Kueh et al. (2016).

$I_{\mathit{pump}}=\frac{I_{\mathit{PumpMax}}}{1+exp\left(\frac{{\left[Na\right]}_{ih}-{\left[Na\right]}_i}{{\left[Na\right]}_{\mathit{is}}}\right)}$ , where $I_{\mathit{PumpMax}}$ is the maximal pump current, ${\left[Na\right]}_{ih}$ is the internal concentration of Na⁺ at which the steady state rate of the pump is one half of the maximal value; ${\left[Na\right]}_{\mathit{is}}$ determines the responsivity of the pump to changes of ${\left[Na\right]}_i$ : the smaller the value of ${\left[Na\right]}_{\mathit{is}}$ is, the steeper is the pump steady state activation curve and the more responsive is the pump current.

The two sodium currents (I_NaF and I_P), I_h, I_leak and the pump current determine ${\left[Na\right]}_i$ . I_h is a mixed cation current carried by Na⁺ and K⁺. To model contributions of I_h and I_leak to ${\left[Na\right]}_i$ , we split these currents into Na⁺ and K⁺ components. Concerning the model of I_pump, every cycle of the pump moves 3 Na⁺ ions out and 2 K⁺ ions into the neuron. The dynamics of ${\left[Na\right]}_i$ is governed by Equation 1.

$\frac{d{\left[Na\right]}_i}{dt}=\frac{-\left(\left(\begin{array}{l}{\overline{g}}_P{m}_P+{\overline{g}}_{NaF}{m}_{NaF}^3{h}_{NaF}\\ +{\overline{g}}_{hNa}{m}_h^2+g_{\mathit{leakNa}}\end{array}\right)\left(V_m-E_{Na}\right)+3I_{\mathit{pump}}\right)}{\mathrm{v}\mathrm{F}}\text{ (1)}$

where $E_{Na}$ is the reversal potential for sodium, ${\overline{g}}_P$ , ${\overline{g}}_{NaF}$ , and ${\overline{g}}_{hNa}$ are the maximal conductances of I_P and I_Na, and the Na⁺ component of I_h, correspondingly, and $m_P$ , $m_{NaF}$ $h_{NaF}$ , $m_h$ are their gating variables, $g_{\mathit{leakNa}}$ is the conductance of Na⁺ component of I_leak, $V_m$ is the membrane potential, v is the volume of the cell, and F = 96,485 C/mol is the Faraday constant. In Equation 1, we used vF = 410 Cnℓ/mol, which gives us realistic size of the neuron: volume 4.2 pℓ and radius 10.05 μm. ${\left[Na\right]}_i$ determines the Na⁺ reversal potential: $E_{Na}=\frac{RT}{F}ln\left(\frac{{\left[Na\right]}_o}{{\left[Na\right]}_i}\right)$ . ${\left[Na\right]}_o$ is constant 0.115 M.

Simulation along the line of coregulation was accomplished by first setting parameter values ( ${\overline{g}}_h$ , I_PumpMax) = (0,0.548) and initializing with the standard set of initial conditions (Supplementary Table S1). Then, after simulation, ${\overline{g}}_h$ is increased by 0.05 nS, and I_PumpMax is decreased in accordance with the coregulation relation established in Ellingson et al. (2021):

$\begin{array}{ll}{I}_{\mathit{PumpMax}}\left({\overline{g}}_h\right)=0.36+\frac{0.16}{\left({\overline{g}}_h+0.85\right)}& (2)\end{array}$

After each simulation, the values of all state variables were then used as initial conditions in the next simulation. This allowed us to smoothly simulate the action of myomodulin.

2.3 Burst analysis

Model activity was analyzed with in-house scripts written in the MATLAB software (The MathWorks, Inc.) in a similar manner as presented in Ellingson et al. (2021).

Depolarized phases of bursting activity were identified as intervals in which the cell rose above-45 mV and remained above this threshold for at least 0.5 s. A burst was defined as an interval in which the cell depolarized and demonstrated continuous spiking activity. Spikes were detected as peaks in the voltage trace above −30 mV, and timing between spikes was used to distinguish full bursts from depolarized phases containing plateau potentials. The beginning of a spike train was identified as a spike which had no predecessor for at least 0.4 s, and the end of a spike train was identified as a spike which had no successor for 0.4 s. If multiple spike trains were detected within a depolarized phase or if no spikes were identified within 0.4 s of the end of a depolarized phase, this depolarized phase would be tagged as a plateau-containing seizure-like oscillation. Depolarized phases which contained an uninterrupted spike-train were identified as bursts. Depolarized phase duration was computed as the interval of time the cell remained above −45 mV. If the depolarized phase was identified as a burst, its burst duration was computed as the interval of time between the first and last spike in the continuous spike train. These measures were nearly identical in bursts (± one interspike interval). Cycle period was computed as the interval between the beginning of a depolarized phase and the beginning of the next depolarized phase. Intraburst spike frequency was computed as the mean of the inverse of interspike intervals within the burst.

In our parameter sweeps of the single HN and HCO models, we observed a variety of different types of activity (Figure 2), which we generally classified as one type of functional activity (functional bursting), or two types of dysfunctional activities: asymmetric bursting, and a plateau-containing activity. Functional bursting was defined as symmetrical or near symmetrical (depolarized phase of neither cell ≥55% of the cycle period) alternation of spiking activity. During a burst, a cell’s membrane potential rises above −45 mV and spikes continuously until it drops below −45 mV. Functional bursting is further differentiated into high-spike-frequency bursting and low-spike-frequency bursting in this article by the intraburst spike frequency (threshold of 18 Hz) and the minimum voltage (threshold of −70 mV) during the post-burst hyperpolarization. Asymmetric bursting activity was defined similarly, but in the asymmetric case bursts of one cell were noticeably longer than those of the other cell. Cases in which the depolarized phase duration of side of the half-center oscillator ≥55% of the period were considered asymmetric. Plateau-like oscillations were characterized by the appearance of plateaus in the depolarized phase of activity, indicating a failure to spike due to a depolarization block.

Figure 1

Figure 1. Map of single-cell HN model activity regimes demonstrates the presence of dysfunctional regimes in a parameter domain which is functional in the HCO model. Blue, yellow, and red indicate low-spike-frequency bursting [LSFB, (B1)], high-spike-frequency bursting [HSFB, (B2)], and plateau-containing oscillations (B3,B4), respectively. Superimposed in black is the myomodulin coregulation path representing experimental data, and the vertical black line denotes the transition into multistability identified in the HCO model along this path. At the control conditions,   $({\overline{g}}_h\text{,}$  I_PumpMax) = (1.6,0.429), the model displays low-spike-frequency bursting (B1). Reduction in I_PumpMax to  $({\overline{g}}_h\text{,}$ I_PumpMax) = (1.6,0.4) results in high-spike-frequency bursting (B2). There is a small region of plateau-containing oscillations where high-spike-frequency bursting and plateau regimes appear mixed at the border of the plateau oscillations region (B3)  $({\overline{g}}_h\text{,}$ I_PumpMax) = (2.6,0.406). This border nearly coincides with the entry into bistability in the HCO (vertical black line). The single HN model representing experimental condition Myomodulin 10 uM exhibits only plateau oscillations (B4) at  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385) (triangle), whereas HCO supports two functional bursting regimes (Figures 2B, 3).

Using these definitions of burst and plateau events, regimes of HN activity could be labeled (1) functional low frequency, (2) functional high frequency, (3) asymmetric, or (4) plateau-containing. Single cell activity cannot be asymmetric, so this regime is not considered in analysis of the single cell model.

2.4 Spike averaging and phase-wise analysis of currents

A spike averaging algorithm was used to smooth out the spiking contribution to bursts to focus on the average contribution of currents to the underlying oscillatory dynamics of bursting. For each spike, state variables, currents, and fluxes were averaged over its duration in the following manner: Once spike times were identified using previously described methods (Ellingson et al., 2021), spike floors were identified as the minimum value of membrane potential between spikes. These points were then used to compute spike duration and to smooth time series of state variables, currents, and Na⁺ fluxes. Between spike floors, all values for a trace were set to the average value between spike floors. To compute average current and flux contributions over a cycle, smoothed depolarized and hyperpolarized phases were sliced out of a cycle, and traces were integrated using the trapezoidal Riemann sum method (function trapz in MATLAB) and normalized to the duration of the phase for comparison at different cycle periods.

3 Results

We developed a model of a leech heart interneuron half-center oscillator (HCO) and described a neuromodulation-governed mechanism for robust control of cycle period (Ellingson et al., 2021). We also identified specific pairs of model parameter values for $I_{\mathit{PumpMax}}$ and ${\overline{g}}_h$ that represent experimentally recorded activities in HCOs under variation of myomodulin concentration. We captured the effects of myomodulin by a simple coregulation path. This was achieved by curve-fitting the identified parameter pairs to an inverse relationship (Ellingson et al., 2021) expressed by Equation 2.

Based on our modeling analysis, we demonstrated that by following this myomodulin coregulation path, the HCO avoids a dysfunctional asymmetric bursting regime at higher levels of $I_{\mathit{PumpMax}}$ and ${\overline{g}}_h$ as well as a dysfunctional seizure-like plateau-containing regime at lower levels of $I_{\mathit{PumpMax}}$ and in doing so achieves larger ranges of burst duration and bursting cycle period.

3.1 Co-existence of low- and high-spike frequency bursting regimes in the HCO model

A two-parameter map of spike frequency indicated emergence of a high-spike-frequency regime within the boundaries of the functional activity of the HCO model (Figure 2A). It’s typical mean intra-burst spike frequency (i.e., spike frequency) was roughly two times higher compared with the low-frequency regime; 24 Hz versus 12 Hz, respectively. Given the sharp transition from the low-spike-frequency bursting (LSFB) regime to the high frequency bursting (HSFB) regime and occasional appearance of the low frequency regime within the domain of the high frequency regime on the map, we hypothesized that these two regimes may coexist within some domain along the experimentally-substantiated myomodulin coregulation curve. This growth of the spike frequency was qualitatively consistent with the trend recorded experimentally (Tobin and Calabrese, 2005) that spike frequency increases under the application of myomodulin by 10–30%.

Figure 2

Figure 2. A coregulation path traverses a domain of co-existent low-spike-frequency and high-spike-frequency bursting regimes with roughly the same cycle period while navigating between dysfunctional regimes. (A) Map of spike frequency in the functional regime overlaid with the experimentally validated path of coregulation in black; dysfunctional regimes are marked in white; points representing control and Myomodulin 10 uM conditions are marked by a black circle or square, respectively. (B1) Bursting activity under control conditions (green trace) is obtained at  $({\overline{g}}_h\text{,}$ I_PumpMax) = (1.6,0.429). (B2) Two functional bursting regimes co-exist at the parameter coordinates  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385) corresponding to experimental Myomodulin 10 uM conditions: low-spike-frequency bursting (blue traces on the top) and the high-spike-frequency bursting (yellow traces on the bottom). Burst duration (C), burst period (D), and spike frequency (E) from simulations along the path of coregulation. Blue points represent the low-spike-frequency regime, and yellow points represent the high-spike-frequency regime, the shaded gray area indicates the multistable domain of parameters, where either regime can manifest given appropriate initial conditions or targeted perturbations.

At the parameter pair ( ${\overline{g}}_h$ , I_PumpMax) = (1.6, 0.429), which represented activity observed in control experiments, model activity is monostable (Figure 2B1). At the parameter pair ( ${\overline{g}}_h$ , I_PumpMax) = (5.4,0.385), which corresponded to application of 10 uM myomodulin, the HCO can produce either the low-spike-frequency regime or the high-spike-frequency (Figure 2B2) depending on the set of initial conditions. To test for multistability along the curve of coregulation, we performed continuous simulation along the path of coregulation by incrementing ${\overline{g}}_h$ and decrementing I_PumpMax in accordance with the coregulation relation Equation 2 with three methods of initializing model state variables (Methods). We noticed the appearance of the high-spike-frequency regime at ( ${\overline{g}}_h$ , I_PumpMax) = (2.5,0.408). Along with the myomodulin coregulation path toward high myomodulin concentration, the HCO system is capable of producing this high-spike-frequency regime until ( ${\overline{g}}_h$ , I_PumpMax) = (8, 0.378), after which the system becomes monostable on the low-frequency regime once again. When the initial conditions were appropriately chosen within the specified parameter boundaries, the HCO neurons were able to exhibit either the low-frequency regime or the high-frequency regime. This was confirmed by simulation for a duration of at least 1,600 s of model time. In both regimes, we observed a decrease in bursting cycle period and burst duration at very similar values as modulatory coregulation proceeds to higher ${\overline{g}}_h$ and lower I_PumpMax (Figures 2C,D). However, the high-spike-frequency regime decreases in spike frequency from 24 Hz to 20 Hz as the coregulation proceeds, while the low spike frequency regime increases in spike frequency from 12 Hz to 15 Hz throughout the domain of bistability (Figure 2E). Thus, along the curve of coregulation, there is a region of multistability in the HCO model between ( ${\overline{g}}_h$ , I_PumpMax) = (2.5,0.408) and ( ${\overline{g}}_h$ , I_PumpMax) = (8, 0.378) in which the typical bursting, e.g., LSFB, coexists with the HSFB regime.

3.2 Regime switching

While both regimes are consistent with functional activity as measured in experiments with HCOs (Kueh et al., 2016; Wenning et al., 2018; Erazo-Toscano et al., 2023), a major question arises as to how the nervous system can control transitions between regimes in order to utilize the potential benefits of the high-spike frequency bursting regime. To investigate this, we mimicked the effect of post-synaptic currents (PSCs) originating from outside the HCO – perhaps as a descending central command or as an input from elsewhere within the full circuit of the leech heartbeat central pattern generator. We represented excitatory and inhibitory PSCs by pulses of conductance for a current with corresponding reversal potential. We were able to produce switches from LSFB to HSFB and vice versa with either excitatory or inhibitory pulses. We found that the success of the switch depended on when within a burst cycle the pulse was delivered and on the magnitude of the conductance pulse.

To investigate the ability to switch between these co-existing regimes we chose a point on the coregulation path at ( ${\overline{g}}_h$ , I_PumpMax) = (5.4,0.385), representing experimental data with 10 μM Myomodulin applied (Ellingson et al., 2021). Simulation of LSFB and HSFB were accomplished by setting parameter values (${\overline{g}}_h$, I_PumpMax) = (5.4,0.385) and initializing with LSFB and HSFB initial conditions, respectively (Supplementary Table 2). The cycle periods for low-and high-spike frequency bursting were notably close, with the former being 5.769 s and the latter slightly longer at 6.184 s. For LSFB, the burst durations of the HCO neurons V₁ and V₂ were 2.42 s and 3.01 s, respectively, exhibiting an asymmetry in the pattern (Assym = 0.22), evaluated following (Ellingson et al., 2021). Conversely, for HSFB, the burst durations for V₁ and V₂ were both 2.75 s, showing no pattern asymmetry (Assym = 0).

By delivering individual pulses with varying magnitudes at different phases of the burst cycle, we investigated the properties of the pulses that triggered a switch. We applied 30 ms pulses of conductance with either an inhibitory (−62.5 mV) or an excitatory (0 mV) reversal potential to one neuron in our HCO model (Figure 3). Since the HCO model is symmetric, for the sake of clarity, we delivered the pulse to a neuron shown at the top, referred to as V₁. Marking the beginning of a cycle with the first spike in the burst of this neuron, we defined the phase of the pulse delivery as the moment in time within a burst cycle, which we normalized by dividing it by the cycle’s duration, and expressed it as a percentage (Figures 3A1,A2). We methodically varied the conductance of these pulses in a range from 1 nS to 101 nS, increasing the conductance in increments of 1 nS for each subsequent pulse. We also delivered pulses at different phases. For HSFB, we applied pulses at 200 different phases, with each phase separated by a time interval of 0.03092 s. We employed the same time interval for LSFB, which resulted in a total of 189 phases, note that the cycle period for LSFB is shorter compared to HSFB. Sweeping systematically pulse amplitude and the phase we mapped 20,000 and 18,900 points for HSFB and LSFB perturbations, respectively. To test whether the switch occurred, we established a threshold frequency of 18 Hz to distinguish between LSFB and HSFB. Counting the switch cases on the maps (Figures 3B1,B2,E1,E2), we found that it is much easier to switch from LSFB to HSFB than vice versa by either an inhibitory or an excitatory pulse with success in 33.7 and 59.5% cases, respectively. Accordingly, switches from HSFB to LSFB were less common. However, our findings indicate that excitatory and inhibitory pulses still facilitated these switches in 4.1 and 9.9% of cases, respectively. This data implies that using an excitatory pulse is more than twice as effective in inducing a switch compared to an inhibitory pulse.

Figure 3

Figure 3. Under modulation, perturbation with either an excitatory or an inhibitory conductance pulse can trigger a switch between the two coexisting HCO activity regimes. At  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385), the HCO model can be switched from the low-spike-frequency bursting (LSFB, blue) to the high-spike-frequency bursting (HSFB, yellow) and vice versa by a 30 ms inhibitory (B1,B2,C1,C2,D1,D2) or an excitatory (E1,E2,F1,F2,G1,G2) pulse of conductance with reversal potentials −0.0625 V and 0.0 V, respectively. The pulse was applied to one neuron, marked as V₁. The success of the switch depends on the phase (A1,A2), starting with the first spike of a burst of V₁, and the amplitude of the pulse conductance (B1,B2,E1,E2). Examples of switches from LSFB to HSFB (C1) and vice versa (C2) with an inhibitory pulse. F1 and F2 exhibit similar examples of switches triggered by excitatory pulses. Panels D1,D2,G1,G2 exhibit the transitions across the threshold frequency 18 Hz.

3.3 Single-cell HN model under modulation

When surgically or pharmacologically isolated and thus without the synaptic input from the rest of the CPG, HN cells can burst endogenously (Cymbalyuk et al., 2002; Erazo-Toscano et al., 2023). Likewise, our single HN model bursts under control conditions, ( ${\overline{g}}_h$ , I_PumpMax) = (1.6,0.429). Based on the prevalence of seizure-like plateau-containing oscillations at low values of I_PumpMax in the HCO model (Ellingson et al., 2021), we hypothesized that the mutually inhibitory interactions of the HCO may be preventing a depolarization block in at least some of the parameter regions of the functional HCO bursting regime. The high-spike frequency regime seen in the HCO model at lower levels of I_PumpMax could be an intermediate regime between low-spike frequency activity and a depolarization block leading to plateaus under higher levels of modulation. To test this hypothesis, we first performed a two-dimensional sweep of $I_{\mathit{PumpMax}}$ and ${\overline{g}}_h$ on the single-cell HN model, investigating the same modulatory parameter domain as in our HCO simulations.

The most apparent difference between oscillatory activity produced by a decoupled HN neuron versus being incorporated in an HCO was in the control of cycle period. In the HCO system, the gradient of period was most closely associated with changes in ${\overline{g}}_h$ (Ellingson et al., 2021). In the single cell, however, throughout most of the tested parameter domain, changes in ${\overline{g}}_h$ have little effect; $I_{\mathit{PumpMax}}$ plays a much bigger role in regulation of the period (Figure 4A). As a demonstration, consider any arbitrary vertical path through the parameter domain shown in Figure 4A. At high levels of $I_{\mathit{PumpMax}}$ , cycle period is around 8 s, and drops as low as 5 s as $I_{\mathit{PumpMax}}$ is reduced. At middling values of $I_{\mathit{PumpMax}}$ this trend abruptly reverses. Cycle period jumps to nearly 10 s and increases to as much as 13.5 s as $I_{\mathit{PumpMax}}$ is continuously reduced. The parameter ${\overline{g}}_h$ appears to play a role in determining the level of $I_{\mathit{PumpMax}}$ which will induce this trend reversal. At ${\overline{g}}_h$ = 0 nS, this reversal occurs at I_PumpMax = 0.35 nA, but at ${\overline{g}}_h$ = 10 nS, this reversal occurs at I_PumpMax = 0.435 nA. This border also appears in maps of depolarized phase duration (Figure 4B) and spike frequency (Figure 4C) suggesting that the single cell undergoes at least one major regime change under modulation. Using our previous methods for identifying plateau potentials in membrane activity (Ellingson et al., 2021), we found that beyond this border at which the period trend reverses, 100% of membrane potential oscillations contain plateau potentials (Figure 1).

Figure 4

Figure 4. Characteristics of the oscillatory activities of a single-cell HN model under two parameter variation  $({\overline{g}}_h\text{,}$ I_PumpMax). (A) Average cycle period of single cell activity. Warmer colors indicate a faster rhythm. Cycle period is more sensitive to parameter I_PumpMax than to ${\overline{g}}_h$ . (B) Average duration of the depolarized phase of bursting or plateau-containing oscillations. (C) Average spike frequency within the depolarized phase. The boundaries of the HCO functional activity regime are superimposed in white.

3.4 Single cell HN model displays several regimes of activity

Similar to the HCO model, based on the above burst characteristics, we were able to categorize single cell activity into three types of regimes: low-spike-frequency bursting activity, high-spike-frequency bursting activity, and seizure-like plateau-containing activity (Figure 1, blue, yellow, and red, respectively). Between the regions of high-spike-frequency bursting and plateaus, we also identified a small domain of hybrid regimes consisting of high-spike-frequency bursting and plateaus (Figure 1B3). The high-spike-frequency regime is qualitatively very similar to the high-spike-frequency regime identified in the HCO model, in the same way that the low-spike-frequency regime is qualitatively similar to the low-spike-frequency regime in the HCO model. Without synaptic inhibition, the HN neurons are able to recover more quickly to initiate the next burst leading to a shorter interburst interval.

We mapped the experimentally derived HCO domain of functional bursting (Figure 2A) with two white dashed curves, marking borders (Figure 4A), onto the corresponding parameter space of a single HN activities (Figure 4A). Then comparing this map to the established map of the activity regimes (Figure 1A), we noticed that the transition from bursting to plateau-like oscillations in the HN model along the coregulation path occurs at ( ${\overline{g}}_h$ , I_PumpMax) = (2.6,0.406). Mutual synaptic inhibition, at middling to high levels of ${\overline{g}}_h$ is sufficient to prevent this depolarization block, allowing for HCO functional activity in parameter spaces which would not be functional without the synaptic inhibition. In fact, 67% of parameter pairs identified as functional in the HCO (1,532 points out of 2,278, Figure 2A), are dysfunctional in a single cell (Figures 1A, 4A). The point at which the single cell transitions to seizure-like plateau-containing activity coincides with the beginning of the domain of multistability in the HCO map. At the same level of modulation, which transitions the single cell into plateaus, the HCO still maintains functional activity, and also gains coexistence of low and high-spike-frequency regimes. In addition, when I_PumpMax is decreased from control values in the single HN, the cell transitions to the high-spike-frequency bursting regime at low levels of ${\overline{g}}_h$ .

3.5 Description of regimes and ionic mechanisms

In this section, we compare trajectories and contributions of the key ionic currents supporting the bursting regimes of the HCO and single cell (Figures 5–7). To investigate the roles of different currents in the underlying dynamics of large oscillations of the membrane potential which include bursts and plateaus, we applied a spike averaging algorithm to the depolarized phase of model voltage traces and currents. This smooths out the contribution of individual spikes to highlight the underlying oscillations which bring the cell to spiking threshold. As can be seen from Figure 7, a notable difference in the characteristics of the low-spike-frequency regime (blue) and the high-spike-frequency regime (yellow) is that the amplitude of intracellular Na⁺ concentration oscillations, which controls the level of pump current, is larger in the high-spike-frequency one. Remarkably, in the HCO it is nearly twice as high in the high-spike-frequency regimes as in the low spike frequency regimes (Figure 7). This results in a higher pump current throughout the burst as well. In low-spike-frequency regimes, the pump current rises and falls during the burst (Figure 5). In high-spike-frequency and plateau-containing regimes, the pump current rises monotonically during the depolarized phase and remains high until the minimum membrane potential is reached in the hyperpolarized phase. This suggests that pump current plays important role in the dynamics of burst termination in these regimes. The pump current remains high in these regimes at the burst termination adding momentum to the outward currents which hyperpolarize the membrane resulting in a much larger afterhyperpolarization. This in turn pushes the membrane potential far below the synaptic reversal potential (−62.5 mV). This means that at the beginning of the hyperpolarized phase, the synaptic current is depolarizing during some time interval, despite it normally functioning as inhibitory.

Figure 5

Figure 5. Current and Na⁺ flux time series of the two bursting regimes produced by the HCO and single HN models under control and modulated conditions. All traces have been averaged over spikes to provide a look at the oscillatory dynamics. Current and Na⁺ flux under control parameters  $({\overline{g}}_h\text{,}$ I_PumpMax) = (1.6,0.429) for HCO (A) and single HN (B), both of which are low frequency bursting regimes. Current and Na⁺ flux under modulated parameters  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385) for HCO, which is bistable – capable of producing either a low-spike-frequency (C) regime or a high-spike-frequency regime (E) and single HN (D) which is a plateau oscillations regime.

Figure 6

Figure 6. Contributions of currents to bursting phases normalized to phase duration in HCO model regimes. Average current for the monostable control (green) regime at  $({\overline{g}}_h\text{,}$ I_PumpMax) = (1.6, 0.429), and the high (gold) and low (blue) spike-frequency regimes at bistable point  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385) during the depolarized (A) and hyperpolarized phase (B). Average Na⁺ flux for each Na⁺ carrying current for the control (green) regime at  $({\overline{g}}_h\text{,}$ I_PumpMax) = (1.6, 0.429), high (gold) and low (blue) spike-frequency regimes at bistable point  $({\overline{g}}_h\text{,}$ I_PumpMax) = (5.4,0.385) during the depolarized (C) and hyperpolarized phase (D).

Figure 7

Figure 7. Comparison of characteristics of the HCO and single HN model oscillatory activities under two parameter variation along the coregulation path. The low-spike-frequency bursting, the high-spike-frequency bursting, and plateau-containing oscillations are marked in blue, yellow, and red-brown, respectively. The HCO model exhibits functional bursting over the whole range of parameter variation with intra-burst spike frequency, minimal membrane potential, and amplitude of the intracellular Na⁺ concentration oscillations being notably distinct. In the single HN model, the low-spike-frequency and high-spike-frequency bursting share a trend and transition into pathological seizure-like plateau oscillations in all presented characteristics. Upper left three panels repeat Figures 2C–E for comparison purposes.

Since the intracellular Na⁺ concentration governs the pump current and is in turn affected by the pump current, we assessed the average contribution of currents to Na⁺ flux during the depolarized and hyperpolarized phases of activity for the HCO regimes of interest (Figure 6). Analysis of individual currents shows that I_NaF and I_P are the major drivers of depolarization and intracellular Na⁺ concentration in all regimes of activity. The higher spike-frequency in the high-spike-frequency regime compared to the low-spike-frequency regime brings in much more Na⁺ through I_NaF in the depolarized phase. Higher levels of intracellular Na⁺ increase pump current activation, so the pump is more active and contributes more heavily to the removal of Na⁺.

In the HCO high-spike-frequency regime, the majority of currents have a higher average magnitude during both the depolarized and hyperpolarized phases when compared to low-spike-frequency regime. The high-spike-frequency regime does not appear above an I_PumpMax of 0.42 nA, but the low-spike-frequency regime does appear at high maximal pump currents (i.e., at control conditions). Yet despite only appearing at lower values of I_PumpMax, the HCO high-spike-frequency regime has a higher average pump current throughout the burst than in the control regime (Figure 7). Through coregulation, a reduction in I_PumpMax of 10% from control leads to an 8% increase in average pump current over the burst when the HCO transitions into the high-spike-frequency regime. There are several major exceptions, however, to this idea of higher total current in high-spike-frequency regimes. Persistent Na⁺ current, I_P, is smaller during the hyperpolarized phase of the high-spike-frequency regime than the low spike-frequency regime. Leak current is an inward current on average during the hyperpolarized phase of the high-spike-frequency regime, largely due to the massive afterhyperpolarization associated with this regime. Synaptic current, for the same reasons, appears to be smaller on average during the hyperpolarized phase of the high-spike-frequency regime. In the high-spike-frequency regime, the afterhyperpolarization is great enough to cross the reversal potentials for the leak and synaptic currents and thus briefly reverses their contributions to the current balance on average.

In all regimes, the pump current is the sole source of the outward Na⁺ flux, and so is unsurprisingly critical in maintaining the Na⁺ electrochemical gradient. The major inward contributors of Na⁺ flux during the depolarized phase are I_NaF and I_P, while the major contributors during the hyperpolarized phase are I_P, I_h, and I_Leak. What accounts for the major difference in the amplitude of Na⁺ oscillations between low and high-spike-frequency regimes is principally the I_NaF Na⁺ entry driven by the higher spike frequency.

We next performed a continuous simulation changing parameters ${\overline{g}}_h$and I_PumpMax along the line of coregulation to assess the effects of modulation on a single HN cell (Figure 7) directly. We found that the cell transitions between the low-spike-frequency and high-spike-frequency regimes at ( ${\overline{g}}_h$ , I_PumpMax) = (2.15,0.413), and to the hybrid plateau containing regime at ( ${\overline{g}}_h$ , I_PumpMax) = (2.6,0.406), and then to full plateaus at ( ${\overline{g}}_h$ , I_PumpMax) = (3.0,0.402). The entrance into the domain of multistability by the HCO model approximately coincides with the transition of the single cell into the area producing hybrid plateau-containing activity like that shown in Figure 1B3. Under increasing levels of modulation in the single HN model, period and depolarized phase duration decreased until the transition into plateaus. The HCO model, however, continues to monotonically decrease in period and depolarized phase duration throughout the entire tested modulatory parameter space. Single HN spike frequency increases monotonically through both the low and high-spike-frequency regimes, suggesting that there could be a smooth transition between regimes in this case. This is unlike the HCO model, in which the transition into the high spike-frequency regime is abrupt, and then displays a decrease in spike frequency under progressively higher levels of modulation, i.e., higher values of ${\overline{g}}_h$ and lower I_PumpMax. Yet even in the single cell, the apparent trends of minimum membrane potential in the hyperpolarized interval are different between the low-spike-frequency and high-spike-frequency regimes. This suggests that the two regimes are actually operating using distinct bursting mechanisms, and this is even clearer in the HCO model, in which minimum membrane potential is constant around −65 mV for the low spike-frequency regime, but the minimum membrane potential of the high-spike frequency regime is much lower, and trends upward through the bistable domain from −85 mV to −73 mV under higher levels of modulation. The measured burst characteristics which most clearly differentiate the two regimes are the spike frequency, the minimum membrane potential in the hyperpolarized interval, and the amplitude of [Na]_i oscillations. Thus, results presented in Figure 7 confirms that the major distinctions between the regimes are observed over the ranges of the coregulated parameters.

4 Discussion

Neuromodulation of rhythmic circuits adjusts their patterns of activity to environmental challenges or altering behavioral goals by arranging orchestrated changes of ionic currents. Previously we have shown that application of endogenous leech neuromodulator myomodulin to leech heartbeat half-center oscillators coordinates a decrease of pump current and an increase of h-current and thus adjusts the alternating bursting pattern in a wide range of the cycle period and burst duration (Masino and Calabrese, 2002; Tobin and Calabrese, 2005). In a biophysical model, this coregulation of two currents expands the functional range of the governed biophysical parameters, maximal pump activity and maximal conductance of h-current (Ellingson et al., 2021). Here, we demonstrate that this myomodulin coregulation leads to emergence of a new bursting regime with high spike frequency and roughly the same cycle period. This regime co-exists with a low spike frequency bursting regime along the coregulation of the two currents. Comparing the activity regimes of the HCO and single HN models, we show that neuromodulation induces plateau-containing seizure-like oscillations in the isolated HN neuron. The point of the breakdown into this dysfunctional activity is very close to the point at which modulation gives rise to the coexistent regimes in the HCO model. Thus, the synaptic inhibition harnesses the seizure-like regime into the two functional bursting regimes with disparate spike frequencies. The new high-spike-frequency bursting regime provides an opportunity for a higher level of synaptic input to motor neurons during modulation. Counterintuitively, in the higher spike frequency regime, lowering $I_{\mathit{PumpMax}}$by neuromodulation leads to an increase of the pump generated outward current.

5 Neuromodulation can cause multi-stability of activity regimes

Changes of the biophysical properties that govern neuronal dynamics can result in the co-existence of activity regimes near the transitions between these regimes (Bos et al., 2021). The coexistence is discernible through the identification of state-dependent activity, the observation of hysteresis associated with transitions between different regimes, and the ability to trigger a lasting or permanent switch between the regimes using short pulse perturbations. A number of neuronal models have been shown to exhibit the coexistence of different activity regimes, including silence and spiking, silence and bursting, spiking and bursting, as well as the coexistence of bursting regimes with distinct properties (Lechner et al., 1996; Newman and Butera, 2010; Malashchenko et al., 2011a,b; Marin et al., 2013).

Consistently, neuromodulation-induced changes to neurons can lead to state-dependent transitions (Doi and Ramirez, 2010; Marder et al., 2014; Sharples and Whelan, 2017) and multistability (Hounsgaard and Kiehn, 1989; Lechner et al., 1996). One prominent example is the induction of multistability in neuron R15 within the abdominal ganglia of Aplysia by serotonin (5-HT) (Canavier et al., 1993, 1994; Butera et al., 1995; Butera, 1998). In a biophysical model of R15, under the effects of modeled 5-HT application, the neuron exhibits multistability, producing either co-existing bursting and beating patterns or multiple co-existing bursting patterns, yet under control conditions this neuron exhibits only one bursting regime (Butera et al., 1995). Notably, up to seven bursting regimes were reported to co-exist in the model under the influence of neuromodulation. Model predictions were corroborated in follow-up experiments demonstrating the co-existence of beating and bursting regimes (Lechner et al., 1996) and a series of co-existing bursting regimes differing in the number of spikes by one (Newman and Butera, 2010). The functional role of this multistability in R15 is not well understood.

5.1 Merits of multistability

To generate various rhythmic patterns, a central pattern generator circuit can implement two conceptually different mechanisms: either it moves through a reversible smooth transition between patterns or through state-dependent abrupt transitions exhibiting hysteresis. The latter implies the coexistence of neighboring regimes, enabling controlled switching between them. Central pattern generators in both vertebrates and invertebrates can employ either of these mechanisms, exhibiting gradual changes of the speed of their patterns or transitioning between rhythms characterized by orders-of-magnitude differences in time scales. Examples include in vertebrates dopamine-induced continuous and episodic locomotor bursting regimes of postnatal mouse spinal cord preparation (Sharples and Whelan, 2017; Sharples et al., 2021), intact cat’s paw-shaking and locomotion (Parker et al., 2021), as well as various mammalian breathing patterns (Lieske et al., 2000), and in invertebrates, swimming and crawling behaviors of medicinal leeches and sea slugs (Popescu and Frost, 2002; Briggman and Kristan, 2006), fast and slow swimming of jellyfish (Mackie and Meech, 1985), two distinct gastric mill patterns of the crab stomatogastric nervous system that are induced either by ventral cardiac neurons or by postoesophageal commissure neurons (White and Nusbaum, 2011). These CPG alternative configurations either share a core circuit or share some subset of circuitry and engage other circuit modules in the process of modulation toward transition (Briggman and Kristan, 2008). It is harder to distinguish these basic mechanisms in the vertebrate preparations with larger numbers of neurons. Recently, in zebrafish preparation, a modular organization has been established where the speed of locomotion increases through sequential engagement of corresponding-to-speed modules along with the increase of neuromodulation (Pallucchi et al., 2024).

It is notable that in a CPG with multiple co-existing functional rhythms, no ongoing or permanent changes to the network or the neurons would be required to switch a pattern. A single pulse of conductance representing a single post-synaptic potential can be sufficient to change the activity pattern or induce a transient response on a smaller time-scale like locomotion and paw-shaking in cats (Parker et al., 2018, 2021). The ability to switch between functional patterns with a single pulse is energy efficient, although it requires the neuronal system to be already dynamically capable of producing multiple regimes. With the type of multistability identified here, with the smaller distributed heartbeat CPG (assembled of bilateral pairs of HN cells), our results suggest that a rhythmic neural system could increase spike frequency without changing the underlying cycle period of bursting merely by a single pulse of post-synaptic current. This transition would allow for a more effective synaptic input to motor neurons or muscles controlled by a CPG while maintaining the functional rhythm. Such a regime could exploit muscle dynamics for more efficient or more powerful motor activation depending on the behavioral goals of the animal. While the advantage of greater inhibition of leech heart motoneurons through the HN inhibitory synapses is unclear currently, spike frequency in premotor neurons other motor CPGs has documented effects on the motor pattern (Katz and Frost, 1997; Brezina et al., 2003a,b). We can speculate that high spike frequency bursting will produce more prominent rebound and would be more effective in shaping the segmental motoneuron output pattern. In experiments where 1 μM myomodulin was applied to leech heart interneurons, there an increase in spike frequency ranging between 10 and 30% and reduction of the variability of the pattern (Tobin and Calabrese, 2005) was observed. While the difference in the co-existing regimes in the modulated HCO model is closer to 100%, the mere co-existence of these two particular regimes is novel. They burst on the same time scale, but spike very differently. To our knowledge, no multistability of this specific type has been identified in neural models. This unique form of multistabilty has interesting implications for the study of motor control by CPGs.

5.2 Downregulation of a current can result in larger average current

Another interesting result from this computational study is the observation that a cell model with a lower maximal pump current displayed a larger overall average pump current during the bursting cycle. This means that using an incomplete pharmacological inhibitor does not guarantee a reduction in the current, especially if steady state activity changes. It entirely depends on the dynamics of the other channels involved, which is often unknown in pharmacological experiments on neural circuits. This understanding should be applied to other systems in which the pump is modified. Downregulation or partial blocking of the pump does not guarantee that there will be less pump current. In fact, under certain conditions, it can increase the average pump current.

5.3 Neuroprotective effects of synaptic inhibition

While seizure activity in vertebrates cannot be assigned to a simple change in excitatory/inhibitory balance, general lack of network inhibition is often implicated (Tasker and Dudek, 1991; Calcagnotto et al., 2005; Cossart et al., 2005; Righes Marafiga et al., 2021). There can be no doubt, moreover, that bicuculline and other pharmacological agents that disrupt GABA-A transmission can cause seizures in vivo and in vitro (see for example Chitolina et al., 2023). In leech heart interneurons in intact ganglia, prolonged bicuculline exposure causes seizure like bursts (Cymbalyuk et al., 2002), and blocking of all chemical synaptic transmission (but not electrical) with Co⁺² causes widespread synchronous seizure like bursts including heart interneurons (Angstadt and Friesen, 1991). Here we show that modulation involving the down-regulation of the Na⁺/K⁺ pump can lead to seizure-like activity in HN single cell and HCO models. This is not wholly surprising since the pump contributes outward current at all membrane potentials which is functionally equivalent to synaptic inhibition. Indeed Angstadt and Friesen (1991), provide evidence that the pump is critical to termination of seizure-like bursts induced by Co²⁺ and others have noted the importance of the pump during seizure activity (see for example Krishnan et al., 2015) Along the functional path of modulation of the HCO leading to smooth changes in rhythm period, modulation of the single component cells leads to seizure-like dysfunctional activity. Downmodulation of the pump current ultimately leads to a depolarization block in the single cell, but this is compensated through the synaptic inhibition inherent in the half-center circuit configuration. Cells which always receive some level of inhibition are tuned to function with it, and this dependence on inhibition for stability is one mechanism, among many, which can explain aberrant activity when inhibition is removed. Where one cell is incapable of bursting in a biological parameter range, a network of cells can robustly produce proper activity. A similar idea is discussed in prior work from our research group where the half-center oscillator structure can prevent the effects of an enhanced leak current from disrupting the pattern of activity in leech HN interneurons (Cymbalyuk et al., 2002). Thus, inhibitory circuits may be inherently more robust to modulation than their components neurons.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

PE: Conceptualization, Formal analysis, Investigation, Methodology, Visualization, Writing – original draft, Writing – review & editing, Software. YS: Writing – review & editing, Formal analysis, Investigation, Visualization. JP: Formal analysis, Investigation, Visualization, Writing – review & editing. RC: Funding acquisition, Methodology, Writing – original draft, Writing – review & editing. GC: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. We acknowledge the NIH for funding this work through Grant 1 R21 NS111355 to GC and RC. We acknowledge the Georgia State University Brains and Behavior Program for funding for PE and JP and grant to GC.

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.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

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.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fncel.2024.1395026/full#supplementary-material

References

Alonso, L. M., and Marder, E. (2020). Temperature compensation in a small rhythmic circuit. eLife 9. doi: 10.7554/eLife.55470

Crossref Full Text | Google Scholar

Angstadt, J. D., and Friesen, W. O. (1991). Synchronized oscillatory activity in leech neurons induced by calcium channel blockers. J. Neurophysiol. 66, 1858–1873. doi: 10.1152/jn.1991.66.6.1858

PubMed Abstract | Crossref Full Text | Google Scholar

Bos, R., Drouillas, B., Bouhadfane, M., Pecchi, E., Trouplin, V., Korogod, S. M., et al. (2021). Trpm5 channels encode bistability of spinal motoneurons and ensure motor control of hindlimbs in mice. Nat. Commun. 12:6815. doi: 10.1038/s41467-021-27113-x

PubMed Abstract | Crossref Full Text | Google Scholar

Brezina, V., Horn, C. C., and Weiss, K. R. (2005). Modeling neuromuscular modulation in Aplysia. III. Interaction of central motor commands and peripheral modulatory state for optimal behavior. J. Neurophysiol. 93, 1523–1556. doi: 10.1152/jn.00475.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Brezina, V., Orekhova, I. V., and Weiss, K. R. (2003a). Neuromuscular modulation in Aplysia. I. Dynamic model. J. Neurophysiol. 90, 2592–2612. doi: 10.1152/jn.01091.2002

PubMed Abstract | Crossref Full Text | Google Scholar

Brezina, V., Orekhova, I. V., and Weiss, K. R. (2003b). Neuromuscular modulation in Aplysia. II. Modulation of the neuromuscular transform in behavior. J. Neurophysiol. 90, 2613–2628. doi: 10.1152/jn.01093.2002

PubMed Abstract | Crossref Full Text | Google Scholar

Briggman, K. L., and Kristan, W. B. Jr. (2006). Imaging dedicated and multifunctional neural circuits generating distinct behaviors. J. Neurosci. 26, 10925–10933. doi: 10.1523/JNEUROSCI.3265-06.2006

PubMed Abstract | Crossref Full Text | Google Scholar

Briggman, K. L., and Kristan, W. B. (2008). Multifunctional pattern-generating circuits. Annu. Rev. Neurosci. 31, 271–294. doi: 10.1146/annurev.neuro.31.060407.125552

Crossref Full Text | Google Scholar

Butera, R. J. (1998). Multirhythmic bursting. Chaos 8, 274–284. doi: 10.1063/1.166358

PubMed Abstract | Crossref Full Text | Google Scholar

Butera, R. J., Clark, J. W. Jr., Canavier, C. C., Baxter, D. A., and Byrne, J. H. (1995). Analysis of the effects of modulatory agents on a modeled bursting neuron: dynamic interactions between voltage and calcium dependent systems. J. Comput. Neurosci. 2, 19–44. doi: 10.1007/BF00962706

PubMed Abstract | Crossref Full Text | Google Scholar

Calabrese, R. L. (1998). Cellular, synaptic, network, and modulatory mechanisms involved in rhythm generation. Curr. Opin. Neurobiol. 8, 710–717. doi: 10.1016/S0959-4388(98)80112-8

PubMed Abstract | Crossref Full Text | Google Scholar

Calcagnotto, M. E., Paredes, M. F., Tihan, T., Barbaro, N. M., and Baraban, S. C. (2005). Dysfunction of synaptic inhibition in epilepsy associated with focal cortical dysplasia. J. Neurosci. 25, 9649–9657. doi: 10.1523/JNEUROSCI.2687-05.2005

PubMed Abstract | Crossref Full Text | Google Scholar

Canavier, C. C., Baxter, D. A., Clark, J. W., and Byrne, J. H. (1993). Nonlinear dynamics in a model neuron provide a novel mechanism for transient synaptic inputs to produce long-term alterations of postsynaptic activity. J. Neurophysiol. 69, 2252–2257. doi: 10.1152/jn.1993.69.6.2252

Crossref Full Text | Google Scholar

Canavier, C. C., Baxter, D. A., Clark, J. W., and Byrne, J. H. (1994). Multiple modes of activity in a model neuron suggest a novel mechanism for the effects of neuromodulators. J. Neurophysiol. 72, 872–882. doi: 10.1152/jn.1994.72.2.872

PubMed Abstract | Crossref Full Text | Google Scholar

Chitolina, R., Gallas-Lopes, M., Reis, C. G., Benvenutti, R., Stahlhofer-Buss, T., Calcagnotto, M. E., et al. (2023). Chemically-induced epileptic seizures in zebrafish: a systematic review. Epilepsy Res. 197:107236. doi: 10.1016/j.eplepsyres.2023.107236

PubMed Abstract | Crossref Full Text | Google Scholar

Cossart, R., Bernard, C., and Ben-Ari, Y. (2005). Multiple facets of GABAergic neurons and synapses: multiple fates of GABA signalling in epilepsies. Trends Neurosci. 28, 108–115. doi: 10.1016/j.tins.2004.11.011

PubMed Abstract | Crossref Full Text | Google Scholar

Cressman, J. R. Jr., Ullah, G., Ziburkus, J., Schiff, S. J., and Barreto, E. (2009). The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states: I. Single neuron dynamics. J. Comput. Neurosci. 26, 159–170. doi: 10.1007/s10827-008-0132-4

PubMed Abstract | Crossref Full Text | Google Scholar

Cymbalyuk, G. S., and Calabrese, R. L. (2000). Oscillatory behaviors in pharmacologically isolated heart interneurons from the medicinal leech. Neurocomputing 32, 97–104. doi: 10.1016/S0925-2312(00)00149-1

Crossref Full Text | Google Scholar

Cymbalyuk, G. S., Gaudry, Q., Masino, M. A., and Calabrese, R. L. (2002). Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J. Neurosci. Off. J. Soc. Neurosci. 22, 10580–10592. doi: 10.1523/JNEUROSCI.22-24-10580.2002

PubMed Abstract | Crossref Full Text | Google Scholar

Doi, A., and Ramirez, J. M. (2008). Neuromodulation and the orchestration of the respiratory rhythm. Respir. Physiol. Neurobiol. 164, 96–104. doi: 10.1016/j.resp.2008.06.007

PubMed Abstract | Crossref Full Text | Google Scholar

Doi, A., and Ramirez, J. M. (2010). State-dependent interactions between excitatory neuromodulators in the neuronal control of breathing. J. Neurosci. 30, 8251–8262. doi: 10.1523/JNEUROSCI.5361-09.2010

PubMed Abstract | Crossref Full Text | Google Scholar

Ellingson, P. J., Barnett, W. H., Kueh, D., Vargas, A., Calabrese, R. L., and Cymbalyuk, G. S. (2021). Comodulation of h- and Na(+)/K(+) pump currents expands the range of functional bursting in a central pattern generator by navigating between dysfunctional regimes. J. Neurosci. 41, 6468–6483. doi: 10.1523/JNEUROSCI.0158-21.2021

PubMed Abstract | Crossref Full Text | Google Scholar

Erazo-Toscano, R., Fomenko, M., Core, S., Calabrese, R. L., and Cymbalyuk, G. (2023). Bursting dynamics based on the persistent Na(+) and Na(+)/K(+) pump currents: a dynamic clamp approach. eNeuro 10:331. doi: 10.1523/ENEURO.0331-22.2023

PubMed Abstract | Crossref Full Text | Google Scholar

Fall, C. P., Lewis, T. J., and Rinzel, J. (2005). Background-activity-dependent properties of a network model for working memory that incorporates cellular bistability. Biol. Cybern. 93, 109–118. doi: 10.1007/s00422-005-0543-5

PubMed Abstract | Crossref Full Text | Google Scholar

Forrest, M. D., Wall, M., Press, D. P., and Feng, J. F. (2012). The sodium-potassium pump controls the intrinsic firing of the cerebellar Purkinje neuron. Plo S one 7:e51169. doi: 10.1371/journal.pone.0051169

PubMed Abstract | Crossref Full Text | Google Scholar

Frohlich, F., and Bazhenov, M. (2006). Coexistence of tonic firing and bursting in cortical neurons. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 74:031922. doi: 10.1103/PhysRevE.74.031922

PubMed Abstract | Crossref Full Text | Google Scholar

Goldman, M. S., Golowasch, J., Marder, E., and Abbott, L. F. (2001). Global structure, robustness, and modulation of neuronal models. J. Neurosci. 21, 5229–5238. doi: 10.1523/JNEUROSCI.21-14-05229.2001

PubMed Abstract | Crossref Full Text | Google Scholar

Guttman, R., Lewis, S., and Rinzel, J. (1980). Control of repetitive firing in squid axon membrane as a model for a neuroneoscillator. J. Physiol. 305, 377–395. doi: 10.1113/jphysiol.1980.sp013370

PubMed Abstract | Crossref Full Text | Google Scholar

Hahn, P. J., and Durand, D. M. (2001). Bistability dynamics in simulations of neural activity in high-extracellular-potassium conditions. J. Comput. Neurosci. 11, 5–18. doi: 10.1023/A:1011250329341

Crossref Full Text | Google Scholar

Harris-Warrick, R. M., and Johnson, B. R. (2010). Checks and balances in neuromodulation. Front. Behav. Neurosci. 4:47. doi: 10.3389/fnbeh.2010.00047

Crossref Full Text | Google Scholar

Hounsgaard, J., Hultborn, H., Jespersen, B., and Kiehn, O. (1984). Intrinsic membrane properties causing a bistable behaviour of alpha-motoneurones. Exp. Brain Res. 55, 391–394. doi: 10.1007/BF00237290

Crossref Full Text | Google Scholar

Hounsgaard, J., and Kiehn, O. (1989). Serotonin-induced bistability of turtle motoneurones caused by a nifedipine-sensitive calcium plateau potential. J. Physiol. 414, 265–282. doi: 10.1113/jphysiol.1989.sp017687

PubMed Abstract | Crossref Full Text | Google Scholar

Katz, P. S., and Frost, W. N. (1997). Removal of spike frequency adaptation via neuromodulation intrinsic to the Tritonia escape swim central pattern generator. J. Neurosci. 17, 7703–7713. doi: 10.1523/JNEUROSCI.17-20-07703.1997

PubMed Abstract | Crossref Full Text | Google Scholar

Krishnan, G. P., Filatov, G., Shilnikov, A., and Bazhenov, M. (2015). Electrogenic properties of the Na(+)/K(+) ATPase control transitions between normal and pathological brain states. J. Neurophysiol. 113, 3356–3374. doi: 10.1152/jn.00460.2014

PubMed Abstract | Crossref Full Text | Google Scholar

Kueh, D., Barnett, W. H., Cymbalyuk, G. S., and Calabrese, R. L. (2016). Na(+)/K(+) pump interacts with the h-current to control bursting activity in central pattern generator neurons of leeches. eLife 5:322. doi: 10.7554/eLife.19322

PubMed Abstract | Crossref Full Text | Google Scholar

Lechner, H. A., Baxter, D. A., Clark, J. W., and Byrne, J. H. (1996). Bistability and its regulation by serotonin in the endogenously bursting neuron R15 in Aplysia. J. Neurophysiol. 75, 957–962. doi: 10.1152/jn.1996.75.2.957

PubMed Abstract | Crossref Full Text | Google Scholar

Lieske, S. P., Thoby-Brisson, M., Telgkamp, P., and Ramirez, J. M. (2000). Reconfiguration of the neural network controlling multiple breathing patterns: eupnea, sighs and gasps. Nat. Neurosci. 3, 600–607. doi: 10.1038/75776

PubMed Abstract | Crossref Full Text | Google Scholar

Loewenstein, Y., Mahon, S., Chadderton, P., Kitamura, K., Sompolinsky, H., Yarom, Y., et al. (2005). Bistability of cerebellar Purkinje cells modulated by sensory stimulation. Nat. Neurosci. 8, 202–211. doi: 10.1038/nn1393

PubMed Abstract | Crossref Full Text | Google Scholar

Mackie, G. O., and Meech, R. W. (1985). Separate sodium and calcium spikes in the same axon. Nature 313, 791–793. doi: 10.1038/313791a0

PubMed Abstract | Crossref Full Text | Google Scholar

Malashchenko, T., Shilnikov, A., and Cymbalyuk, G. (2011a). Bistability of bursting and silence regimes in a model of a leech heart interneuron. Phys. Rev. E Stat. Nonlinear Soft Matter Phys. 84:041910. doi: 10.1103/PhysRevE.84.041910

Crossref Full Text | Google Scholar

Malashchenko, T., Shilnikov, A., and Cymbalyuk, G. (2011b). Six types of multistability in a neuronal model based on slow calcium current. PLoS One 6:e21782. doi: 10.1371/journal.pone.0021782

PubMed Abstract | Crossref Full Text | Google Scholar

Marder, E. (2012). Neuromodulation of neuronal circuits: back to the future. Neuron 76, 1–11. doi: 10.1016/j.neuron.2012.09.010

PubMed Abstract | Crossref Full Text | Google Scholar

Marder, E., and Calabrese, R. L. (1996). Principles of rhythmic motor pattern generation. Physiol. Rev. 76, 687–717. doi: 10.1152/physrev.1996.76.3.687

Crossref Full Text | Google Scholar

Marder, E., O'Leary, T., and Shruti, S. (2014). Neuromodulation of circuits with variable parameters: single neurons and small circuits reveal principles of state-dependent and robust neuromodulation. Annu. Rev. Neurosci. 37, 329–346. doi: 10.1146/annurev-neuro-071013-013958

PubMed Abstract | Crossref Full Text | Google Scholar

Marin, B., Barnett, W. H., Doloc-Mihu, A., Calabrese, R. L., and Cymbalyuk, G. S. (2013). High prevalence of multistability of rest states and bursting in a database of a model neuron. PLoS Comput. Biol. 9:e1002930. doi: 10.1371/journal.pcbi.1002930

PubMed Abstract | Crossref Full Text | Google Scholar

Masino, M. A., and Calabrese, R. L. (2002). Period differences between segmental oscillators produce intersegmental phase differences in the leech heartbeat timing network. J. Neurophysiol. 87, 1603–1615. doi: 10.1152/jn.00338.2001

PubMed Abstract | Crossref Full Text | Google Scholar

Nadim, F., Brezina, V., Destexhe, A., and Linster, C. (2008). State dependence of network output: modeling and experiments. J. Neurosci. 28, 11806–11813. doi: 10.1523/JNEUROSCI.3796-08.2008

PubMed Abstract | Crossref Full Text | Google Scholar

Newman, J. P., and Butera, R. J. (2010). Mechanism, dynamics, and biological existence of multistability in a large class of bursting neurons. Chaos 20:023118. doi: 10.1063/1.3413995

PubMed Abstract | Crossref Full Text | Google Scholar

Olsen, O. H., and Calabrese, R. L. (1996). Activation of intrinsic and synaptic currents in leech heart interneurons by realistic waveforms. J. Neurosci. Off. J. Soc. Neurosci. 16, 4958–4970. doi: 10.1523/JNEUROSCI.16-16-04958.1996

PubMed Abstract | Crossref Full Text | Google Scholar

Olypher, A., Cymbalyuk, G., and Calabrese, R. L. (2006). Hybrid systems analysis of the control of burst duration by low-voltage-activated calcium current in leech heart interneurons. J. Neurophysiol. 96, 2857–2867. doi: 10.1152/jn.00582.2006

PubMed Abstract | Crossref Full Text | Google Scholar

Opdyke, C. A., and Calabrese, R. L. (1994). A persistent sodium current contributes to oscillatory activity in heart interneurons of the medicinal leech. J. Comp. Physiol. A Sens. Neural Behav. Physiol. 175, 781–789. doi: 10.1007/BF00191850

PubMed Abstract | Crossref Full Text | Google Scholar

Pallucchi, I., Bertuzzi, M., Madrid, D., Fontanel, P., Higashijima, S. I., and El Manira, A. (2024). Molecular blueprints for spinal circuit modules controlling locomotor speed in zebrafish. Nat. Neurosci. 27, 78–89. doi: 10.1038/s41593-023-01479-1

PubMed Abstract | Crossref Full Text | Google Scholar

Parker, J., Bondy, B., Prilutsky, B. I., and Cymbalyuk, G. (2018). Control of transitions between locomotor-like and paw shake-like rhythms in a model of a multistable central pattern generator. J. Neurophysiol. 120, 1074–1089. doi: 10.1152/jn.00696.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Parker, J. R., Klishko, A. N., Prilutsky, B. I., and Cymbalyuk, G. S. (2021). Asymmetric and transient properties of reciprocal activity of antagonists during the paw-shake response in the cat. PLoS Comput. Biol. 17:e1009677. doi: 10.1371/journal.pcbi.1009677

PubMed Abstract | Crossref Full Text | Google Scholar

Paydarfar, D., Forger, D. B., and Clay, J. R. (2006). Noisy inputs and the induction of on-off switching behavior in a neuronal pacemaker. J. Neurophysiol. 96, 3338–3348. doi: 10.1152/jn.00486.2006

Crossref Full Text | Google Scholar

Perrier, J. F., and Hounsgaard, J. (2000). Development and regulation of response properties in spinal cord motoneurons. Brain Res. Bull. 53, 529–535. doi: 10.1016/S0361-9230(00)00386-5

PubMed Abstract | Crossref Full Text | Google Scholar

Popescu, I. R., and Frost, W. N. (2002). Highly dissimilar behaviors mediated by a multifunctional network in the marine mollusk Tritonia diomedea. J. Neurosci. 22, 1985–1993. doi: 10.1523/JNEUROSCI.22-05-01985.2002

PubMed Abstract | Crossref Full Text | Google Scholar

Prinz, A. A., Billimoria, C. P., and Marder, E. (2003). Alternative to hand-tuning conductance-based models: construction and analysis of databases of model neurons. J. Neurophysiol. 90, 3998–4015. doi: 10.1152/jn.00641.2003

PubMed Abstract | Crossref Full Text | Google Scholar

Prinz, A. A., Bucher, D., and Marder, E. (2004). Similar network activity from disparate circuit parameters. Nat. Neurosci. 7, 1345–1352. doi: 10.1038/nn1352

PubMed Abstract | Crossref Full Text | Google Scholar

Righes Marafiga, J., Vendramin Pasquetti, M., and Calcagnotto, M. E. (2021). GABAergic interneurons in epilepsy: more than a simple change in inhibition. Epilepsy Behav. 121:106935. doi: 10.1016/j.yebeh.2020.106935

PubMed Abstract | Crossref Full Text | Google Scholar

Sharples, S. A., Parker, J., Vargas, A., Milla-Cruz, J. J., Lognon, A. P., Cheng, N., et al. (2021). Contributions of h- and Na(+)/K(+) pump currents to the generation of episodic and continuous rhythmic activities. Front. Cell. Neurosci. 15:715427. doi: 10.3389/fncel.2021.715427

Crossref Full Text | Google Scholar

Sharples, S. A., and Whelan, P. J. (2017). Modulation of rhythmic activity in mammalian spinal networks is dependent on excitability state. eNeuro 4:368. doi: 10.1523/ENEURO.0368-16.2017

PubMed Abstract | Crossref Full Text | Google Scholar

Sorensen, M., DeWeerth, S., Cymbalyuk, G., and Calabrese, R. L. (2004). Using a hybrid neural system to reveal regulation of neuronal network activity by an intrinsic current. J. Neurosci. Off. J. Soc. Neurosci. 24, 5427–5438. doi: 10.1523/JNEUROSCI.4449-03.2004

PubMed Abstract | Crossref Full Text | Google Scholar

Stadele, C., and Stein, W. (2022). Neuromodulation enables temperature robustness and coupling between fast and slow oscillator circuits. Front. Cell. Neurosci. 16:849160. doi: 10.3389/fncel.2022.849160

PubMed Abstract | Crossref Full Text | Google Scholar

Tasker, J. G., and Dudek, F. E. (1991). Electrophysiology of GABA-mediated synaptic transmission and possible roles in epilepsy. Neurochem. Res. 16, 251–262. doi: 10.1007/BF00966088

PubMed Abstract | Crossref Full Text | Google Scholar

Tobin, A.-E., and Calabrese, R. L. (2005). Myomodulin increases Ih and inhibits the NA/K pump to modulate bursting in leech heart interneurons. J. Neurophysiol. 94, 3938–3950. doi: 10.1152/jn.00340.2005

PubMed Abstract | Crossref Full Text | Google Scholar

Tobin, A. E., and Calabrese, R. L. (2006). Endogenous and half-center bursting in morphologically inspired models of leech heart interneurons. J. Neurophysiol. 96, 2089–2106. doi: 10.1152/jn.00025.2006

PubMed Abstract | Crossref Full Text | Google Scholar

Toscano, R. J. E., Ellingson, P. J., Calabrese, R. L., and Cymbalyuk, G. S. (2021). Contribution of the Na+/K+ pump to rhythmic bursting, explored with modeling and dynamic clamp analyses. J. Vis. Exp. 171:1473. doi: 10.3791/61473

Crossref Full Text | Google Scholar

Turrigiano, G. G., Marder, E., and Abbott, L. F. (1996). Cellular short-term memory from a slow potassium conductance. J. Neurophysiol. 75, 963–966. doi: 10.1152/jn.1996.75.2.963

PubMed Abstract | Crossref Full Text | Google Scholar

Ullah, G., Cressman, J. R., Barreto, E., and Schiff, S. J. (2009). The influence of sodium and potassium dynamics on excitability, seizures, and the stability of persistent states. II. Network and glial dynamics. J. Comput. Neurosci. 26, 171–183. doi: 10.1007/s10827-008-0130-6

PubMed Abstract | Crossref Full Text | Google Scholar

Weaver, A. L., Roffman, R. C., Norris, B. J., and Calabrese, R. L. (2010). A role for compromise: synaptic inhibition and electrical coupling interact to control phasing in the leech heartbeat CpG. Front. Behav. Neurosci. 4:38. doi: 10.3389/fnbeh.2010.00038

Crossref Full Text | Google Scholar

Wenning, A., Norris, B. J., Gunay, C., Kueh, D., and Calabrese, R. L. (2018). Output variability across animals and levels in a motor system. eLife 7:31123. doi: 10.7554/eLife.31123

PubMed Abstract | Crossref Full Text | Google Scholar

White, R. S., and Nusbaum, M. P. (2011). The same core rhythm generator underlies different rhythmic motor patterns. J. Neurosci. 31, 11484–11494. doi: 10.1523/JNEUROSCI.1885-11.2011

PubMed Abstract | Crossref Full Text | Google Scholar

Williams, S. R., Christensen, S. R., Stuart, G. J., and Hausser, M. (2002). Membrane potential bistability is controlled by the hyperpolarization-activated current I (H) in rat cerebellar Purkinje neurons in vitro. J. Physiol. 539, 469–483. doi: 10.1113/jphysiol.2001.013136

PubMed Abstract | Crossref Full Text | Google Scholar

Ziburkus, J., Cressman, J. R., Barreto, E., and Schiff, S. J. (2006). Interneuron and pyramidal cell interplay during in vitro seizure-like events. J. Neurophysiol. 95, 3948–3954. doi: 10.1152/jn.01378.2005

PubMed Abstract | Crossref Full Text | Google Scholar