Featured ArticleResearch Articles, Systems/Circuits

Comodulation of h- and Na+/K+ Pump Currents Expands the Range of Functional Bursting in a Central Pattern Generator by Navigating between Dysfunctional Regimes

Parker J. Ellingson, William H. Barnett, Daniel Kueh, Alex Vargas, Ronald L. Calabrese and Gennady S. Cymbalyuk
Journal of Neuroscience 28 July 2021, 41 (30) 6468-6483; https://doi.org/10.1523/JNEUROSCI.0158-21.2021

Abstract

Central pattern generators (CPGs), specialized oscillatory neuronal networks controlling rhythmic motor behaviors such as breathing and locomotion, must adjust their patterns of activity to a variable environment and changing behavioral goals. Neuromodulation adjusts these patterns by orchestrating changes in multiple ionic currents. In the medicinal leech, the endogenous neuromodulator myomodulin speeds up the heartbeat CPG by reducing the electrogenic Na+/K+ pump current and increasing h-current in pairs of mutually inhibitory leech heart interneurons (HNs), which form half-center oscillators (HN HCOs). Here we investigate whether the comodulation of two currents could have advantages over a single current in the control of functional bursting patterns of a CPG. We use a conductance-based biophysical model of an HN HCO to explain the experimental effects of myomodulin. We demonstrate that, in the model, comodulation of the Na+/K+ pump current and h-current expands the range of functional bursting activity by avoiding transitions into nonfunctional regimes, such as asymmetric bursting and plateau-containing seizure-like activity. We validate the model by finding parameters that reproduce temporal bursting characteristics matching experimental recordings from HN HCOs under control, three different myomodulin concentrations, and Cs+ treated conditions. The matching cases are located along the border of an asymmetric regime away from the border with more dangerous seizure-like activity. We found a simple comodulation mechanism with an inverse relation between the pump and h-currents makes a good fit of the matching cases and comprises a general mechanism for the robust and flexible control of oscillatory neuronal networks.

SIGNIFICANCE STATEMENT Rhythm-generating neuronal circuits adjust their oscillatory patterns to accommodate a changing environment through neuromodulation. In different species, chemical messengers participating in such processes may target two or more membrane currents. In medicinal leeches, the neuromodulator myomodulin speeds up the heartbeat central pattern generator by reducing Na+/K+ pump current and increasing h-current. In a computational model, we show that this comodulation expands the range of central pattern generator's functional activity by navigating the circuit between dysfunctional regimes resulting in a much wider range of cycle period. This control would not be attainable by modulating only one current, emphasizing the synergy of combined effects. Given the prevalence of h-current and Na+/K+ pump current in neurons, similar comodulation mechanisms may exist across species.

Introduction

To achieve behavioral flexibility necessary for survival in a variable environment, neuronal circuits must produce, control, and scale functional activity over a broad range of underlying biophysical properties (Marder and Calabrese, 1996; Calabrese, 1998; Doi and Ramirez, 2008; Marder, 2012). Central pattern generators (CPGs) and other rhythmic neuronal circuits are continuously adjusted through neuromodulation, which is orchestrated by changes in multiple ionic currents, including the electrogenic Na+/K+ pump (Marder and Calabrese, 1996; Doi and Ramirez, 2008; Grashow et al., 2010; Harris-Warrick, 2011; Marder, 2012; Marder et al., 2014; Sharples and Whelan, 2017). Yet we are only beginning to understand how the phenomenal robustness of neurons and networks is attained given the wide range of functional burst properties produced by modulation (Grashow et al., 2009; Tang et al., 2012; Marder et al., 2015; Dashevskiy and Cymbalyuk, 2018; Parker et al., 2018). Because neuromodulators or cocktails of neuromodulators may have multiple cellular targets, it is often difficult to determine the combined effects of modulation by sampling effects of single-factor variation, such as dose–response experiments. For example, some modulators, such as dopamine or serotonin, commonly have opposing effects (Seamans et al., 2001; Teshiba et al., 2001; Sharples and Whelan, 2017). Opposing effects may be distinguished by timescales and dose dependence, producing unique transient and steady-state effects (Harris-Warrick et al., 1998; Brezina et al., 2003; Spitzer et al., 2008) and may enhance flexibility through state dependence or protect against overmodulation by incorporating negative feedback control mechanisms (Harris-Warrick and Johnson, 2010). In contrast, effects of the endogenous neuromodulator myomodulin on the leech heartbeat CPG are caused by changes in two membrane currents with superficially similar effects; it decreases burst period by increasing h-current (Ih) and decreasing Na+/K+ pump current (IPump) (Tobin and Calabrese, 2005). How modulators coordinate these effects on activity patterns is poorly understood. Evaluating the advantages offered by the comodulation of two currents with similar effects and the special role of IPump in modulation requires comprehensive accounting using accurate biophysical modeling.

IPump is a unique modulation target (Bertorello and Aperia, 1990; Bertorello et al., 1990; Catarsi and Brunelli, 1991; Scuri et al., 2007; Hazelwood et al., 2008; H. Y. Zhang and Sillar, 2012; L. N. Zhang et al., 2012, 2013; H. Y. Zhang et al., 2015). While maintaining the proper gradients of Na+ and K+ concentrations across the membrane, it produces an outward current with activation kinetics, which do not depend on the membrane potential but instead are governed by the internal Na+ concentration ([Na+]i). Thus, IPump can play an important role in the bursting activity (Li et al., 1996; Frohlich et al., 2006; Arganda et al., 2007; Pulver and Griffith, 2010; Barreto and Cressman, 2011; Krishnan and Bazhenov, 2011; Yu et al., 2012; Jasinski et al., 2013). IPump naturally interacts with Ih: it is most active during depolarized phase because of raised [Na+]i and remains active during the hyperpolarized phase, where IPump is opposed by Ih (Forrest et al., 2012; Y. Zhang et al., 2017).

The interaction of IPump and Ih plays a prominent role in the leech heartbeat CPG. This CPG encompasses two heart interneuron half-center oscillators (HN HCOs): pairs of mutually inhibitory HN interneurons producing alternating bursting activity. IPump and Ih interactions were revealed using the H+/Na+ antiporter monensin, which stimulates IPump by increasing [Na+]i (Kueh et al., 2016). In the presence of Ih, application of monensin dramatically decreases the burst duration (BD) and interburst interval (IBI). If Ih is blocked, then monensin decreases BD but lengthens IBI. This experiment demonstrates how each phase of bursting, BD and IBI, is independently affected by IPump depending on its interaction with Ih. The speedup of the HN HCO burst period by myomodulin (Tobin and Calabrese, 2005), by increasing Ih and decreasing IPump in HNs, suggests coordinated participation of IPump and Ih and demonstrates their importance as a joint modulatory target controlling burst characteristics.

Here we investigated how comodulation of IPump and Ih contributes to flexible control of functional bursting activity.

Materials and Methods

Model

We optimized a single-compartment Hodgkin-Huxley style (Hodgkin and Huxley, 1952) model of a single leech heart interneuron (HN) and an HCO from Kueh et al. (2016). The HN HCO comprises two identical HN neurons mutually inhibiting each other (see Fig. 1). Each neuron contains a leak current (ILeak), a IPump governed by intracellular [Na+]i,, and eight voltage-gated currents: fast Na+ (INaF), persistent Na+ (IP), fast low-voltage activated Ca2+ (ICaF), low-voltage activated slowly inactivating Ca2+ (ICaS), hyperpolarization-activated (h-) (Ih), slowly inactivating delayed-rectifier K+ (IK1), persistent K+ (IK2), and fast A-type K+ (IKA) currents. The leak current and Ih each have two components: one carrying Na+, ILeak,Na, and Ih,Na, and one carrying K+, ILeak,K, and Ih,K. The current conservation equation governing membrane potential dynamics used for each neuron is as follows: CmdVdt=−(INaF + IP + ILeak,Na + Ih,Na + ICaF + ICaS + IK1 + IK2 + IKA  + ILeak,K + Ih,K + IPump + Isyn)(1)

The individual currents were computed in a conductance-based manner as follows: Ix=g¯xmxaxhxbx(V−Ex)(2) dmxdt=mx∞(V)−mτmx(V)(3) dhxdt=hx∞(V)−hτhx(V)(4) where g¯x and Ex are the maximal conductance and reversal potential, respectively, of x current where x∈{NaF,P,CaF,CaS,h,K1,K2,KA}. The state variables mx and hx are the current's activation and inactivation voltage-gating variables, respectively; mx∞(V) and hx∞(V) are its steady-state functions; and τmx(V) and τhx(V) are its time constant functions. Current parameters and equations for these functions are available in Tables 1–3.

The Na+ and K+ components of the leak current and Ih were separately evaluated and Na+ components were taken into account in the dynamics of [Na+]i as follows: Ih,Na=37g¯hmh2(V−ENa)(5) Ih,K=47g¯hmh2(V−EK)(6) ILeak,Na=gLeak,Na(V−ENa)(7) ILeak,K=gLeak,K(V−EK)(8)

Components of the leak conductance were determined in terms of the conventional equivalent leak conductance and leak reversal potential, referred to as reference values, into two separate specific Na+ and K+ leak currents so that: ILeak=gLeak(V−ELeak,Ref)=ILeak,Na+ILeak,K(9)

The resulting partial conductances are as follows: gLeak,Na=gLeakELeak,Ref−EKENa,Ref−EK(10) gLeak,K=gLeakELeak,Ref−ENa,RefEK−ENa,Ref(11)

The ENa,Ref was fixed at 0.045 V, and these partial conductances were fixed before simulation and were kept the same for the whole study.

Intracellular [Na+]i is computed dynamically, taking into account the Na+ fluxes generated by two Na+ currents, INaF and IP, Na+ components of leak and Ih, and the IPump as follows: d[Na]idt=−INaF + IP + Ih,Na + ILeak,Na + 3IPumpvF(12)

[Na+]i also determines the Na+ reversal potential as follows: ENa=RTFln([Na]o[Na]i)(13)

Above, F is Faraday's constant, R is the gas constant, T is the temperature in Kelvin, v is the volume of the intracellular Na+ compartment, and [Na]o is the extracellular [Na+].

Pump current is governed by [Na]i as follows: IPump=IPumpMax1+exp([Na]ih−[Na]i[Na]is)(14) where [Na]ih is [Na]i of half-activation of the pump current and [Na]is is the pump sensitivity to [Na]i.

Synaptic Cl– currents have spike-dependent and graded components as follows: Isyn=ISynS + ISynG(15) ISynS=g¯SynSYPreMPre(V−ESyn)(16) ISynG=g¯SynGPPre3C + PPre3(V−ESyn)(17) where C=10−32 C3 and ESyn=−0.0625V.

Spike-dependent synapse activation was computed as follows: dYdt=X−Y0.011(18) dXdt=X∞(V)−X0.002(19) dMdt=M∞(V)−M0.2(20) X∞(V)=11 + exp(−1000(V + 0.01))(21) M∞(V)=0.1 + 0.91 + exp(−1000(V + 0.04))(22)

Graded synaptic activation was computed using a proxy of presynaptic Ca2+ concentration in the following manner: dPdt=ICa−BP(23) where B is a Ca2+ buffering rate constant (B=10s−1), and ICa is computed as follows: ICa=max(0,(−ICaF−ICaS−A)/109) dAdt=A∞(V)−A0.2(25) A∞(V)=10−101 + exp(–100(V + 0.02))(26)

We integrated the differential equations describing the model dynamics using the Embedded Runge-Kutta Prince-Dormand (8, 9) method from the gsl_odeiv module in the GNU Scientific Library version 2.6 (https://www.gnu.org/software/gsl/). Absolute tolerance was set at 10−9, relative tolerance was set at 10−10, and maximal time step was set to 10−3 s. Models were implemented in the C programming language and compiled using the GNU Compiler Collection version 7.5.0 (https://gcc.gnu.org/).

View this table:
Table 1.

Cell constants and current parameters conductance g¯x(nS), activation exponential ax, inactivation exponential bx, and reversal potential Ex (V)

Parameter sweeps

We characterized and mapped activity regimes of the model under variation of two parameters, the maximal Ih conductance (g¯h) and the maximal IPump (IPumpMax). For this, we fix g¯h to a certain value in the range starting with 0 nS, and sweep IPumpMax values, starting from 0.5 nA and decreasing it by a step value of 0.001 nA each next simulation until a value of 0.3 nA was reached. Initial conditions for each simulation for a given g¯h along decreasing IPumpMax were taken from the end of the previous simulation at the previous value of IPumpMax. Then we incremented g¯h by a step value of 0.2 nS and fixed g¯h for the next sweep of IPumpMax. For each new value of g¯h, a standard set of initial conditions was used (Table 4). This procedure was repeated until all parameter values with the aforementioned step sizes were taken for [0-10] nS and [0.3-0.5] nA ranges for g¯h and IPumpMax, respectively. At each parameter pair, the system was simulated for 1600 s. The first 1300 s containing the transient part of the trace was removed, and only the last 300 s was used in the analysis to characterize the steady-state activity.

Analysis of temporal characteristics of model activity regimes

Analysis of the activity of the model HCO was conducted using custom scripts developed in the MATLAB software (The MathWorks, Inc.). In our parameter sweeps, we observed a variety of different types of activity, 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. We defined functional bursting as symmetrical activity in which cells alternate in generating bursts of spikes; during the burst, cell's membrane potential would rise above −45 mV, and continuously spike until its membrane potential dropped below −45 mV. Asymmetric bursting activity was defined as a similar activity, but with alternating bursts in which bursts of one cell were noticeably longer than those of the other cell. Finally, plateau-containing activity was characterized by the appearance of plateau events in the HCO's activity. We defined plateau-containing activity as an oscillatory activity in which at least one cell in the HCO exhibited at least one depolarized interval in which the cell's membrane potential would rise above −45 mV, but spiking within such interval failed. At some parameter domains, plateau activity could be periodic and consist purely of alternating plateau events, while at other domains, activity could be a mixture of plateau events with conventional bursts.

For automated detection and analysis of the three activity regimes defined above, we specify two terms. The depolarized phase of activity of a neuron is the time interval in which membrane potential rose above −45 mV and remained there for at least 0.5 s. Correspondingly, the hyperpolarized phase of the activity was the interval in which membrane potential was below −45 mV. The duration of these phases provides a measure of the envelope of bursting in functional and asymmetrical regimes and also allows us to uniformly measure the temporal characteristics of plateau events using the same method of analysis.

We define a burst event as a continuous spiking interval during a depolarized phase of activity. Spike times were identified as moments of time at which the peaks in the voltage trace were above a threshold of −30 mV. Spike times were then used to distinguish bursts from plateaus. Spikes that had no predecessor or occurred at least 0.4 s after their predecessor were marked as the beginning of a spike-train. Then, the end of the spike-train was determined by spikes that had no successive spikes for at least 0.4 s. If multiple spike-trains were detected within a single depolarized phase or if the end of the depolarized phase did not have a preceding spike within 0.4 s, that depolarized phase would be tagged as a plateau event. A depolarized phase, which contained a single uninterrupted spike-train from the beginning to the end of the depolarized interval, was then tagged as a burst event. In depolarized phases that were classified as bursts, BD was computed as the time between first and last spike detected, burst period was computed as the time between the first spike in a burst and the first spike in the next burst, and IBI was computed as the time between the last spike in a burst and the first spike in the next burst. In cases where the depolarized phase was defined as a burst, the depolarized phase duration was nearly identical (±1 interspike interval) to BD, making this a good metric for describing the temporal pattern of either bursts or plateaus.

We define a measure of asymmetry between the depolarized phases of each side (HN(R) and HN(L)) of the HCO using the following expression: Asymm=2|UD¯HN,R−UD¯HN,L|UD¯HN,R + UD¯HN,L(27) where UD¯HN,R and UD¯HN,L are the mean depolarized phase durations of model cells HN(R) and HN(L), respectively. Asymmetric bursting regimes were defined as those bursting regimes which did not exhibit any plateau events and for which the asymmetry measure was >0.2. This threshold measure is represented by a case in which the BD of one model neuron is at least 55% of the period and the BD of the other model neuron is at most 45% of the period. Correspondingly, a bursting regime, which did not contain any plateau events, had an asymmetry measure <0.2, and had cycle period within the current and previous experimental ranges (Tobin and Calabrese, 2005; Kueh et al., 2016; Wenning et al., 2018) were labeled as functional.

Dimensions of the area supporting functional bursting

We investigated the functional regime along each of the g¯h and IPumpMax axes. The n × m (51 × 201) sweep of parameter space was divided into continuous 1 × m and 1 × n slices in which each slice represents the variation of either g¯h or IPumpMax while the other is held constant. At each slice, the boundaries of the functional regime were identified, and their difference was calculated as the local range. Then, for each axis, the maximal possible parameter range retaining functional activity was determined. BDs, burst periods, and IBIs at and within these boundaries were also investigated, and their maximal possible ranges were determined. These maximal ranges were then compared with the maximal range achievable through comodulation, where the axis of comodulation was an inverse function fitted to all parameter pairs in the functional regime (see Curve fitting).

Experimental design and statistical analysis

We conducted two sets of experiments, each with a within-subjects (repeated-measures) design, whereby we would determine the dose-dependent effects of myomodulin with and without 2 mm Cs+ (Ih blocker) applied. Each set of experiments was conducted with 5 medicinal leeches, Hirudo spp., using the electrophysiological protocol outlined below. For both experiments, we analyzed the mean burst period, coefficient of variation of burst period, and mean BD of the HN(3,4) neurons using one-way repeated-measures ANOVAs with treatment as the independent variable. Post hoc pairwise comparisons were then conducted using Tukey's Honest Significant Difference test to determine which of the treatment conditions were significantly different from control or Cs+ alone and whether there was a dose effect between treatments in each experiment. All burst statistics for experiments are reported with the grand mean and SEM. p < 0.05 was regarded as statistically significant for all statistical tests.

Electrophysiology

Experimental protocols were similar to those described in Kueh et al. (2016). After surgically isolating the mid-body ganglion 3 or 4 of the nerve cord of medicinal leeches (n = 5) and removing the perineurium, extracellular suction electrodes were used to record from left and right HN(3,4) neurons. After allowing the preparation to settle, control recordings were taken. Subsequently, three concentrations of myomodulin (1, 10, and 100 μm) were applied with a 10 min wash-out between each application. Each application lasted between 5 and 10 min, and short recordings (containing 10-14 bursting cycles) were taken at each concentration once the preparation settled into a regular rhythm. We also performed a second set of experiments on a separate group of animals (n = 5) in which Cs+ (2 mm) was used as an Ih blocker. The preferred vertebrate Ih blocker, ZD 7288, does not block leech Ih (multiple unpublished observations). In this set of experiments, ganglia were prepared in the same manner and the same protocol was used; but before the addition of each myomodulin dose, Cs+ was applied for 5-10 min. When analyzing experimental data, burst period was computed as the time between median spikes of adjacent bursts, BD was computed as the time between the first and last spike of each burst, and IBI was computed as the time between the last spike of a burst and the first spike of the subsequent burst.

Mapping experimental data

For each case from a set of experimental conditions (control, Cs+, 1 μm myomodulin, 10 μm myomodulin, 100 μm myomodulin), we determined the parameter pair from within the boundaries of the functional regime which best produced model activity matching the averages of the burst period and the BD of the experimental case. For this, we introduced a simple measure of distance between the model and experimental bursting activity Dist(g¯h,IPumpMax) as follows: Dist(g¯h,IPumpMax)=(BP¯HN−BP¯Exp)2 + (BD¯HN,R−BD¯Exp)22(28) where BP¯HN and BD¯HN,R are the average burst period of both cells and the average BD of the R cell in the HCO model at a given parameter set, respectively; and BP¯Exp and BD¯Exp are the corresponding average values from the experimental condition. We performed a search of model parameters that minimize Dist(g¯h,IPumpMax) to be used as representatives of each experimental condition.

Once a parameter pair was selected for the Cs+ condition, we restricted the choice of parameter pairs for the Cs+ with myomodulin conditions (Cs+ + 1 μm myomodulin, Cs+ + 10 μm myomodulin, Cs+ + 100 μm myomodulin) to the same value of g¯h chosen for the Cs+ only condition.

Curve fitting

Curve fitting of (g¯h,IPumpMax) parameter pairs was accomplished using the MATLAB fit() function from the Curve Fitting Toolbox. A custom inverse function prototype was defined using the fittype() function as follows: IPumpMax(g¯h)=c1 + c2(g¯h−c3)(29)

Initial values for the constants [c1, c2, c3] = [0.47, 0.2, −0.1].

Code accessibility

C simulation files with compilation instructions are available on ModelDB.

Results

Myomodulin decreases the burst period in leech heartbeat HCOs by increasing Ih and reducing IPump (Tobin and Calabrese, 2005). Here, we investigated advantages that such comodulation might provide by considering a model of a leech heartbeat HCO (Fig. 1). The living HN HCOs exhibits alternating, symmetric bursting with a period of 6-9 s. During bursting, the membrane potentials oscillate between −65 mV at the lowest point between bursts, −40 mV baseline potential during a burst, and up to 10 mV at the peaks of action potentials during the burst. Ih and persistent Na+ current, IP, drive the cell toward a certain membrane potential where the two low-threshold Ca2+ currents, ICaS and ICaF, activate and initiate a burst. Then, IP and ICaS support the burst. IPump is always outward, has no reversal potential, and [Na]i determines its magnitude. IPump and IP interact throughout the burst and between bursts (Fig. 2).

Figure 1.
Figure 1.

Leech HN HCO produces bursting activity alternating between two heart interneurons (HNs). A, Schematic representation of the HCO, which is comprised of two neurons located on the left and right sides of a ganglion. Two neurons (large circles) have mutually inhibitory synaptic connections (small circles). Model neurons and synaptic connections are identical. Neurons exhibit bursting activity of membrane potentials (B), which reproduces temporal properties of an extracellular recording of the living heartbeat HCO under control conditions (C).

Figure 2.
Figure 2.

HCO model activity of HN(R) and HN(L) cells represented by membrane potentials, select currents, and intracellular sodium concentration. The cells burst in antiphase (Vm,R and Vm,L) mutually inhibiting one another. Between bursts, Ih, which is activated by hyperpolarization, and persistent sodium current (IP), which is deactivated but still notable, work together to promote the transition to the depolarized spiking phase of bursting. The pump current (IPump) contributes to hyperpolarization and is most active during bursts, where it opposes IP. It is also quite active during the hyperpolarized phase of bursting where it opposes Ih and IP. The bursts are sustained by activated IP. IP is a major contributor to intracellular sodium concentration ([Na]i) during both phases of bursting. The pump current (IPump) rises quickly with [Na]i during the spiking phase of bursting and follows [Na]i between bursts.

Two-parameter map exhibits several activity regimes

To investigate the effects of comodulation of h- and pump currents, we explored model HCO activity in a large domain of two parameters, g¯h and IPumpMax (see Materials and Methods). A two-dimensional parameter sweep revealed a strong influence of these currents on the model activity (Fig. 3). We found that average period of both cells generally decreased as g¯h was increased (Fig. 3A). Coefficient of variation of cycle period was generally low (Fig. 3B); but at low levels of IPumpMax, there were some domains of higher variability. The entire tested parameter space supported oscillatory activity, which had periods ranging between 5 and 12 s. However, not all observed activity regimes were functional. We classified them into three major categories: functional bursting, asymmetric bursting, and plateau-containing regimes (see Materials and Methods) and mapped them along with color-coded temporal characteristics.

Figure 3.
Figure 3.

Two-parameter maps of HCO model activity characteristics obtained under variation of parameters g¯h and IPumpMax. A, Average period of both cells from the HCO model. HCO period ranges from 5 to 12 s where warmer colors represent shorter periods. The general trend suggests that g¯h is an effective parameter for controlling period of activity. B, Coefficient of variation of cycle period, where warmer colors represent higher variability. C, A measure of the asymmetry of depolarized phase durations between the two cells in HCO. Warmer colors represent more asymmetric activity. Asymmetric activity appeared in conditions where g¯h and IPumpMax were either both high or both low. D, Average duration of the depolarized phase of HN(R) activity (D1) and HN(L) activity (D2). Warmer colors represent shorter depolarized phase durations. Distinct complementary “checkerboard” patterns appear in these maps for larger values of both g¯h and IPumpMax indicating asymmetry and, thus, bistability. Coefficient of variation of depolarized phase duration in HN(R) model activity (E1) and HN(L) model activity (E2). F, Ratio of depolarized phases in model cell HN(R) activity (F1) and HN(L) activity (F2) at a given parameter set, which were classified as plateau-like rather than functional bursts. Warmer colors represent traces with more plateau-containing depolarized phases. B, C, E, F, White dashed lines indicate the borders of the functional regime.

Comparison of the depolarized phase duration of each cell in the model HCO revealed a “checkerboard” asymmetric pattern at higher values of g¯h and IPumpMax (Fig. 3D1,D2). In this region, cycle period changes smoothly under variation of the parameters, but the steady-state BDs of the cells are notably different and their proportions in the period flip back and forth between longer and shorter from point to point on the map. Since in our model HCO, both neurons and their connections are identical, it follows that either cell could be the one with longer period. An asymmetric bursting regime with longer BD of one cell must be accompanied by another coexisting regime in which the other cell has a longer BD. The apparent randomness of the dominating side on the map, visualized as complementary checkerboards in depolarized phase duration, is likely a result of high sensitivity of the asymmetric regimes to the choice of the initial states of the neurons. We had attempted to trace regimes during a sweep; each time we tested a new pair of parameters, we set the initial conditions to those of the last time point from the previous simulation. This last point belonged to the steady-state bursting activity attained from the previous simulation, which is slightly different from the steady state which the HCO with the new parameter set will produce and thus the trace is slightly perturbed versus the old steady-state bursting activity. This difference in initial states explains the presence of both steady-state asymmetric regimes in our data. On the map (Fig. 3C), asymmetry measures (see Materials and Methods) ranged from 0 to 1.85, although 80% of evaluated cases throughout the entire map were <1. At higher values of g¯h and IPumpMax, asymmetry was pronounced. In order to operationally label regimes as asymmetric, we chose a conservative asymmetry threshold of 0.2; this represented the case in which the duration of the depolarized phase of one cell in the HCO was ≥55% of the period. In the model, the period was always equal to the sum of the two BDs of both HCO neurons. Although we considered asymmetric regimes as dysfunctional, we suggest that, under relatively low measures of asymmetry, the bursting asymmetric regime would present only a mild risk for animal viability.

Another dysfunctional activity regime was noted by inspection of the coefficient of variation of the cycle period of model activity, which detected an area of high period variability at low values of IPumpMax. Further investigation of patterns in this area on the map revealed a large domain supporting complex activity consisting of plateau events, single spikes, and bursts. While the exact patterns varied widely, they all exhibited depolarized plateau events in one or both cells' activity. Plateaus in this system in vivo would be considered dysfunctional activities, so we characterized these regimes by the proportion of plateaus events among the depolarized phases of the pattern period (see Materials and Methods) (Fig. 3F1,F2). The major predictor of the presence of plateaus was low IPumpMax, as no plateaus were identified in simulations in which IPumpMax was >0.41 nA (Fig. 3F1,F2). Regimes in which either cell exhibited any plateaus were considered dysfunctional and labeled as plateau-containing. Those regimes at low values of IPumpMax, which were both asymmetric and contained plateaus, we generally marked as plateau-containing as this property represents a more serious pathologic condition. These plateaus appear to be similar to seizure-like plateau activity induced by Co2+ or K+ blockers in leech neurons (Angstadt and Friesen, 1991; Opdyke and Calabrese, 1994).

Functional bursting was revealed within a domain of functional regimes on the two-parameter map (Fig. 4A). Its borders are marked in white in Figure 3B, C, E, F. It exhibited a range of periods from 4.9 to 11.8 s. Representative voltage traces from these regimes are shown in Figure 4B.

Figure 4.
Figure 4.

A, Functional and nonfunctional regimes mapped under variation of parameters g¯h and IPumpMax. B, Representative examples of each regime. A, The regimes have been categorized as asymmetric bursting (green), functional bursting (blue), or plateau-containing activity (red). B, Examples of membrane potential traces representing the asymmetric regime (B1) at g¯h = 3.6 nS, IPumpMax = 0.46 nA, the functional regime (B2) at g¯h = 3.6 nS, IPumpMax = 0.4 nA, and the plateau-containing regime (B3) at g¯h = 3.6 nS, IPumpMax = 0.36 nA. This map unveils a functional pathway through parameter space.

Supporting functional bursting through comodulation

By applying a mask of the area labeled functional to the map of burst periods, we obtained our primary finding: comodulation of the two currents can expand the temporal burst characteristics while retaining the functional pattern. Starting from anywhere in this functional comodulation pathway, shaped like a turning stripe, it is clear that by variation of one parameter along either the g¯h or the IPumpMax axes, burst period can only be adjusted in a small range before encountering a transition to a dysfunctional regime (Fig. 4A). By modulating both pump and Ih together, however, bursting cycle period of the HCO can be smoothly controlled over a wide range (Fig. 5A).

Figure 5.
Figure 5.

A color map of burst period within the area supporting the functional regime (A) demonstrates that comodulation of g¯h and IPumpMax expands the range of functional bursting. Overlaid in black is a curve fit of all points in the functional parameter space. This curve represents a general path of smooth comodulation introduced by the application of myomodulin, and it has the following form: IPumpMax(g¯h)=0.3 + 0.85(g¯h + 5.7). Arrows indicate the direction traveled by the HCO model through parameter space to represent increased modulation by myomodulin. B1, Minimum (blue) and maximum (black) achievable period by modulating only g¯hat all test values of IPumpMax. B2, Maximal possible change in burst period modulating only g¯h (black) is less than the achievable change in burst period through comodulation (single value marked in red). C1, Minimum (blue) and maximum (black) achievable period by modulating only IPumpMax at all test values of g¯h. C2, Maximal possible change in burst period modulating only IPumpMax (black) is less than the achievable change in burst period through comodulation (single value marked in red). Global maxima in possible shifts in parameter, burst period, and BD within the functional regime are included in Table 6.

View this table:
Table 2.

Parameters for special currents ILeak and IPump

We searched the functional area of the map and determined the largest viable parameter difference along the g¯h or the IPumpMax axes. Then, we compared burst characteristics at those boundary parameter values to the maximal difference of parameter and burst characteristics projected along a fitted line of comodulation through the entire functional pathway (see Materials and Methods). By searching along the IPumpMax axis, we found that the greatest continuous range of g¯h values occurred at IPumpMax = 0.385 nA, where g¯h could be modulated in the range [3.2, 10] nS to produce a range of burst periods between 7.1 and 4.9 s, a difference of 2.2 s. The greatest achievable range in burst period by modulation of g¯h alone through the functional region occurred at IPumpMax = 0.425 nA where g¯h could be modulated in the range [0, 2.2] nS to produce a range of burst periods between 7.9 and 11.0 s. So, by only modulating g¯h, a range of burst period of 3.1 s is achievable (Fig. 5B). Searching along the g¯h axis, we found that the greatest continuous range of IPumpMax values occurred at g¯h = 0 nS, where IPumpMax could be modulated in the range [0.368, 0.500] nA. These two boundary values of IPumpMax produce burst periods of 8.2 and 11.7 s respectively, a difference of 3.5 s. However, the greatest achievable range in burst periods within this parameter variation occurred at g¯h = 0 nS between IPumpMax = 0.398 nA and IPumpMax = 0.444 nA where the burst periods were 7.0 and 11.8 s, respectively (Fig. 5C). So, by changing only IPumpMax, a change in burst period of 4.8 s is achievable.

To infer a comodulation mechanism coordinating the changes of the pump and Ih, we first hypothesized that comodulation would be most robust to perturbation if it traversed the center of the functional regime area, navigating away from the borders, and that the simplest comodulation mechanism would be the basic reciprocal relation of the two parameters. We fitted a simple curve representing reciprocal relations between g¯h and IPumpMax to all points (i.e., parameter pairs) in the functional parameter area on the map (R2 = 0.7) as follows: IPumpMax(g¯h)=0.3 + 0.85(g¯h + 5.7)(30)

Maximal parameter distance along this curve was between points in parameter space (g¯h, IPumpMax) = (0.0, 0.446) and (g¯h, IPumpMax) = (10, 0.354). Between these two points in parameter space along the curve of comodulation, we found that burst period could be varied between 11.8 and 5.6 s for a maximum achievable range of 6.2 s. This occurred between the parameter points (g¯h, IPumpMax) = (0.0, 0.446) and (g¯h, IPumpMax) = (9.8, 0.355). These results demonstrate that the burst period could be modified over a larger range under comodulation (6.2 s) than under single-parameter variation (4.8 s). This mapping shows that comodulation is a viable mechanism for expanding the range of achievable temporal characteristics. Maximal distances between parameter pairs and maximal ranges in burst statistics for each search condition are provided in Table 5.

Figure 6.
Figure 6.

Mapping experimental results with myomodulin onto the model parameter plane suggest a specific comodulation curve. A, Statistical comparisons of myomodulin dose–response of experimentally measured burst period and BD are provided. In addition, burst period and BD means of individual animals (n = 5) are included (gray). Significant differences in sample means were found between control and all myomodulin conditions (*p < 0.05) as well as between 1 μm myomodulin and 100 μm myomodulin (#p < 0.05), suggesting a dose effect. Error bars on the experimental series indicate SEM. Model means for burst period and BD for best representative cases of each experimental condition are included for comparison. In experiments, asymmetry measures (Table 7), coefficient of variation of burst period (Table 9), and coefficient of variation of BD (Table 11) were low. Each experimental condition was mapped (B) to a single point on the diagram according to burst period and BD associated with a particular dose of myomodulin: control (circle), 1 μm myomodulin (inverted triangle), 10 μm myomodulin (diamond), and 100 μm myomodulin (upright triangle). The relationship between IPumpMax and g¯his then fitted to the positions in parameter space chosen for control and myomodulin conditions producing the following curve:IPumpMax(g¯h)=0.36 + 0.16(g¯h + 0.85). Chosen parameter sets in the model match well with experimental data. C, Model voltage traces are compared directly with extracellular recordings of bilateral HN cells from corresponding conditions.

View this table:
Table 3.

Equations used for calculating the steady-state activation mx∞(V) and inactivation hx∞(V), and time constants of activation τmx(V)(s)and inactivation τhx(V) (s)

Model validation and experimental results

In concert, we performed extracellular recordings of pairs of HN(3) or HN(4) neurons, treated equivalently as there have been no observed consistent differences in these HCOs (Kueh et al., 2016), under control conditions and for three different doses of myomodulin (1, 10, 100 μm) (Fig. 6A). Across experiments (n = 5), myomodulin concentration accounted for the variance in burst period (p < 0.0001) and BD (p < 0.0001). The general trend supported our previous studies, which had found that the addition of 1 μm myomodulin decreased burst period from 12.1 ± 1.2 s to 6.8 ± 1.0 s (Tobin and Calabrese, 2005). These new experiments also showed that myomodulin affects the burst period in a dose-dependent manner (control: 9.0 ± 0.7 s; myomodulin 1 μm: 6.4 ± 0.3 s; myomodulin 10 μm: 5.6 ± 0.3 s; myomodulin 100 μm: 4.2 ± 0.1 s) (Fig. 6A1). BD is also impacted in a dose-dependent manner (control: 5.3 ± 0.4 s; myomodulin 1 μm: 3.8 ± 0.2 s; myomodulin 10 μm: 3.3 ± 0.2 s; myomodulin 100 μm: 2.5 ± 0.1 s) (Fig. 6A2). Significant differences in burst period and BD were found between control conditions and each myomodulin condition (* in Fig. 6A1,A2; p < 0.05 for each), and between the 1 μm and 100 μm conditions (# in Fig. 6A1,A2; p < 0.05). Higher doses of myomodulin caused a greater decrease in period.

A second set of experiments was performed with Cs+ (2 mm), an Ih blocker, to evaluate whether the model matches the action of myomodulin on the pump current alone. Extracellular recordings were performed on HN(3) or HN(4) neurons under control conditions, Cs+ conditions, and the same three doses of myomodulin now with Cs+ introduced before application of each dose of myomodulin (Fig. 7A). Here, the period is also affected in a dose-dependent manner (control: 8.4 ± 0.3 s; Cs+: 9.6 ± 0.6 s; myomodulin 1 μm + Cs+: 7.8 ± 0.5 s; myomodulin 10 μm + Cs+: 6.6 ± 0.3 s; myomodulin 100 μm + Cs+: 4.4 ± 0.2 s) (Fig. 7A1). BD is affected in a dose-dependent manner as well (control: 4.8 ± 0.2 s; Cs+: 4.9 ± 0.3 s; myomodulin 1 μm + Cs+: 3.8 ± 0.3 s; myomodulin 10 μm + Cs+: 2.9 ± 0.2 s; myomodulin 100 μm + Cs+: 1.9 ± 0.1 s) (Fig. 7A2). Across five experiments, we found that the variances in burst period and BD were accounted for by treatment effects (p < 0.0001 and p < 0.0001, respectively). Significant differences in burst period and BD were found between the control condition and the myomodulin 10 μm + Cs+ and myomodulin 100 μm + Cs+ conditions (* in Fig. 7A1,A2; p < 0.05 for each). All combined treatment conditions demonstrated significant differences in burst period and BD from the Cs+ alone condition († in Fig. 7A1,A2; p < 0.05 for each). Burst period was significantly different between all combined treatment conditions (# in Fig. 7A1; p < 0.05 for each), but BD was only significant between myomodulin 10 μm + Cs+ and myomodulin 100 μm + Cs+ conditions (# in Fig. 7A2; p < 0.05 for each).

Figure 7.
Figure 7.

The model matches results from experiments with 2 mm Cs+ (blocks Ih) and myomodulin. A, Statistical comparisons of Cs+ with myomodulin dose–response of experimentally measured burst period (A1) and BD (A2) are provided. In addition, burst period and BD means of individual animals (n = 5) are included (gray). Significant differences in sample means of BD and burst period were found between Cs+-only conditions and all Cs+ with myomodulin conditions (†p < 0.05) and also between control conditions and the Cs+ with 1 μm myomodulin and Cs+ with 100 μm myomodulin (*p < 0.05) conditions. Significant differences in burst period and BD between Cs+ with myomodulin conditions are labeled (#p < 0.05). Error bars on the experimental series indicate SEM. Model means for burst period and BD for best representative cases of each experimental condition are included for comparison. In experiments, asymmetry measures (Table 8), coefficient of variation of burst period (Table 10), and coefficient of variation of BD (Table 12) were low. Each experimental condition is mapped (B) to a single point in the map according to burst period and BD associated with control and application of Cs+ with a particular dose of myomodulin: control (circle), Cs+ (star), Cs+ with 1 μm myomodulin (inverted triangle), Cs+ with 10 μm myomodulin (diamond), and Cs+ with 100 μm myomodulin (upright triangle). The 10 μm and 100 μm myomodulin conditions mapped to the same location. C, Model voltage traces are compared with extracellular recordings of contralateral HN cells from corresponding conditions. HN(R) is the top and HN(L) the bottom trace for both the model and the living neurons.

View this table:
Table 4.

Initial conditions for the 40-dimensional HN HCO modela

View this table:
Table 5.

Maximal continuous parameter and burst statistic shifts along each axis compared with maximal parameter and burst statistic shifts along a curve fitted to the midline of the functional regime in model parameter space (Fig. 5A) and along a curve fitted to the experimental (myomodulin) points (Fig. 6A)a

View this table:
Table 6.

Points in parameter space selected to represent experimental conditions with burst period and BD at those points in parameter space

View this table:
Table 7.

Asymmetry measures from myomodulin only experiments: Experimental Set 1

To verify our model, for each experimental condition, we identified a single “best fit” point within the functional area of the map (Figs. 6B, 7B) at which differences between model values of BD and burst period and corresponding experimental values were minimal (Eq. 27; Table 6). These selected points represent the neuromodulation induced by the different doses of myomodulin; experimental conditions with higher concentrations of myomodulin map to model parameter domains at higher g¯h and lower IPumpMax (Fig. 6B), consistent with the findings of Tobin and Calabrese (2005). The point of best fit between model and Cs+ conditions occurred at g¯h = 1 nS and IPumpMax = 0.448 nA. This point is near the control conditions, but with a lower g¯h value, which is consistent with the blocked Ih experimental condition. The remaining conditions with myomodulin and Cs+ were then mapped to points in the functional area of the map with the same g¯h values under the presumption that in this model g¯h = 1 nS represents an Ih blocked by Cs+ (Fig. 7B). There was good agreement between model and experiment for the Cs+ with 1 μm myomodulin and Cs+ with 10 μm myomodulin conditions, while the Cs+ with 100 μm myomodulin was mapped to the same point in parameter space as the Cs+ with 10 μm myomodulin.

In conclusion, model BD and burst period have good agreement for experimental values under all conditions except Cs+ with 100 μm myomodulin. As the dose of myomodulin increases, burst period and BD decrease and correspond to higher g¯h values and lower IPumpMax values on the map. Direct comparisons between model and experimental voltage traces for both sides of the HCO are provided in Figures 6C and 7C.

To test whether reciprocal relation of the two currents could describe experimental data, we then fit a comodulation curve over to the selected parameter value representations of experimental activity, and suggest that the reciprocal relationship may describe the comodulation of Ih and IPump by myomodulin (R2 = 0.97) (Fig. 6B) as follows: IPumpMax(g¯h)=0.36+0.16(g¯h+0.85)(31)

This experimentally justified comodulation curve has an even greater range of achievable periods than the midline comodulation curve. Maximal parameter distance here occurred between the boundary points (g¯h, IPumpMax) = (0.4, 0.491) and (g¯h, IPumpMax) = (10, 0.38) where the burst period was 11.7 and 4.9 s, respectively. In this case, burst period could be changed by up 6.8 s compared with the midline comodulation curve, which could achieve a 6.2 s range; the modulation of only IPumpMax, which could achieve a 4.8 s range; and the modulation of g¯h alone, which could achieve a 3.1 s range. Remarkably, this comodulation curve and the experimental points were not in the center of the functional regime, but close to the borders of the asymmetric regime. In retrospect, this made sense, as a perturbation resulting in mild asymmetry would not be nearly as catastrophic as a perturbation resulting in a depolarization block.

HN HCOs do not typically burst in asymmetric fashion (Kueh et al., 2016; Wenning et al., 2018). When we analyzed asymmetry in the experimental dataset, we saw that, across the 5 subjects in the myomodulin experiments and also across the 5 subjects in the combined myomodulin with Cs+ experiments, average levels of asymmetry, computed in the same manner as the model data, were all <0.2, our selected cutoff for model asymmetry (Tables 7 and 8). The largest average level of asymmetry was observed in the Cs+ alone condition with an average asymmetry measure of 0.105 and an SE of 0.03 across 5 subjects.

View this table:
Table 8.

Asymmetry measures from combined treatment experiments: Experimental Set 2

The variability in the burst period and BD in experiments was relatively small. Means and SEs for individual leeches are shown alongside grand means and model comparisons in Figures 6A and 7A. Coefficients of variation of burst period averaged across animals did not exceed 0.1 in either the myomodulin only experiments or the Cs+ plus myomodulin experiments (Tables 9 and 10). One subject under control conditions had a coefficient of variation of burst period of 0.128. Across the 5 subjects in the Cs+ plus myomodulin experiments, variability in the coefficient of variation of burst period between treatment conditions was significant (p < 0.05), but average coefficients of variation were <0.1 in all cases. Coefficients of variation of BD averaged across animals were slightly higher but did not exceed 0.13 in either the Cs+ plus myomodulin experiments or the myomodulin only experiments (Tables 11 and 12). The highest coefficient of variation of BD observed in any one preparation was 0.24 in the 100 μm myomodulin with Cs+ condition. The effect of the treatment could not account for the variation in the mean coefficient of variation of BD (p > 0.05). For a full analysis of the burst characteristics of the HN CPG from an extensive dataset, see Wenning et al. (2018).

View this table:
Table 9.

Coefficient of variation of burst period from myomodulin only experiments: Experimental Set 1

View this table:
Table 10.

Coefficient of variation of burst period from combined treatment experiments: Experimental Set 2

View this table:
Table 11.

Coefficient of variation of BD from myomodulin only experiments: Experimental Set 1

View this table:
Table 12.

Coefficient of variation of BD from combined treatment experiments: Experimental Set 2

Discussion

Neuromodulation allows neuronal networks to flexibly adapt to changing behavioral goals and variable environments. An important question posed by experiments with CPGs is why neuromodulators target more than one membrane current. Also intriguing: why is IPump among targeted currents? The leech heartbeat CPG is an ideal system for exploring these questions, offering well-described cellular and circuit dynamics and a biophysical model rooted in experimental data. We investigated activities of a biophysical model of the leech heartbeat HN HCO under parameter variation of the maximal conductance of the Ih (g¯h) and the maximal pump current (IPumpMax). This model of comodulation in leech HNs revealed a mechanism that allows an HCO to regulate its period in a range that is unattainable through modulation of a single current. The two-parameter map of activities revealed a stripe of parameters supporting functional bursting. It is located between large areas of two dysfunctional patterns: asymmetric bursting and plateau-containing seizure-like activity. Within the functional range, we found points that readily matched results from myomodulin dose-dependence experiments. The curve of modulation fitted to these experimental data tracks very close to the border of the asymmetric regime, and far from the more pathologic plateau-containing regime. The mechanism of comodulation of the pump and Ih, identified as their reciprocal relation, expands the functional ranges of g¯h and IPumpMax by navigating between the areas producing dysfunctional patterns. The comodulation offers more than simple synergistic augmentation of the cycle period range, it expands the range of functional bursting periods by 42% compared with modulation of IPumpMax alone or 119% compared with modulation of g¯h alone. Thus, investigating the modulatory space of models can reveal complex interactions that are difficult to isolate in experiments. This comodulation mechanism allows the HCO circuit to produce functional bursting, while regulating its temporal properties in a wide range.

The pump current contributes directly to the dynamics of rhythm generation

The role of the Na+/K+ pump in the dynamics of rhythm generating circuits is not well understood. By executing a vital housekeeping function of maintaining Na+ and K+ gradients across the membrane, it generates an outward current which must be accounted for in the checks and balances of oscillatory pattern generating dynamics. Multiple studies, computational and experimental, have shown how its electrogenic nature is fundamental to the excitability of neurons (Forrest et al., 2012; Krishnan et al., 2015; Kueh et al., 2016; Picton et al., 2017). Here, we show that the pump's contribution is critical for robust yet flexible control of bursting. Significantly, in this HN HCO circuit, modifications to IPump have immediate consequences to the bursting rhythm. Blockage or large reductions in pump activity affect Na+ and K+ ionic gradients with consequent effects on bursting rhythms over the time scale of several bursts. For example, blocking the pump with strophanthidin gradually disrupts the HN HCO's rhythm as the cells reach a depolarization block and are unable to fire action potentials (Kueh et al., 2016). Furthermore, in HN HCOs, IPump is a large, active component contributing to bursting dynamics through interaction with Na+ and other currents, and its downmodulation can be used as an effective mechanism to regulate temporal burst characteristics. The pump current interacts with Ih during the IBI both activating it by causing hyperpolarization and opposing it by being an outward current. Imbalance between these currents leads to plateau-containing, highly variable activity or asymmetric activity, which would be maladaptive for the animal. IPump is thus an active contributor to the dynamics of rhythm generation on both long and short timescales. The ability of IPump to act on short timescales supports its potential role in information processing as well (Forrest, 2014).

Variability and relationship to homeostatic mechanisms

With the exception of the asymmetric cases, the duty cycle of the model is very close to 50%, while experimental traces typically display a duty cycle of ∼60% in control and myomodulin alone. This has been a consistent discrepancy with the HN HCO model since its inception (Hill et al., 2001). The onset of spiking in the model is abrupt at the end of inhibition because of the single-compartment nature of the model, a choice that enhances dynamical analysis but sacrifices the nuances of spiking activity. Despite this issue, these HN HCO models have had predictive value and reproduced the core properties of burst generation in a variety of pharmacological, dynamic clamp, and hybrid-system experiments (Hill et al., 2001; Cymbalyuk et al., 2002; Sorensen et al., 2004; Olypher et al., 2006; García et al., 2008; Weaver et al., 2010; Kueh et al., 2016). We note that, in our experiments, the mean duty cycle across animals decreased to 50% in Cs+ and as low as 40% in Cs+ with 100 μm myomodulin.

We found higher levels of animal-to-animal variability in the myomodulin and Cs+ experiments, suggesting that the modulatory effects on the pump alone may vary between individual animals. We propose that Na+/K+ pump expression levels may be variable between animals, but that at least some other channel expressions, including the Ih, may be coregulated with pump to preserve the described mechanism of comodulation. It is well established that ion channel expression is homeostatically coregulated to bound the types of produced activity (Khorkova and Golowasch, 2007; Marder, 2011; O'Leary et al., 2013; Golowasch, 2019; Goaillard and Marder, 2021). Such homeostatic mechanisms prevent significant deviations from a functional pattern through compensation, a negative feedback loop. While the mechanisms that regulate the proportion of ion channels in the membrane may be in place to bound activity, comodulation of the currents that are already present is a separate problem. We emphasize with this system that comodulation is an active measure, which controls the activity in a wide range while simultaneously avoiding catastrophic dysfunction. From a modeling perspective, the methodology for investigating comodulation and coregulation is similar, but these processes serve different biological functions. While coregulation preserves the activity produced in the face of perturbation and long-term changes, comodulation changes or scales the activity, in this case the period.

Proximity to less dangerous boundary increases robustness

Model investigation of the modulation mechanisms should consider the multiparameter space of modulatory action (Goldman et al., 2001). Here, by covarying two currents implicated in comodulation by myomodulin, we found that the experimentally verified comodulation curve tracks very closely to the edge of our defined functional regime (Fig. 6B), which borders the asymmetric regimes. The transition at this edge from functional activity to asymmetric activity is a smooth one, not a catastrophic transition, like the border between functional and plateau activity; and this positioning may contribute robustness to the HCO. Robust functional bursting in this circuit is critical for the animal's survival, as outputs of HN neurons drive the motor neurons which control the tubular leech hearts. Plateauing seizure-like activity of these cells alone would be deleterious for heart function and could also disrupt ionic gradients, so avoidance of the plateau-containing regime would be critical. This regime is reminiscent of seizure-like plateau activity induced in leech neurons, including HNs by substitution of Ca2+ with Co2+ or by application of K+ channel blockers (Angstadt and Friesen, 1991; Opdyke and Calabrese, 1994). Small variations in symmetry of the activity, however, would have little to no consequences for the animal's welfare. Since two asymmetric regimes coexist, in the living animal, biological noise would likely switch the HCO back and forth between them, resulting in nearly symmetric bursting patterns with BDs equal on average.

Multiple targets of neuromodulation

In some circuits, neuromodulators that produce multiple effects and target different currents in opposite directions may create a mechanism protecting circuits against overmodulation, which would lead a circuit into a dysfunctional regime (Harris-Warrick and Johnson, 2010). Modulators with apparently opposing effects can also underlie complex transient responses because they act on different timescales. Studies on the Aplysia radular neuromuscular system demonstrate that two different classes of modulators are coreleased with acetylcholine (Brezina et al., 2003): myomodulins and small cardiac peptides. These modulators have apparently opposing effects on the muscle and occur on different timescales. Experimentally, myomodulins decrease muscle contractile responses to acetylcholine, whereas small cardiac peptides increase the contractile response of the muscle. The long timescale effects of these modulators on the enhancement of calcium current and muscle relaxation rate allow the system to have a “memory” of the previous behavioral state (activity pattern) when changes occur and increase efficient muscle performance throughout multiple cycles. The fast effects on enhancement of potassium current in the muscle allow the system to adapt when the motor pattern changes with minimal losses in performance. These modulatory systems also allow the muscle contractions to be robust in the face of irregular firing patterns often seen in the behaving animal (Brezina et al., 2005) and to navigate through the parameter space to regions unreachable in the steady state to optimize muscle response throughout different motor patterns. From a neuronal circuit's perspective, use of multiple modulators with multiple cellular targets is functionally similar to our case, where a single modulator affects multiple targets within the cell. We found that, even in the steady state, the effect of modulating multiple ionic currents has tangible benefits for navigating the complex control space of the motor rhythm output. Furthermore, while overmodulation implies that the magnitude of modulation can exceed a functional range and lead a circuit into a dysfunctional regime, in the leech heartbeat HCO (comodulation with similar effects) and in the Aplysia radular neuromuscular system (comodulation with opposing effects), control is not achieved by mitigation of an effect, but rather by the complex combination of effects. These findings enforce that studies on neuromodulation should be using specific, system-relevant combinations of neuromodulator-targeted factors rather than single factors, one at a time. Computational modeling helps bridge this gap by reproducing the emergent properties of a neuronal circuit when physiological studies become less feasible. With their ability to capture these emergent properties from the combination of more basic mechanisms, models can be studied as a dynamical system, which then can be evaluated with simpler, targeted multifactor physiological experiments.

Footnotes

  • The authors declare no competing financial interests.

  • This work was supported by National Institutes of Health Grant 1 R21 NS111355 to G.S.C. and R.L.C.; and Georgia State University Brains and Behavior Fellowship Program to P.J.E. and W.H.B.

  • Correspondence should be addressed to Gennady S. Cymbalyuk at gcymbalyuk{at}gsu.edu

References

  1. ↵
    2. Angstadt JD,
    3. Friesen WO
    (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 pmid:1812221
    OpenUrlCrossRefPubMed
  2. ↵
    2. Arganda S,
    3. Guantes R,
    4. de Polavieja GG
    (2007) Sodium pumps adapt spike bursting to stimulus statistics. Nat Neurosci 10:1467–1473. doi:10.1038/nn1982
    OpenUrlCrossRefPubMed
  3. ↵
    2. Barreto E,
    3. Cressman JR
    (2011) Ion concentration dynamics as a mechanism for neuronal bursting. J Biol Phys 37:361–373. doi:10.1007/s10867-010-9212-6 pmid:22654181
    OpenUrlCrossRefPubMed
  4. ↵
    2. Bertorello A,
    3. Aperia A
    (1990) Short-term regulation of Na(+),K(+)-ATPase activity by dopamine. Am J Hypertens 3:51S–54S. doi:10.1093/ajh/3.6.51s pmid:2166534
    OpenUrlCrossRefPubMed
  5. ↵
    2. Bertorello AM,
    3. Hopfield JF,
    4. Aperia A,
    5. Greengard P
    (1990) Inhibition by dopamine of (Na(+)+K+)ATPase activity in neostriatal neurons through D1 and D2 dopamine receptor synergism. Nature 347:386–388. doi:10.1038/347386a0 pmid:1977083
    OpenUrlCrossRefPubMed
  6. ↵
    2. Brezina V,
    3. Orekhova IV,
    4. Weiss KR
    (2003) Neuromuscular modulation in Aplysia: I. Dynamic model. J Neurophysiol 90:2592–2612. doi:10.1152/jn.01091.2002 pmid:12853443
    OpenUrlCrossRefPubMed
  7. ↵
    2. Brezina V,
    3. Horn CC,
    4. Weiss KR
    (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 pmid:15469963
    OpenUrlCrossRefPubMed
  8. ↵
    2. Calabrese RL
    (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 pmid:9914244
    OpenUrlCrossRefPubMed
  9. ↵
    2. Catarsi S,
    3. Brunelli M
    (1991) Serotonin depresses the after-hyperpolarization through the inhibition of the Na+/K+ electrogenic pump in T sensory neurones of the leech. J Exp Biol 155:261–273. doi:10.1242/jeb.155.1.261
    OpenUrlAbstract/FREE Full Text
  10. ↵
    2. Cymbalyuk GS,
    3. Gaudry Q,
    4. Masino MA,
    5. Calabrese RL
    (2002) Bursting in leech heart interneurons: cell-autonomous and network-based mechanisms. J Neurosci 22:10580–10592. doi:10.1523/JNEUROSCI.22-24-10580.2002 pmid:12486150
    OpenUrlAbstract/FREE Full Text
  11. ↵
    2. Dashevskiy T,
    3. Cymbalyuk G
    (2018) Propensity for bistability of bursting and silence in the leech heart interneuron. Front Comput Neurosci 12:5. doi:10.3389/fncom.2018.00005 pmid:29467641
    OpenUrlCrossRefPubMed
  12. ↵
    2. Doi A,
    3. Ramirez JM
    (2008) Neuromodulation and the orchestration of the respiratory rhythm. Respir Physiol Neurobiol 164:96–104. doi:10.1016/j.resp.2008.06.007 pmid:18602029
    OpenUrlCrossRefPubMed
  13. ↵
    2. Forrest MD
    (2014) The sodium-potassium pump is an information processing element in brain computation. Front Physiol 5:472.
    OpenUrl
  14. ↵
    2. Forrest MD,
    3. Wall MJ,
    4. Press DA,
    5. Feng J
    (2012) The sodium-potassium pump controls the intrinsic firing of the cerebellar Purkinje neuron. PLoS One 7:e51169. doi:10.1371/journal.pone.0051169
    OpenUrlCrossRefPubMed
  15. ↵
    2. Frohlich F,
    3. Bazhenov M,
    4. Timofeev I,
    5. Steriade M,
    6. Sejnowski TJ
    (2006) Slow state transitions of sustained neural oscillations by activity-dependent modulation of intrinsic excitability. J Neurosci 26:6153–6162. doi:10.1523/JNEUROSCI.5509-05.2006 pmid:16763023
    OpenUrlAbstract/FREE Full Text
  16. ↵
    2. García PS,
    3. Wright TM,
    4. Cunningham IR,
    5. Calabrese RL
    (2008) Using a model to assess the role of the spatiotemporal pattern of inhibitory input and intrasegmental electrical coupling in the intersegmental and side-to-side coordination of motor neurons by the leech heartbeat central pattern generator. J Neurophysiol 100:1354–1371. doi:10.1152/jn.90579.2008 pmid:18579654
    OpenUrlCrossRefPubMed
  17. ↵
    2. Goaillard JM,
    3. Marder E
    (2021) Ion channel degeneracy, variability, and covariation in neuron and circuit resilience. Annu Rev Neurosci. Advance online publication. Retrieved Mar 26, 2021. doi: 10.1146/annurev-neuro-092920-121538. doi:10.1146/annurev-neuro-092920-121538
    OpenUrlCrossRef
  18. ↵
    2. Goldman MS,
    3. Golowasch J,
    4. Marder E,
    5. Abbott LF
    (2001) Global structure, robustness, and modulation of neuronal models. J Neurosci 21:5229–5238. pmid:11438598
    OpenUrlAbstract/FREE Full Text
  19. ↵
    2. Golowasch J
    (2019) Neuronal homeostasis: voltage brings it all together. Curr Biol 29:R641–R644. doi:10.1016/j.cub.2019.05.029 pmid:31287985
    OpenUrlCrossRefPubMed
  20. ↵
    2. Grashow R,
    3. Brookings T,
    4. Marder E
    (2009) Reliable neuromodulation from circuits with variable underlying structure. Proc Natl Acad Sci USA 106:11742–11746. doi:10.1073/pnas.0905614106 pmid:19553211
    OpenUrlAbstract/FREE Full Text
  21. ↵
    2. Grashow R,
    3. Brookings T,
    4. Marder E
    (2010) Compensation for variable intrinsic neuronal excitability by circuit-synaptic interactions. J Neurosci 30:9145–9156. doi:10.1523/JNEUROSCI.0980-10.2010 pmid:20610748
    OpenUrlAbstract/FREE Full Text
  22. ↵
    2. Harris-Warrick RM
    (2011) Neuromodulation and flexibility in central pattern generator networks. Curr Opin Neurobiol 21:685–692. doi:10.1016/j.conb.2011.05.011 pmid:21646013
    OpenUrlCrossRefPubMed
  23. ↵
    2. Harris-Warrick RM,
    3. Johnson BR
    (2010) Checks and balances in neuromodulation. Front Behav Neurosci 4:47. doi:10.3389/fnbeh.2010.00047
    OpenUrlCrossRefPubMed
  24. ↵
    2. Harris-Warrick RM,
    3. Johnson BR,
    4. Peck JH,
    5. Kloppenburg P,
    6. Ayali A,
    7. Skarbinski J
    (1998) Distributed effects of dopamine modulation in the crustacean pyloric network. Ann NY Acad Sci 860:155–167. doi:10.1111/j.1749-6632.1998.tb09046.x pmid:9928309
    OpenUrlCrossRefPubMed
  25. ↵
    2. Hazelwood LA,
    3. Free RB,
    4. Cabrera DM,
    5. Skinbjerg M,
    6. Sibley DR
    (2008) Reciprocal modulation of function between the D1 and D2 dopamine receptors and the Na+,K+-ATPase. J Biol Chem 283:36441–36453. doi:10.1074/jbc.M805520200 pmid:18984584
    OpenUrlAbstract/FREE Full Text
  26. ↵
    2. Hill AA,
    3. Lu J,
    4. Masino MA,
    5. Olsen OH,
    6. Calabrese RL
    (2001) A model of a segmental oscillator in the leech heartbeat neuronal network. J Comput Neurosci 10:281–302. doi:10.1023/A:1011216131638
    OpenUrlCrossRefPubMed
  27. ↵
    2. Hodgkin AL,
    3. Huxley AF
    (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 117:500–544. doi:10.1113/jphysiol.1952.sp004764 pmid:12991237
    OpenUrlCrossRefPubMed
  28. ↵
    2. Jasinski PE,
    3. Molkov YI,
    4. Shevtsova NA,
    5. Smith JC,
    6. Rybak IA
    (2013) Sodium and calcium mechanisms of rhythmic bursting in excitatory neural networks of the pre-Botzinger complex: a computational modelling study. Eur J Neurosci 37:212–230. doi:10.1111/ejn.12042 pmid:23121313
    OpenUrlCrossRefPubMed
  29. ↵
    2. Khorkova O,
    3. Golowasch J
    (2007) Neuromodulators, not activity, control coordinated expression of ionic currents. J Neurosci 27:8709–8718. doi:10.1523/JNEUROSCI.1274-07.2007 pmid:17687048
    OpenUrlAbstract/FREE Full Text
  30. ↵
    2. Krishnan GP,
    3. Bazhenov M
    (2011) Ionic dynamics mediate spontaneous termination of seizures and postictal depression state. J Neurosci 31:8870–8882. doi:10.1523/JNEUROSCI.6200-10.2011 pmid:21677171
    OpenUrlAbstract/FREE Full Text
  31. ↵
    2. Krishnan GP,
    3. Filatov G,
    4. Shilnikov A,
    5. 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 pmid:25589588
    OpenUrlCrossRefPubMed
  32. ↵
    2. Kueh D,
    3. Barnett WH,
    4. Cymbalyuk GS,
    5. Calabrese RL
    (2016) Na+/K+ pump interacts with the h-current to control bursting activity in central pattern generator neurons of leeches. Elife 5:e19322. doi:10.7554/Elife.19322
    OpenUrlCrossRefPubMed
  33. ↵
    2. Li YX,
    3. Bertram R,
    4. Rinzel J
    (1996) Modeling N-methyl-D-aspartate-induced bursting in dopamine neurons. Neuroscience 71:397–410. doi:10.1016/0306-4522(95)00483-1 pmid:9053795
    OpenUrlCrossRefPubMed
  34. ↵
    2. Marder E
    (2011) Variability, compensation, and modulation in neurons and circuits. Proc Natl Acad Sci USA 108:15542–15548. doi:10.1073/pnas.1010674108
    OpenUrlAbstract/FREE Full Text
  35. ↵
    2. Marder E
    (2012) Neuromodulation of neuronal circuits: back to the future. Neuron 76:1–11. doi:10.1016/j.neuron.2012.09.010 pmid:23040802
    OpenUrlCrossRefPubMed
  36. ↵
    2. Marder E,
    3. Calabrese RL
    (1996) Principles of rhythmic motor pattern generation. Physiol Rev 76:687–717. doi:10.1152/physrev.1996.76.3.687 pmid:8757786
    OpenUrlCrossRefPubMed
  37. ↵
    2. Marder E,
    3. O'Leary T,
    4. 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 pmid:25032499
    OpenUrlCrossRefPubMed
  38. ↵
    2. Marder E,
    3. Goeritz ML,
    4. Otopalik AG
    (2015) Robust circuit rhythms in small circuits arise from variable circuit components and mechanisms. Curr Opin Neurobiol 31:156–163. doi:10.1016/j.conb.2014.10.012 pmid:25460072
    OpenUrlCrossRefPubMed
  39. ↵
    2. O'Leary T,
    3. Williams AH,
    4. Caplan JS,
    5. Marder E
    (2013) Correlations in ion channel expression emerge from homeostatic tuning rules. Proc Natl Acad Sci USA 110:E2645–E2654. doi:10.1073/pnas.1309966110 pmid:23798391
    OpenUrlAbstract/FREE Full Text
  40. ↵
    2. Opdyke CA,
    3. Calabrese RL
    (1994) A persistent sodium current contributes to oscillatory activity in heart interneurons of the medicinal leech. J Comp Physiol A Neuroethol Sens Neural Behav Physiol 175:781–789. doi:10.1007/BF00191850 pmid:7807419
    OpenUrlCrossRefPubMed
  41. ↵
    2. Olypher A,
    3. Cymbalyuk G,
    4. Calabrese RL
    (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 pmid:16943313
    OpenUrlCrossRefPubMed
  42. ↵
    2. Parker JR,
    3. Bondy B,
    4. Prilutsky BI,
    5. Cymbalyuk GS
    (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 pmid:29766765
    OpenUrlCrossRefPubMed
  43. ↵
    2. Picton LD,
    3. Zhang H,
    4. Sillar KT
    (2017) Sodium pump regulation of locomotor control circuits. J Neurophysiol 118:1070–1081. doi:10.1152/jn.00066.2017 pmid:28539392
    OpenUrlCrossRefPubMed
  44. ↵
    2. Pulver SR,
    3. Griffith LC
    (2010) Spike integration and cellular memory in a rhythmic network from Na+/K+ pump current dynamics. Nat Neurosci 13:53–59.
    OpenUrlCrossRefPubMed
  45. ↵
    2. Seamans JK,
    3. Gorelova N,
    4. Durstewitz D,
    5. Yang CR
    (2001) Bidirectional dopamine modulation of GABAergic inhibition in prefrontal cortical pyramidal neurons. J Neurosci 21:3628–3638. pmid:11331392
    OpenUrlAbstract/FREE Full Text
  46. ↵
    2. Scuri R,
    3. Lombardo P,
    4. Cataldo E,
    5. Ristori C,
    6. Brunelli M
    (2007) Inhibition of Na+/K+ ATPase potentiates synaptic transmission in tactile sensory neurons of the leech. Eur J Neurosci 25:159–167. doi:10.1111/j.1460-9568.2006.05257.x pmid:17241277
    OpenUrlCrossRefPubMed
  47. ↵
    2. Sharples SA,
    3. Whelan PJ
    (2017) Modulation of rhythmic activity in mammalian spinal networks is dependent on excitability state. eNeuro 4:ENEURO.0368-16.2017. doi:10.1523/ENEURO.0368-16.2017
    OpenUrlCrossRef
  48. ↵
    2. Sorensen M,
    3. DeWeerth S,
    4. Cymbalyuk G,
    5. Calabrese RL
    (2004) Using a hybrid neural system to reveal regulation of neuronal network activity by an intrinsic current. J Neurosci 24:5427–5438. doi:10.1523/JNEUROSCI.4449-03.2004 pmid:15190116
    OpenUrlAbstract/FREE Full Text
  49. ↵
    2. Spitzer N,
    3. Cymbalyuk G,
    4. Zhang H,
    5. Edwards DH,
    6. Baro DJ
    (2008) Serotonin transduction cascades mediate variable changes in pyloric network cycle frequency in response to the same modulatory challenge. J Neurophysiol 99:2844–2863. doi:10.1152/jn.00986.2007 pmid:18400960
    OpenUrlCrossRefPubMed
  50. ↵
    2. Tang LS,
    3. Taylor AL,
    4. Rinberg A,
    5. Marder E
    (2012) Robustness of a rhythmic circuit to short- and long-term temperature changes. J Neurosci 32:10075–10085. doi:10.1523/JNEUROSCI.1443-12.2012 pmid:22815521
    OpenUrlAbstract/FREE Full Text
  51. ↵
    2. Teshiba T,
    3. Shamsian A,
    4. Yashar B,
    5. Yeh SR,
    6. Edwards DH,
    7. Krasne FB
    (2001) Dual and opposing modulatory effects of serotonin on crayfish lateral giant escape command neurons. J Neurosci 21:4523–4529. pmid:11404440
    OpenUrlAbstract/FREE Full Text
  52. ↵
    2. Tobin AE,
    3. Calabrese RL
    (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 pmid:16093342
    OpenUrlCrossRefPubMed
  53. ↵
    2. Weaver AL,
    3. Roffman RC,
    4. Norris BJ,
    5. Calabrese RL
    (2010) A role for compromise: synaptic inhibition and electrical coupling interact to control phasing in the leech heartbeat CPG. Front Behav Neurosci 4:38.
    OpenUrlCrossRefPubMed
  54. ↵
    2. Wenning A,
    3. Norris BJ,
    4. Günay C,
    5. Kueh D,
    6. Calabrese RL
    (2018) Output variability across animals and levels in a motor system. Elife 7:e31123. doi:10.7554/Elife.31123
    OpenUrlCrossRef
  55. ↵
    2. Yu N,
    3. Morris CE,
    4. Joos B,
    5. Longtin A
    (2012) Spontaneous excitation patterns computed for axons with injury-like impairments of sodium channels and Na/K pumps. PLoS Comput Biol 8:e1002664. doi:10.1371/journal.pcbi.1002664 pmid:23028273
    OpenUrlCrossRefPubMed
  56. ↵
    2. Zhang HY,
    3. Sillar KT
    (2012) Short-term memory of motor network performance via activity-dependent potentiation of Na+/K+ pump function. Curr Biol 22:526–531. doi:10.1016/j.cub.2012.01.058 pmid:22405867
    OpenUrlCrossRefPubMed
  57. ↵
    2. Zhang HY,
    3. Picton L,
    4. Li WC,
    5. Sillar KT
    (2015) Mechanisms underlying the activity-dependent regulation of locomotor network performance by the Na+ pump. Sci Rep 5:16188. doi:10.1038/srep16188 pmid:26541477
    OpenUrlCrossRefPubMed
  58. ↵
    2. Zhang LN,
    3. Su SW,
    4. Guo F,
    5. Guo HC,
    6. Shi XL,
    7. Li WY,
    8. Liu X,
    9. Wang YL
    (2012) Serotonin-mediated modulation of Na+/K+ pump current in rat hippocampal CA1 pyramidal neurons. BMC Neurosci 13:10. doi:10.1186/1471-2202-13-10 pmid:22257758
    OpenUrlCrossRefPubMed
  59. ↵
    2. Zhang LN,
    3. Li JX,
    4. Hao L,
    5. Sun YJ,
    6. Xie YH,
    7. Wu SM,
    8. Liu L,
    9. Chen XL,
    10. Gao ZB
    (2013) Crosstalk between dopamine receptors and the Na(+)/K(+)-ATPase (review). Mol Med Rep 8:1291–1299. doi:10.3892/mmr.2013.1697 pmid:24065247
    OpenUrlCrossRefPubMed
  60. ↵
    2. Zhang Y,
    3. Bucher D,
    4. Nadim F
    (2017) Ionic mechanisms underlying history-dependence of conduction delay in an unmyelinated axon. Elife 6:e25382. doi:10.7554/Elife.25382
    OpenUrlCrossRef
Back to top

In this issue

Email
Print
View Full Page PDF
Citation Tools
Respond to this article
Comodulation of h- and Na+/K+ Pump Currents Expands the Range of Functional Bursting in a Central Pattern Generator by Navigating between Dysfunctional Regimes
Parker J. Ellingson, William H. Barnett, Daniel Kueh, Alex Vargas, Ronald L. Calabrese, Gennady S. Cymbalyuk
Journal of Neuroscience 28 July 2021, 41 (30) 6468-6483; DOI: 10.1523/JNEUROSCI.0158-21.2021
Twitter logo Facebook logo Mendeley logo

Jump to section

Keywords

Responses to this article

Respond to this article

Jump to comment:

No eLetters have been published for this article.