Bursting Dynamics Based on the Persistent Na+ and Na+/K+ Pump Currents: A Dynamic Clamp Approach

Abstract Life-supporting rhythmic motor functions like heart-beating in invertebrates and breathing in vertebrates require an indefatigable generation of a robust rhythm by specialized oscillatory circuits, central pattern generators (CPGs). These CPGs should be sufficiently flexible to adjust to environmental changes and behavioral goals. Continuous self-sustained operation of bursting neurons requires intracellular Na+ concentration to remain in a functional range and to have checks and balances of the Na+ fluxes met on a cycle-to-cycle basis during bursting. We hypothesize that at a high excitability state, the interaction of the Na+/K+ pump current, Ipump, and persistent Na+ current, INaP, produces a mechanism supporting functional bursting. INaP is a low voltage-activated inward current that initiates and supports the bursting phase. This current does not inactivate and is a significant source of Na+ influx. Ipump is an outward current activated by [Na+]i and is the major source of Na+ efflux. Both currents are active and counteract each other between and during bursts. We apply a combination of electrophysiology, computational modeling, and dynamic clamp to investigate the role of Ipump and INaP in the leech heartbeat CPG interneurons (HN neurons). Applying dynamic clamp to introduce additional Ipump and INaP into the dynamics of living synaptically isolated HN neurons in real time, we show that their joint increase produces transition into a new bursting regime characterized by higher spike frequency and larger amplitude of the membrane potential oscillations. Further increase of Ipump speeds up this rhythm by shortening burst duration (BD) and interburst interval (IBI).


Introduction
Life-sustaining rhythmic motor behaviors such as breathing are controlled by oscillatory neuronal networks known as central pattern generators (CPGs). The oscillatory patterns in CPGs arise from the interaction of endogenous membrane and synaptic currents (Harris-Warrick, 1993;Marder and Calabrese, 1996;Prinz et al., 2004b). CPGs must be flexible and adapt to behavioral goals and environmental conditions by changing the motor output period and spike frequency (Harris-Warrick and Marder, 1991;Katz, 1996). Many studies have demonstrated that membrane and synaptic properties are subject to neuromodulation that implements changes in the rhythmic patterns. Among these properties, neuromodulators such as myomodulin and dopamine target the activity of the Na 1 / K 1 pump (Stimers and Dobretsov, 1998;Dobretsov and Stimers, 2005;Tobin and Calabrese, 2005;Wang et al., 2012;Kueh et al., 2016;Picton et al., 2017b). The Na 1 /K 1 pump sculpts the activity of CPG neurons by generating an outward current reflecting Na 1 entry associated with spiking activity. This feature, taken together with its voltage independence, distinguishes the pump current from typical voltage-gated membrane currents (Dobretsov and Stimers, 2005;Kueh et al., 2016;Picton et al., 2017b;Ellingson et al., 2021;Sharples et al., 2021). The biophysical mechanisms underlying the contribution of pump activity to the generation of the neuronal rhythms are not fully understood.
The persistent Na 1 current (I NaP ) is a noninactivating voltage-gated inward current found in HN neurons; I NaP supports burst formation and controls spike frequency (Opdyke and Calabrese, 1994). I NaP is a target of neuromodulation by neuropeptide FMRFamide and the rhythmicity of HN neurons is highly sensitive to I NaP . For example, low concentrations of FMRFamide speed up the burst rhythm of HN neurons, while high concentrations disrupt neural activity (Simon et al., 1992;Schmidt et al., 1995;Nadim and Calabrese, 1997). I NaP must be balanced by an outward current such as noninactivating K 1 current , leak current [Doloc-Mihu and Calabrese, 2014; e.g., the HN persistent K 1 current (I K2 ); see Materials and Methods], or I pump . However, the dynamics that arise from the interaction between I NaP and either I K2 or I pump are fundamentally different because I K2 represents a typical conductance-based voltage-activated current while I pump is a voltage-independent Na 1 -activated current.
Computational modeling and real-time hybrid systems such as the dynamic clamp are essential tools to unravel neuronal dynamics. We implemented a dynamic clamp to inject I NaP and I pump into single HN neurons based on a published computational model (A.A. Hill et al., 2001;Cymbalyuk et al., 2002;A.A.V. Hill et al., 2002;Olypher et al., 2006;Barnett and Cymbalyuk, 2011;Kueh et al., 2016;Ellingson et al., 2021;Erazo-Toscano et al., 2021). Our methods allow us to investigate how the interaction of these two ionic currents affects the rhythmic activity of isolated HN neurons. Our experiments with the HN neurons revealed a new bursting regime distinguished by a large amplitude of voltage oscillations and high spike frequency (HA bursting), elicited by the interaction of I NaP and I pump . We determined that increasing I pump simultaneously decreased burst duration (BD) and interburst interval (IBI) of this bursting regime. We have developed a reduced computational model that explains the mechanism underlying the speeding up of the neuronal rhythm. We demonstrate that the interaction of I NaP and I pump creates a flexible and robust oscillatory mechanism that supports the HA bursting.

Animals, solutions, and preparation
Leeches Hirudo verbana (https://www.leech.com/ collections/live-leeches) were kept in artificial pond water [0.05% w/v Instant Ocean Sea (Spectrum Brands Inc.) diluted in deionized water] at 16°C on a 12/12 h light/dark cycle. Before each dissection, the animals were immobilized by burying them in a bed of crushed ice. For the dissection, the animal was immersed in a dissecting dish filled with chilled leech saline containing 115 mM NaCl, 4 mM KCl, 1.7 mM CaCl 2 , 10 mM D-glucose, and 10 mM HEPES in deionized water; pH adjusted to 7.4 with 1 M NaOH. The animals were pinned dorsal side up, and a longitudinal incision cut to expose the internal organs. Individual mid-body ganglia 7 were removed from the animal and pinned ventral side up in a 35-mm sylgard-coated (Sylgard, Dow Corning) Petri dish and covered with saline. Directly before the experiment, the glial sheath was removed from the ganglion with microscissors and a scalpel. The preparation was superfused continuously with leech saline (bath volume ;0.5 ml, flow rate ;5 ml/min). All experiments were performed at room temperature (21-22°C).

Electrophysiology
Electrodes were pulled on a Flaming/Brown micropipette puller (P-97, Sutter Instruments) from borosilicate glass (1 mm OD, 0.75 mm ID; A-M Systems). The glass electrodes were filled with 2 M KAcetate and 20 mM KCl; we performed quality control to select only electrodes with a resistance higher than 30 MV. All intracellular traces were acquired using an Axoclamp 2A amplifier (Molecular Devices) in discontinuous current-clamp mode (DCC) at a switching/sampling rate .3 kHz. Intracellular signals were digitized with an Axon Digidata 1440A digitizer (5-kHz acquisition rate) and recorded with the Clampex 10.4 software (Molecular Devices).
After impaling the neuron with the sharp microelectrode, we injected a steady negative current (À0.3 nA) into the neuron for ;3 min, then progressively removed the injected bias current in steps of 0.01 nA until the magnitude of the injected current was reduced to À0.1 nA; the neuron was kept with À0.1-nA steady current for the remainder of the experiment.

Dynamic clamp
HN neurons have been extensively studied using voltage clamp, and the following currents were described: a fast Na 1 current (I NaF ), a persistent Na 1 current (I NaP ; Opdyke and Calabrese, 1994), a hyperpolarization-activated inward current (I h ; Angstadt and Calabrese, 1989), and three K 1 currents (Simon et al., 1992), delayed rectifier-like K 1 current (I K1 ), a persistent K 1 current (I K2 ), and a fast transient K 1 current (I KA ). For each of these currents model, Hodgkin-Huxley style equations have been developed (A.A. Hill et al., 2001;A.A.V. Hill et al., 2002;Kueh et al., 2016;Ellingson et al., 2021;Erazo-Toscano et al., 2021). A pump current has also been isolated in voltage clamp (Tobin and Calabrese, 2005), characterized, and model equations developed (Kueh et al., 2016). The dynamic clamp is implemented in MATLAB and Simulink (MathWorks) on dSpace real-time digital signal processing boards: the DS1104 R&D, and the DS1103 PPC (dSpace Inc.; Olypher et al., 2006;Barnett and Cymbalyuk, 2011). The dynamic clamp system executed the model HN persistent Na 1 current (I NaP ) and the model HN Na 1 /K 1 pump current (I pump ) equations so that these currents can be added to or subtracted from the recorded neuron. The dynamic clamp implementation of these currents has been described previously (Barnett and Cymbalyuk, 2011;Erazo-Toscano et al., 2021). Briefly, the dynamic clamp reads the membrane potential of the neuron from the Axoclamp 2A amplifier and updates the current injected by the amplifier according to model equations evaluated in real time (.3 kHz). The dynamic clamp also estimates the cell's native fast Na 1 current (spike-generating current; I NaF ), native persistent Na 1 current (I NaP native ) and native pump current (I pump native ) along with the injected I NaP and I pump to then estimate the intracellular Na 1 concentration ([Na 1 ] i ) using Equation 1: (1) where vF (vF ¼ 0.0024 nl Â C/mol) is the volume (% 4:4 rl) of the cytosolic Na 1 reservoir multiplied by the Faraday constant (96,485 C/mol). The extracellular Na 1 concentration ([Na 1 ] o ) is assumed to be constant. In the dynamic clamp model, we simplified the evaluation of Na 1 dynamics and did not consider Na 1 fluxes generated by the leak and h-currents (compare to Eq. 4). The calculated intracellular Na 1 ½ i is used to compute the injected I pump and to estimate native pump current: where I pump max , [Na 1 ] ih , [Na 1 ] is determine the maximal value, the intracellular Na 1 concentration for half-activation, and the sensitivity of the Na 1 /K 1 pump current to [Na 1 ] i , respectively. On the other hand, the injected I NaP depends on the recorded HN neuron's membrane potential and the evaluated [Na 1 ] i , which determines Na 1 reversal potential. In our dynamic-clamp experiments, the I pump and I NaP were controlled by I pump max and g NaP , respectively; and where noted, we estimated the native currents I NaP native ( g NaP native ¼ 5 nS) and I pump native (I pump max native ¼ 0.1 nA), which were assumed to be endogenously present in HN neurons.
Dynamic clamp model is described by the following system of ordinary differential equations with the parameters provided in Table 1 with the steady state activation and inactivation functions described by the Boltzmann C m dV dt ¼ ÀðI NaF 1 I NaP 1 I NaP native 1 I K1 1 I K2 1 I KA 1 I h 1 I CaF 1 I CaS 1 I pump 1 I leak Þ

Data analysis of experimental burst characteristics
We computed the experimental burst characteristics using previously described methods . We detected spikes in the membrane potential traces with a Matlab function identifying the peaks and selecting those with a height larger than 20 mV. First, we computed the interspike intervals (ISIs) and the ISIs larger than 800 ms were discriminated as the interburst intervals (IBIs). The burst duration (BD) was the time between the first and last spikes of the burst (Fig. 1), where a burst was defined as a group of five or more spikes. A group with a smaller number of spikes was discarded and its time interval was concatenated with the two large ISIs, surrounding this small group, into a single interburst interval. The spiking activity within bursts was characterized by the total number of spikes, instantaneous spike frequency defined for each spike as 1/ISI of the interval following its peak, and the average intraburst spike frequency which for short we call spike frequency. The interburst interval (IBI) was the time between the last spike of a burst and the first spike potential of the subsequent burst (Fig. 1). The IBI was assigned to the preceding burst. The cycle period (T), used to assess rhythm regularity, was defined as the time difference between the median spikes of two subsequent bursts ( Fig. 1) and was assigned to the first of the two bursts. We required a minimum of eight consecutive bursts that could be detected for a sample recording to be accepted into the analysis.
We characterized a spiking-silence waveform for a bursting cycle by a basic envelope of the membrane potential oscillation with depolarized and hyperpolarized phases. For the depolarized phase, we averaged membrane potential (V m ) over the time interval between the consecutive minima bounding each spike and assigned the value obtained for every point within this interval. Between bursts, V m was harvested as is, without averaging. We used the envelope of the bursting cycle to compare membrane potential waveforms obtained for different parameter sets and to compare with the simple oscillatory waveform of a two-dimensional model describing interaction of the membrane potential and intracellular Na 1 concentration (Fig. 1D).
We computed the amplitudes of the averaged V m , [Na 1 ] i , and I pump oscillations of the bursting envelope by detecting the peaks and troughs of the signal and then subtracting the trough values from the peak values in a pairwise fashion. We also computed the median [Na 1 ] i   Figure 1. Experimental data from a single HN(7) neuron with I NaP and I pump introduced by dynamic clamp to support bursting activity. In A-C, membrane potential V m (mV; top), the instantaneous spike frequency (Hz; middle), and virtual internal Na 1 concentration ([Na 1 ] i ; mM; bottom) estimated by the dynamic clamp model. In this and all subsequent figures [Na 1 ] i is calculated using both estimated native (I pump , I NaP , and I NaF ) and dynamic clamp injected (I pump and I NaP ) currents. Parameter values for g NaP and I pump max are those used for dynamic clamp injected currents. The horizontal bars above the membrane potential traces mark the temporal characteristics of the bursting: bursting cycle period (T), the interburst interval (IBI), and the burst duration (BD). The gray diamonds over each burst tag the median spike used to compute the burst cycle period. The thick colored curves over the voltage traces are the (voltage) envelopes of burst waveforms (Materials and Methods). In this and subsequent figures, high voltage amplitude (HA), high spike frequency bursts are shown in red with the corresponding voltage envelope in yellow, and low voltage amplitude (LA), low spike frequency bursts are shown in blue with the corresponding voltage envelope in yellow as well. A-C, Sample data for different dynamic clamp g NaP and I pump max parameters from the same preparation. A, Coefficient of variation for cycle period: 0.25. B, Coefficient of variation for cycle period: 0.38. C, Coefficient of variation for cycle period: 0.1. D1, D2, Computed voltage envelope with spike minima identified by ovals. These are representative bursts, LA (D1) and HA (D2), respectively, with different dynamic clamp parameters.
Research Article: New Research and I pump , as the median value of the signal between peaks and troughs of the [Na 1 ] i and I pump , respectively; thus, the median was computed for the backside of the waveform (from maximum to minimum), only. We reported the averaged median across recorded cycles of oscillations.
In summary, we gathered the following characteristics of bursting activity: burst duration (BD), interburst interval (IBI), cycle period (T), amplitude of the voltage envelope, instantaneous spike frequency, spike frequency, [Na 1 ] i envelope amplitude, median [Na 1 ] i , I pump envelope amplitude, and median I pump . We performed normalization by dividing the dependent variables (DVs) of interest (BD, IBI, [Na 1 ] i and I pump amplitude values, [Na 1 ] i and I pump median values) by the value of the corresponding DV at I pump max ¼ 0.3 nA, g NaP ¼ 6 nS. These reference values of the parameters were chosen since this combination consistently supported regular rhythmic bursting. Normalization was performed within each individual.

Machine learning algorithms
We implemented a two-threshold algorithm to discriminate between low and high-amplitude oscillations of the membrane potential and [Na 1 ] i in the dynamic clamp experimental data (Fig. 1D1,D2). An expert assisted in discriminating the two regimes with the two-threshold algorithm on 35% of the randomly selected data points by choosing the thresholds. As the result, the classification algorithm assigned a value of "1" to the neuron sample activity if both the membrane potential amplitude and intracellular sodium concentration amplitude were larger than the respective thresholds of 20 mV and 1 mM, indicating high-amplitude bursting activity. Conversely, a value of "0" was assigned if either or both of the amplitudes were below or equal to their respective thresholds, indicating low amplitude bursting activity. The output of the two-threshold algorithm was used to train a Gaussian naive Bayes supervised machine learning (ML) algorithm, which then was applied to classify the remaining 65% of the data. The ML algorithm was implemented in Python 3.6.8 with the Scikit-learn 1.1.3 library obtained from the Anaconda distribution; and further details about the naive Bayes algorithm can be found in the scikitlearn documentation (https://scikit-learn.org/stable/modules/ naive_bayes.html).
We performed linear regression analyses between the amplitude of voltage bursting envelope (voltage amplitude) and [Na 1 ] i envelope amplitude and between the voltage amplitude and spike frequency; the regression analyses were performed in two complementary ways. First, we computed collinearity in an experiment-by-experiment fashion. Second, we computed collinearity on the pooled data. All statistical analyses, including linear regression and nonparametric Mann-Whitney U test were performed using Python 3.6.8 and Scipy 1.8 stats from the Anaconda distribution.

Simplified HN neuron model
We propose that the interaction between I NaP and I pump can explain the dynamical mechanism underlying the phenomena we observed experimentally and have developed a simplified version of the HN model focused on this interaction in the canonical HN model (A.A. Hill et al., 2001;A.A.V. Hill et al., 2002;Kueh et al., 2016;Ellingson et al., 2021;Erazo-Toscano et al., 2021). It has two state variables: membrane potential, V m , and intracellular Na 1 concentration, Na 1 ½ i . The membrane potential variable governs four voltage-gated currents with instantaneous activation/inactivation: persistent sodium current with the two components, I NaP and I NaP native , fast sodium current (I NaF ), hyperpolarization-activated current (I h ), and noninactivating potassium current (I K2 ). The intracellular Na 1 concentration variable governs the pump current (I pump ) and modulates Na 1 currents through the changes in Na 1 reversal potential. Leak and hyperpolarization-activated currents were split into K 1 and Na 1 components (I leak-K and I leak-Na and I h-K and I h-Na ). This model adopted the same parameters values of the two injected currents I NaP and I pump used in experiments. The equations of the simplified model are as follows: Equations 3 and 4 estimate division between Na 1 and K 1 components of h-current as I hÀNa ¼ 3 7 I h and I hÀK ¼ 4 7 I h suggested by the h-current reversal potential À21 mV (Kueh et al., 2016). In contrast to the dynamic clamp, in the 2D model, the evaluation of Na 1 dynamics includes native Na 1 currents in addition to dynamic clamp currents. The individual voltage-gated currents were computed with instantaneous activation and inactivation gating variables: where g i and E i are the maximal conductance and reversal potential, respectively, with the subscript "i" identifying the type of the current ( Table 2). The reversal potential for Na 1 currents E Na was computed using the Nernst The voltage-gated activation and inactivation are instantaneous and denoted as the steady-state functions m 1i ðV m Þ and h 1i ðV m Þ, respectively (Table 3). I pump is determined by Equation 2 with parameters from Table 4. I leak-K and I leak-Na are Ohmic currents (parameters in Table 4), determined by equation where i stands for either K or Na: Data analysis of the 2D model The 2D HN interneuron model produces a simple oscillatory voltage waveform similar to the experimental membrane potential envelope waveform. To divide the cycle period into the phases corresponding to the bursting waveform, we set a threshold at À45 mV and considered the depolarized and hyperpolarized phases as corresponding to the burst and silent phases, respectively following the definition of the phases established previously . Thus, in the 2D model simulated oscillations, BD is the time interval the voltage is above the threshold, and the IBI is the time interval when it is below the threshold. To compute the simulated amplitudes of V m and [Na 1 ] i we subtracted the troughs from the peaks of the signal.

Optimization of the simplified HN model
We used an evolutionary algorithm to find a 2D model parameter set that fits the envelop waveforms and temporal characteristics of bursting activity of an experimental set of used values of I pump max . The evolutionary algorithm optimized the parameters of native currents representing a given recorded living neuron: maximal conductance of native I hg h , the voltage of half-activation (V ½ ) of I h , conductances of K 1 and Na 1 components of leak current, g leak-K and g leak-Na , respectively, maximal conductance, steepness of activation, and voltage of half-activation of the native persistent Na 1 current, and the maximal conductance of persistent K 1 currentg K2 . The algorithm kept the dynamic-clamp currents (injected I NaP and I pump ) parameters intact.
To assess the quality of the model fit, we introduced a cost function that took into account the differences in temporal characteristics and deviations in maxima and minima of membrane potential and [Na 1 ] i waveforms between the experimental and model activities.
The evolutionary algorithm generated a population of an "f" number of parameter sets randomly generated around an initial canonical set. We also refer to these parameter sets as individuals. The outcome of the cost function evaluated all individuals, and the parameter set with the smallest value of the function was passed on to create the next generation. At the next generation, a population of individuals was created by adding small noise (2.5% of the parameter value) to each optimization parameter of the passed individual (Ashlock, 2006;Eiben, 2015). The algorithm iterated the process for "n" number of generations. We optimized the following targeted parameters: g leak-Na , g leak-K , g h , g NaF , native g NaP , the voltage of halfactivation of native I NaP-native , steepness of activation of native I NaP-native , g K2 and the voltage of half-activation of I K2 . These parameters represented cell to cell variability between preparations. The cost function incorporated the outcome of the experimental protocol and kept the experimental values of the parameters g NaP and I pump max the same while varying only the targeted parameters. An example of the code and analysis of a 2D model with a parameter set fitted to the bursting activity of an experimental set could be found in the Extended Data 1.

Results
We investigated the dynamics emerging from the interaction between I NaP and I pump in isolated leech heart interneurons (HNs). We focused on the HN(7) interneurons because, it was essential for this study that in the seventh ganglion the two HNs do not have mutual synaptic connections and are thus isolated as soon as this ganglion is extracted. Thus, this preparation allows us to readily investigate single-cell dynamics without pharmacological intervention (Calabrese, 1977). We introduced the two currents with the dynamic clamp and used two different protocols to investigate the effect of varying g NaP and I pump max on activity of the interneuron. For the first protocol, we kept I pump max constant at a value taken from the range [0.1-0.9 nA] (N ¼ 40) and randomly varied g NaP within the range between 1.0 and 7.0 nS. For the second protocol, we kept g NaP constant at 6 nS (N ¼ 13) and varied I pump max within the range between 0.1 and 0.9 nA in random order. For each step in a protocol, we assessed whether the system had reached steady state by attainment of a [Na 1 ] i Research Article: New Research amplitude that showed small stochastic variability ( Fig.  2A) and analyzed all cycles after this point requiring a minimum of eight for the data to be considered. Some animals were used for both protocols if the recorded interneuron remained functional for a sufficient length of time. In total, 41 animals were used for this study.

Injection of small I NaP and I pump reinstates bursting in HN neurons
The dynamics of a single HN are very sensitive to the shunting leak current introduced by sharp microelectrodes. The isolated HN interneurons exhibit bursting activity, if recorded extracellularly with a loose whole-cell patch pipette, and tonic spiking activity, if recorded intracellularly with a sharp microelectrode . We found that we could reinstate bursting activity in HN neurons recorded with sharp microelectrodes by introducing relatively small currents, I NaP ( g NaP ¼ 1 nS) and I pump (I pump max ¼ 0.1 nA), using dynamic clamp. The reestablished bursting activity has properties typical for HN bursting: low voltage envelope amplitude underlying the slow membrane potential oscillations and the spike frequency similar to that obtained in intracellular recordings from neurons incorporated in a half-center oscillator (HCO; Fig.  1A; Tobin and Calabrese, 2005;Kueh et al., 2016). Below, we call this bursting regime with low amplitude voltage oscillations, LA bursting.
The combination of I pump max and g NaP distinguishes two bursting regimes exhibited by hybrid HN neurons Using the first dynamic clamp protocol, we systematically increased g NaP while keeping I pump max at a constant value (Fig. 1A,B,D2). When we increased g NaP from 1.0 to 3.0 nS, the burst duration and interburst interval became irregular (Fig. 1B, coefficient of variation ¼ 0.38). When we further increased g NaP to 5 nS, the bursting activity became regular again and its waveform changed: the spike frequency and amplitude of the underlying envelope oscillations of V m and [Na 1 ] i became visibly higher (Fig. 1C, D2). These observations suggested the emergence of a new bursting regime as a result of the variation of I pump max and g NaP . Below, we call this regime the high-amplitude (HA) bursting.
We systematically varied dynamic-clamp parameters, I pump max and g NaP ; to test whether these two bursting regimes are quantitatively distinct. We found that the LA bursting possesses the following properties, the amplitude of membrane potential envelope had mean ¼ 4.1 mV, SD ¼ 2.1 mV (n ¼ 15), spike frequency had mean ¼ 4.1 Hz, SD ¼ 1.4 Hz (n ¼ 15), and estimated amplitude of [Na 1 ] i oscillation was relatively low with mean ¼ 0.42 mM, SD ¼ 0.14 mM (n ¼ 15; Figs. 1A, 2A). The other regime, the HA bursting, has high amplitude of envelope of membrane potential with mean ¼ 44.8 mV, SD ¼ 1.16 mV (n ¼ 22), high spike frequency with mean ¼ 17.11 Hz, SD ¼ 0.9 Hz (n ¼ 22), and high amplitude of [Na 1 ] i oscillation with mean ¼ 2 mM, SD ¼ 0.13 mM (n ¼ 22; Fig. 1C, 2A). For control conditions in HN (7) neurons, the average spike frequency in the literature is 11.4 6 2.2 Hz .
In all experiments, a sufficiently large increase of g NaP caused the transition from the LA to HA bursting (Fig. 2). We also observed intermediate activities exhibiting mixed bursts with low and high-spike frequencies that resemble a hybrid of both HA and LA bursting (Figs. 1B, 2A, g NaP ¼ 4 nS, black traces). We discarded such activities along with any other activities with high variability of the cycle period (coefficient of variation ! 0.25) from further analysis. The decision to use this cutoff value was informed by the comparison to the whole population, which showed a low degree of variability in the burst cycle period (T). The coefficient of variation (CV) for the cycle period ranged from 0.01 to 0.43, with a median of 0.06 and first and third quantiles of 0.04 and 0.10, respectively. These statistics indicate that the majority of the population had a relatively low degree of variation. The burst duration CV showed a wider range of variation, spanning from 0.02 to 0.95, with a median of 0.11 and first and third quantiles of 0.07 and 0.17, respectively. Similarly, the interburst interval CV ranged from 0.02 to 0.79, with a median of 0.08 and first and third quantiles of 0.06 and 0.17. In HA bursting, I NaP augmented the depolarization phase of the neuron, facilitated [Na 1 ] i influx, and supported increased spike frequency ( Fig. 2A). The strong positive correlation between voltage envelope and [Na 1 ] i oscillation amplitudes demonstrate that membrane potential depolarization and [Na 1 ] i influx grow together (experiment-by-experiment averaged r ¼ 0.98, p , 0.01, pooled data r ¼ 0.766, p , 0.001, n ¼ 22; Fig. 2B). The voltage amplitude and the spike frequency are also strongly positively correlated (experiment-by-experiment averaged r ¼ 0.96, p 0.01, pooled data r ¼ 0.83, p , 0.001, n ¼ 22; Fig. 2C). The scatterplots (Fig. 2B,C) indicate the differences in the voltage amplitude, the intracellular Na 1 concentration amplitude, and the spike frequency between the LA and HA bursting regimes. We used a twothreshold algorithm to classify our data and train a Gaussian naive Bayes supervised machine learning classification algorithm (ML; Materials and Methods). The quality of the classification is evaluated by a report table showing the accuracy of the classification in terms of precision and recall (Table 5). Precision is the ratio of true positives to the sum of true and false positives. The closer the precision score is to 1.0, the better the expected model performance. Recall is the ratio of true positives to the sum of true positives and false negatives. The closer the recall score is to 1.0, the better the expected model performance. The F1 score is the weighted harmonic mean of precision and recall. The closer the value of the F1 score is to 1.0, the better the expected performance of the model. Support is the number of actual occurrences of the class in the dataset, it diagnoses the performance evaluation process.
Looking at the precision and recall scores, one can see that for both classes (LA and HA), the precision score is above 0.9 and the recall score is above 0.9, which indicates that the model is accurately identifying both classes with high precision and recall. The F1-score, which is a weighted average of precision and recall, is also very high for both classes, indicating a high level of overall accuracy. In addition, the support value for each class is relatively high, indicating that the model was tested on a sufficient amount of data to evaluate its performance accurately. Overall, the weighted average F1-score of 0.99 indicates that the classification model is performing at a very high level of accuracy, and the accuracy of 0.99 further supports this conclusion.
To test the accuracy of the ML classification algorithm and evaluate the meaningfulness of the separation between the LA and HA groups, we performed a one-tailed independent samples t test. The results showed a highly significant difference between the LA and HA groups (t ¼ À13.38, p , 0.001, n ¼ 44), indicating that the algorithm was successful in accurately distinguishing between these two bursting patterns.
We mapped the LA and HA bursting regimes distinguished by the classification algorithm on the parameter plane ( g P ; I pump max ; Fig. 2D). While the LA bursting clusters at the bottom left corner of the scatterplots, the HA rhythm spreads out on the upper area of the plots (Fig. 2B,C). A two-parameter map of g NaP and I pump max shows that g P High CV g NaP = 4 nS g NaP = 5 nS g NaP = 5 nS g NaP = 6 nS g NaP = 2 nS g NaP = 3 nS g NaP = 4 nS  Figure 2. Two bursting regimes observed while varying g NaP and I pump max introduced by dynamic clamp. A, This set of voltage traces exemplify how while keeping I pump max at 0.2 nA and setting g NaP above some critical value, causes a transition of the neuron into the HA bursting regime. In this and all subsequent figures I pump and I NaP are dynamic-clamp injected currents and g NaP and I pump max are their parameter values, respectively. The voltage traces are from the same preparation under similar conditions; the only thing that changes between subsequent traces is the g P parameter, the arrow points at the moment when the experimenter changed g NaP from 4 to 5 nS. B, Scatterplot of waveform envelope voltage amplitude and [Na 1 ] i amplitude. C, Scatterplot of voltage amplitude and average intraburst spike frequency. D, Two-parameter map ( g NaP I pump max ) of color-codes parameter sets exhibiting two different bursting regimes: the LA and HA bursting regimes are separated to the left and right areas of the map, blue and red groups, respectively. There is an area where blue and red data points overlap at the border. The size of each circle indicated the number of accepted experiments with the color-coded outcome.
indicates that there were experiments with high variability in cycle period T (high coefficient of variation: K ! 0.27) that were excluded with this parameter set from the analysis. The threshold for the transition from the LA to the HA bursting varies between neurons and may result from animal-to-animal variability (Wenning et al., 2018). B-D, Output from Gaussian naive Bayes machine learning classification algorithm are color-coded by red and blue; red for the HA bursting and blue for the LA rhythm. Gray lines connect data points from the same experiment (voltage traces not shown). The gold line connects the data points from the sample traces (A). Data points are connected in ascending g NaP order, not in the order the parameters were varied. In most experiments, we varied g NaP in a random order, but for illustration purposes, the parameters for the experimental traces in A are shown in ascending order. determines the switch between LA and the HA bursting (Fig. 2D). LA and HA bursting are segregated to different areas where g NaP is lower (LA) or higher (HA). Because of animal-to-animal variability, there is also an area of overlap where both regimes were observed.

In the HA bursting, increase of I pump max decreases HN burst duration and interburst interval
The Na 1 /K 1 pump current is voltage-independent and can be activated during and between bursts. Because the Na 1 /K 1 pump current can thus regulate neuronal activity throughout the burst cycle, it potentially affects both BD and IBI. The Na 1 /K 1 pump current is activated by increased [Na 1 ] i . In both the LA and HA bursting regimes, [Na 1 ] i increases during bursts, but in the HA, the [Na 1 ] i reaches higher values because of the increased joint contribution of injected I NaP and estimated I NaF (Fig. 1). Then, the Na 1 /K 1 pump remains active after the burst termination throughout the IBI. The interaction between I pump and I NaP controls the IBI by providing mutual feedback. Because of its [Na 1 ] i -dependence, the Na 1 /K 1 pump responds to an increase in intracellular sodium concentration ([Na 1 ] i ) and could terminate bursts. Demonstrating that the peaks of the oscillations of I pump coincide with the last spike of the burst, Figure 3 agrees with this mechanism of burst termination. When [Na 1 ] i declines, the I pump decreases, permitting I NaP to depolarize the membrane potential and initiate a new cycle (Fig. 3). The troughs of the oscillations of I pump coincide with the first spike of each burst.
In the HA regime, increasing I pump max speeds up the period of HN bursting (Figs. 3, 4A Figure 3. Increasing the maximal pump activity (I pump max ) speeds up the HA bursting regime. A set of dynamic clamp traces shows that while keeping the HN(7) neuron in the HA bursting regime with dynamic-clamp injected g NaP ¼ 6 nS, incrementing dynamic-clamp injected I pump max decreased the burst duration and interburst interval, and thus the cycle period. Note that the increasing I pump max changes the pump dynamics, by increasing the amplitude of the I pump oscillation but not its median value. Each experimental data sample comprises a set of four signal traces recorded in real time by our dynamic clamp system, from top to bottom, the traces are membrane potential V m (mV), injected I NaP (nA), injected I pump (nA), virtual [Na 1 ] i (mM  (Fig. 4E). The bursting characteristics were not normally distributed and exhibited a positive skew. To assess statistical differences between groups, we performed a nonparametric Mann-Whitney U test. We compared the low and high points used to compute the percent change for various parameters, such as BD, IBI, baseline [Na 1 ] i , [Na 1 ] i amplitude, baseline I pump , and I pump amplitude. Significant differences were found in their test statistic U-values, associated p-values, and sample sizes: BD (U ¼ 156.0, p ¼ 5.5 Â 10 À6 , n ¼ 44), IBI (U ¼ 156.0, p ¼ 5.5 Â 10 À6 , n ¼ 44), Baseline [Na 1 ] i (U ¼ 156.0, p ¼ 5.5 Â 10 À6 , n ¼ 44), [Na 1 ] i amplitude of oscillations (U ¼ 156.0, p ¼ 5.5 Â 10 À6 , n ¼ 44), Baseline I pump (U ¼ 78.0, p ¼ 1.0, n ¼ 44), I pump amplitude (U ¼ 13.0, p ¼ 0.00015, n ¼ 44).

A simple HN neuron model explains how I pump max controls burst duration and interburst interval
We investigated the oscillatory properties of the interaction of I NaP and I pump in a simple model with two dynamic variables, membrane potential, V m , and intracellular Na 1 concentration, [Na 1 ] i . We curve-fitted this 2D model to match experimental data (Materials and Methods). Our simulations emulate experiments following the protocol 2 (Materials and Methods): we systematically varied I pump max from 0.2 to 0.9 nA but kept all other parameters constant, including conductance of injected I NaP , g NaP ¼ 6 nS. The 2D model has a subset of currents (I h , I NaP-native , I K2 , I NaF , I leak-K , and I leak-Na components) incorporated as steadystate currents with their kinetics taken from the canonical model. They are incorporated to represent cell-to-cell variability in experiments. The injected persistent sodium current (I NaP ) and the injected Na 1 /K 1 pump current (I pump ) had the same parameters as in experiments. The optimization algorithm simultaneously matched the waveforms of the experimental voltage envelope and [Na 1 ] i , BD, and IBI with their simulated counterparts. Figure 5 shows an example of a model optimized to one representative experiment, including eight values of I pump max . The simulated BD and IBI follow the same trend of the dependence on I pump max as observed experimentally and are within range of experimentally observed characteristics (relative error: BD ¼ 0.10, IBI ¼ 0.09; Fig. 5A1,A2).
To relate the 2D model dynamics to experimental data, we outlined its phase portrait. We considered a phase plane of the two phase variables ([Na 1 ] i , V m ), where each point represents a certain state of the neuron described by the two variables. We plotted its oscillatory trajectory projected along with the envelopes of experimentally obtained bursting waveforms (Fig. 5B1-B3,C1-C3). The oscillatory motion of the 2D model on the phase plane is governed by Equations 3 and 4 describing the dependence of the speed of change of the phase variables V m and [Na 1 ] I , respectively, on the position of the system on the plane. The nullclines of the system determine borders between domains on the phase plane where the corresponding variables change the sign of their speed, e.g., reverse directions of motion, since by definition a nullcline of a phase variable is a set of points on the phase plane where the derivative of the variable is zero. The intersection(s) between the two nullclines defines the stationary steady state(s) of the 2D system. The V m -nullcline marks the model's states at which the inward currents are in precise balance with outward currents and, thus, the total current is zero (Eq. 3). It separates the domains with negative and positive derivative of V m . The V m -nullcline has a z-shape, the upper and lower branches of z-shaped nullcline attract the trajectories (stable branches) since below them, the derivative of V m is positive and V m increases, while above them, the derivative of V m is negative, and V m decreases (repolarizes; Fig. 5B2). and with the described directions of the speed of V m . The middle (unstable) branch, connecting the two branches at the knee points, repels the trajectories since the directions of the membrane potential derivative are away from it. It establishes a threshold between depolarized and hyperpolarized states of the neuron or phases of oscillatory activity. Similarly, the [Na 1 ] i -nullcline marks the neuron's states at which the inward Na 1 fluxes are in precise balance with outward Na 1 fluxes and the total Na 1 flux is zero (Eq. 4). On the left side of the [Na 1 ] i -nullcline, the speed of [Na 1 ] i is positive and [Na 1 ] i increases, and on the right side of this nullcline, the speed reverses the sign and [Na 1 ] i decreases. The membrane potential V m changes much faster than the intracellular Na 1 concentration, since the factors determining time scales of the equations are more than three orders of magnitude different, C m ( vF. Thus, trajectories almost immediately reach either the upper or lower branch of the z-shaped V m -nullcline, reaching a state close to the balance of the inward and outward currents without notable change in the intracellular Na 1 concentration. The phase point representing the model's state will almost always be located in the vicinity of the stable branches (Fenichel, 1979). If the [Na 1 ] i -nullcline intersects the V m -nullcline at the middle branch, the intersection point is an unstable stationary state. After reaching one of the branches of V m -nullcline, the trajectory will evolve so, that the [Na 1 ] i will slowly change according to the sign of its derivative, while staying in the close vicinity to the V m -nullcline. These dynamics produces a closed periodic orbit (limit cycle) on the phase plane ( Fig. 5B1-B3) and simultaneous oscillations of V m and [Na 1 ] i (Fig. 5C1-C3).
The limit cycle comprises four sections. Two of them are slow and follow the upper and lower branches of the z-shaped V m -nullcline: the upper slow section of the orbit represents the depolarized, spiking phase and the lower slow section represents hyperpolarized, quiescent phase of the bursting cycle. In contrast, the fast sections represent the rapid transitions between spiking and silent phases roughly at the knee points which serve as the thresholds for these transitions (Figs. 5B2, 6C2). Thus, these knee points of the V m -nullcline determine the range and thus amplitude of the oscillations of [Na 1 ] i . The durations of the depolarized (BD) and hyperpolarized (IBI) phases are similarly determined by the [Na 1 ] i amplitude and distinctly determined by the speed of [Na 1 ] i along the corresponding upper and lower sections of the V m -nullcline. Consideration of these two controlling factors explains the dependence of the BD and IBI on I pump max . We investigated how the phase portrait evolved along with the changes of I pump max from 0.2 to 0.9 nA (Fig. 5B1-B3). The range of the oscillations of [Na 1 ] i shortens with the increase of I pump max as the nullclines shift to the left along the [Na 1 ] i axis and the distance between the knee points shrinks (Fig. 5B1-B3; Table 6). The parameter I pump max scales the sigmoidal activation curve of the Na 1 /K 1 pump, so that with larger values of I pump max , the pump current reaches the same range at the smaller corresponding concentration values and more narrow range of [Na 1 ] i (Fig. 6A1,A2). As we increased I pump max from 0.3 to 0.9 nA, the I pump activation curves became steeper while the median [Na 1 ] i and [Na 1 ] i amplitude decreased; Figure 6A1, A2 shows the activation curves for three levels of I pump max . The vertical and horizontal projections delimit the minimum and maximum I pump and [Na 1 ] i , respectively, and facilitate visual comparison between I pump max 0.3, 0.6, and 0.9 nA. This comparison is consistent with the above analysis ( Fig. 4C-F), where we determined that increasing I pump max decreases the [Na 1 ] i median and amplitude, while the amplitude of I pump increases, and the median of I pump remains unchanged. Altogether these analyses demonstrate that increasing I pump max produces a more significant I pump , which compresses and shifts the [Na 1 ] i oscillations. This scaling effect explains the shift of the knee-points of the V m -nullcline toward smaller intracellular Na 1 concentration and the compression of the distance between them (Fig. 6B,C). This factor of the amplitude of [Na 1 ] i predicts simultaneous shortening of the BD and IBI with the same rate relative to the increase of I pump max . The second controlling factor describes how the speed, d Na 1 ½ i dt , along the slow sections of the limit cycle changes with the increase of I pump max . These changes are caused by  Figure 5. A two-dimensional HN model unravels the dynamics of the interaction of the persistent Na 1 and pump currents in a leech heart interneuron with I NaP and I pump augmented by dynamic clamp. Our model explains the speeding up of the HA bursting observed in the dynamic clamp experiments with increasing I pump max . A, Envelope characteristics of the 2D model (cyan) match corresponding bursting characteristics obtained from a single dynamic clamp experiment: increasing dynamic-clamp injected I pump max decreases BD and IBI. These graphs show the effect of increasing I pump max from 0.3 to 0.9 nA, while keeping dynamic clamp injected g NaP at 6 nS on BD (A1) and IBI (A2). The model mimicked the parameter changes of experimental protocol without change any other model parameters. These graphs demonstrate a good fit between experiment and model (BD relative error ¼ 0.1004, IBI relative error ¼ 0.0913). B1-B3, Phase portraits showing model periodic orbits (limit cycles) and nullclines and projected experimental data, accompanied by timeseries of the experimental and model membrane potential, membrane potential envelope, and [Na 1 ] i (C1-C3) for three different values of injected I NaP . C1-C3, Timeseries of experimental (yellow) and simulated (cyan) V m envelope for three levels of I pump max . B2, C2, The V-nullcline is the Z-shaped blue curve and the [Na 1 ] i -nullcline is the magenta curve. The intersection between these nullclines determines location of the unstable steady state within the limit cycle. Their geometry defines the trajectory of the limit cycle of the model (Fig. 6B). the change of the magnitude of the total Na 1 flux (Eq. 4) and can be different for the depolarized and hyperpolarized phases of the cycle. More depolarized membrane potential produces higher Na 1 influx through the fast Na 1 and persistent Na 1 currents and supports the motion along the cycle. With the increased I pump max ; the obtained Na 1 concentration activates stronger pump current creating stronger Na 1 efflux producing slowing down (opposing) effect on the motion during depolarized phase and speeding up (supporting) effect on the motion during hyperpolarized phase of the cycle. Notably, at the hyperpolarized membrane potentials the fast Na 1 current is fully deactivated and the persistent Na 1 current is partially deactivated. Along with the above described compression (Fig. 6B,C) the hyperpolarized section of the cycle expands toward a more hyperpolarized membrane potential where the persistent Na 1 current is further diminished by deactivation, produces smaller Na 1 influx and thus suggests a steeper dependence of the change of speed, d Na 1 ½ i dt , on I pump max over the hyperpolarized section (Fig. 6C). We evaluated relative contributions of this second factor on BDs and IBIs by considering the slow motion of the phase point between the knee points along with the depolarized and hyperpolarized branches of the V m -nullcline, respectively, with the speed d Na 1 ½ i dt . The speed changes along the sections of the branches and, for our approximations, we evaluated it at 200 points (100 points per branch section) distributed with equal intervals in terms of Na 1 ½ i between the knee points (vertical dashed lines), using Equation 4 (Fig. 6B,C). The computed average speeds for these two sections over these points show that with the increase of I pump max , the average speed over the hyperpolarized section increased with a higher rate than over the depolarized one (Table 6). This factor predicts that IBI shortens faster than BD relative to the increase of I pump max .  Approximation of the BD and IBI by the integration of the slow motion along with the slow branches of the V mnullcline roughly matches the values obtained from the simulated cycle with the relative errors of 18% (Table 7; Figs. 7,8). With I pump max growing from 0.2 to 0.9 nA, the amplitude of oscillations of [Na 1 ] i shrinks by 57%, the average speed, d Na 1 ½ i dt , during the burst increases by roughly 11%, whereas during the interburst interval, it increases by 56% (Table 6; Fig. 7). Altogether, the two factors, the dependence of the intracellular Na 1 concentration amplitude and speed on I pump max , explain the mechanism underlying the decrease of BD and IBI with the increase of I pump max in the dynamic-clamp experiments. The amplitude of [Na 1 ] i is the dominating factor determining the burst duration, while both the [Na 1 ] i amplitude and speed determine the interburst interval.

Discussion
All rhythmically active neuronal networks, like the CPG studied here, require robust neuronal mechanisms for burst formation and termination. Associated with bursts are changes in intracellular Na 1 concentration so that Na 1 influxes and effluxes, while oscillatory, are balanced on the cycle-to-cycle basis. Na 1 efflux is mostly taken care of by the Na 1 /K 1 pump, which is a basic cellular engine maintaining the physiological gradients of Na 1 and K 1 ions across the membrane. Because of unbalanced exchange of 3 Na 1 ions moved out for two K 1 ions returned back into the cell, the pump generates an outward current whenever it operates and because it responds by increasing its activity whenever there is Na 1 buildup in the cell it responds dynamically to spiking and bursting activity. Recent studies show that the Na 1 /K 1 pump current contributes to the dynamics of neurons and is a target of neuromodulation Catarsi and Brunelli, 1991;Tobin and Calabrese, 2005;Arganda et al., 2007;Scuri et al., 2007;Hazelwood et al., 2008;Pulver and Griffith, 2010;Forrest et al., 2012;L.N. Zhang et al., 2012L.N. Zhang et al., , 2013H.Y. Zhang et al., 2015;Picton et al., 2017a,b;Sharples et al., 2021;Hachoumi et al., 2022). For example, the pump has been shown to play an important role in sensory information coding (Arganda et al., 2007;Scuri et al., 2007), and the dynamics of short-term memory in a rhythmic motor activity (Pulver and Griffith, 2010 . We collected the [Na 1 ] i amplitude and speed during the burst and interburst interval computed from the 100 points as depicted in Figure 6B and C. To compute the relative change, we divided each [Na 1 ] i amplitude or speed by its correspondent collected at the initial point (I pump max ¼ 0.2 nA).  Figure 7. The estimated BD and IBI preserve the dependence on I pump max observed experimentally and confirmed by the 2D HN model. The estimated BD and IBI from nullcline analysis preserve the trend observed on the BD and IBI from simulations of the 2D model. Importantly, these estimates were computed by interpolating 100 points from the nullclines' geometry and demonstrating that the [Na 1 ] i speed plays an essential role in determining the BD and IBI. To compute these estimates, we used the points illustrated in Figure 6 (V and [Na 1 ] i coordinates) to compute the currents and [Na 1 ] i speed. We divided the change in [Na 1 ] i amplitude by the [Na 1 ] i speed to get the estimated BD and IBI. The estimations were done based on measurements made at 100 points on the depolarized and hyperpolarized branches of V m -nullcline between the knee-points. The average relative error between the simulated and approximate measurements were 18%. al., 2017;Hachoumi et al., 2022). The Na 1 /K 1 pump current plays an atypical role among membrane currents because it is continuously active across the entire voltage range of neuronal activity, while other currents have distinctive voltage ranges of operation and may inactivate. Moreover, because the pump is activity dependent, adjusting to [Na 1 ] i that builds through spiking and bursting activity, it dynamically adjusts to that activity. Thus, in bursting neurons the pump can potentially contribute to cessation of the burst phase and the prolongation of the interburst interval and must be considered as part of the burst mechanism on the level of a single cell dynamics (Li et al., 1996;Barreto and Cressman, 2011;Krishnan and Bazhenov, 2011;Yu et al., 2012;Jasinski et al., 2013;Sharples et al., 2021). Here, using bursting HN interneurons from the leech heartbeat CPG in conjunction with intracellular electrophysiology, dynamic clamp, and modeling, we identify a key interaction with persistent Na 1 current, a current which supports spiking and burst formation and is also active over a very broad voltage range and thus contributes directly to [Na 1 ] i even below spiking threshold. A quantitative analysis of this interaction was possible with our recently developed dynamic clamp implementation estimating Na 1 fluxes and [Na 1 ] i and injecting simulated pump and persistent Na 1 currents into the cell in real time.
Leech HN interneurons are organized into a functional CPG by strong inhibitory synaptic connections, but each HN interneuron is capable of producing rhythmic bursting autonomously when synaptically isolated with bicuculline  or by ganglionic isolation as in the HN(7) interneurons studied here. The dynamics of the HN neurons are well described; the ionic currents and their inhibitory synaptic interactions have been characterized using voltage clamp and these data have been used to build experimentally useful biophysical models of an HN neuron and an HN mutually inhibitory pair or half-center oscillator (HCO; A.A. Hill et al., 2001). This model has demonstrated predictive power for new experiments unveiling roles of interaction of specific ionic currents in control of bursting properties of HN and HN-HCO (A.A. Hill et al., 2001;Cymbalyuk et al., 2002;Simoni et al., 2004;Sorensen et al., 2004;Olypher et al., 2006;Kueh et al., 2016;Ellingson et al., 2021). In HN interneurons of the leech heartbeat CPG, persistent Na 1 and Na 1 /K 1 pump currents (I NaP and I pump ) are native currents that are targets of neuromodulation by endogenous neuropeptides in the leech (Opdyke and Calabrese, 1994;Schmidt et al., 1995;Nadim and Calabrese, 1997;Tobin and Calabrese, 2005;Kueh et al., 2016;Ellingson et al., 2021). By applying the dynamic clamp technique, we showed that their interaction undergirds bursting activity. While allowing manipulation of neuronal activity, invasive recording tools, such as sharp intracellular electrodes, can unintentionally affect measured neuronal activity Pusuluri et al., 2021). Here, we demonstrated that appropriate dynamic clamp implementation of these inward and outward currents (I NaP and I pump ) can alleviate the sharp electrode membrane damage and restore typical bursting of HN neurons (Low-Voltage Amplitude bursting regime) in which neuronal excitability and ability to repolarize are both compromised by the nonspecific leak current with the reversal potential 0 mV introduced by sharp electrode penetration Fig. 1A). I pump responds dynamically to [Na 1 ] i , differentiating it from a steady (or voltage-gated) outward current. Thus, pump activity (as embodied in the parameter I pump max ) modulates excitability in HN neurons; increasing pump activity diminishes excitability while decreasing the pump activity enhances excitability. Our experiments demonstrated that excitability must be preserved by mutually counterbalancing currents if bursting is to be maintained with sharp electrode recordings. This finding suggests that this interaction creates a mechanism which flexibly support bursting activity across a wide range of dynamical parameters. This result also accentuates that the dynamic clamp technique is uniquely valuable for prototyping and testing brain-machine interface techniques (Prinz et al., 2004a;Erazo-Toscano et al., 2021).
Sufficient increase of I NaP and I pump with dynamic clamp causes transition of HN neurons into an intense bursting regime. In every HN neuron that we studied, we found that joint upregulation of I NaP and I pump , specifically by increasing I pump max to at least 0.1 nA and g NaP .5.0 nS, gave rise to this new bursting regime, the high voltage amplitude (HA) bursting (Figs. 1C, 2). In addition, we showed that once in the HA bursting regime, increasing I pump max decreases the burst duration (BD) and the interburst interval (IBI), effectively speeding up the rhythm of the oscillation; this mechanism allows robust and flexible control of the pace of the neuronal bursting (Figs. 4, 5).
Depending on environmental factors or behavioral goals, the leech heartbeat CPG should be able to adjust the pacing of the heartbeat. Our experiments demonstrate that HN neurons are capable of two bursting regimes differentiated by their level of excitability signified by the intraburst spike frequency. In the LA, excitability is low, and bursts are long with low spike frequency. In contrast, in the HA, excitability is high, bursts are short with high spike frequency. Because HN neurons project several segments/ganglia to engage other HN neurons of the CPG and motor neurons through inhibitory synapses and electrical junctions, spikes are the major mode of communication in the network, thus the differences in spike frequency between the LA and HA bursting regimes suggest these bursting regimes could support different functions.
The Na 1 /K 1 pump monitors cellular activity and provides a negative feedback The Na 1 /K 1 pump is generally recognized as playing a critical role in self-regulating neuronal activity (Skou, 1988;Brodsky and Guidotti, 1990;Luger, 1991;Glynn, 1993), but less often is the underlying [Na 1 ] i that it regulates considered as a measure of neuronal activity. By generating outward current in response to raised spiking activity, the Na 1 /K 1 pump produces negative feedback which could stabilize functional regimes. Malfunctions in dynamics of [Na 1 ] i can lead to dynamical diseases (Belair et al., 1995). For example, functional regimes can make transitions into seizure regimes or coexist with them (Hahn and Durand, 2001;Forger and Paydarfar, 2004;Ziburkus et al., 2006;Cressman et al., 2009;Ullah et al., 2009;Fröhlich et al., 2010;Krishnan and Bazhenov, 2011). During seizures the Na 1 /K 1 pump can play a crucial role in termination and in postictal depression Ullah et al., 2009;Fröhlich et al., 2010;Krishnan and Bazhenov, 2011;Krishnan et al., 2015). Based on our results we hypothesize that when neuromodulation targets Na 1 /K 1 pump activity, co-modulation of additional current(s) could mitigate dangerous consequences for the dynamics of the functional regime. For example, in the leech heartbeat CPG, myomodulin, an endogenous leech neuropeptide, reduces the pump's activity and enhances I h ; in synaptically isolated HN neurons, weak bursting rhythms become robust with myomodulin application (Tobin and Calabrese, 2005). The net effect is that myomodulin supports functional bursting, speeds up bursting pattern, and avoids dysfunctional seizure-like and asymmetrical bursting regimes in HN neurons (Tobin and Calabrese, 2005;Kueh et al., 2016;Ellingson et al., 2021). In contrast, enhancing only one current can disrupt functional activity in neurons. The predicted new HA bursting might be invoked in the cases of high neuromodulation where the Na 1 -carrying currents are activated, and pump current is reduced or is challenged with Na 1 influx. This regime has not been demonstrated yet in experiments with the leech heartbeat HN neurons, and further investigation is needed to explore the full range of conditions that can lead to the generation of HA bursting activity in these neurons.
The importance of [Na 1 ] i as the measure of spiking activity in production of the negative feedback is also evident from the studies of the Na 1 -dependent K 1 (I KNa ) currents which read out this proxy and synergize with the pump as suggested by Hage and Salkoff (2012). In the leech heartbeat neuronal model, a similar interaction is reported with the noninactivating K 1 (I K2 ) current providing negative feedback based on integration of membrane potential during spiking. This mechanism also requires coordination of conductance values of the persistent Na 1 and K 1 currents and the leak current (Doloc-Mihu and Calabrese, 2014. Further, [Na 1 ] i activates the Na 1 / K 1 pump, and because the pump activity consumes ATP and oxygen, [Na 1 ] i may act as an yet undescribed homeostatic mechanism. The impact of these [Na 1 ] i fluctuations on cellular activity remains unexplored mainly, but our work suggests that it may be a fruitful avenue of future research. Our work estimated underlying oscillations of [Na 1 ] i that accompany bursting cycles. The amplitude of these rhythmic fluctuations is determined by burst intensity (burst duration, voltage envelope amplitude, and spike frequency; Fig. 4) and thus serve as a proxy for the intensity of the bursting regime. The dynamics of the intracellular Na 1 concentration along with the membrane potential explains the properties of the HA bursting. To describe the core of this mechanisms, we assembled a biophysical model with just two state variables (2D model), i.e., the membrane potential and the intracellular Na 1 concentration ([Na 1 ] i ; Fig. 5). This allowed us to apply phase plane analysis and to investigate self-sustained oscillations. The 2D model captured and explained the changes in each of the two burst characteristics (burst duration, and interburst interval), and two waveform characteristics (voltage envelope and intracellular sodium oscillation amplitudes) within the dynamic-clamp experiments injecting the Na 1 /K 1 pump current scaled by I pump max . The 2D model quantitively reproduced the dependence of the burst characteristics on I pump max measured experimentally. Our model shows that one can estimate the burst duration and interburst interval based on the geometry of the nullclines. These results suggest that during and between the bursts the inward and outward currents are almost in perfect balance and depolarization and repolarization is driven by the relatively slow change in the intracellular Na 1 concentration according to fast dynamics of the membrane potential and the slow dynamics of Na 1 concentration and is well described as dynamics of a relaxation oscillator (Bogoliubov and Mitropolsky, 1961;Rinzel, 1986). This analysis emphasizes that the range of the intracellular Na 1 concentration as the major factor responsible for the simultaneous shortening of the burst duration and interburst interval of HA bursting with the up-scaling of Na 1 / K 1 pump activity, and the disbalance of the Na 1 influx and efflux as the second factor controlling the speed of concentration change during and between the bursts. In HA bursting, I NaP and I pump are near balance during the burst and the interburst interval, driving burst initiation and burst termination. The conventional understanding of control over bursting characteristics by conductance-based currents is fundamentally different; for example, activation of outward I K or inactivation of inward I Ca may terminate a burst, determined by activation/inactivation variables that change over time with a time constant (A.A. Hill et al., 2001;Olypher et al., 2006). In contrast, in the mechanism formed by the interaction of I NaP and I pump , the two currents are balanced, supporting the burst, and the burst terminates through the saddle-node bifurcation in the slow subsystem. Considering Na 1 concentration allowed us to thoroughly characterize the dynamics governing burst characteristics in HN neurons, and with this we predict transition of these neurons into a high excitability bursting regime under the dynamical challenge posed by neuromodulation. We propose that there are at least two different bursting patterns, and that the roles of ionic currents changes with the switch between LA and HA potentially by neuromodulation. We dissected this new HA bursting mechanism down to interaction I NaP and I pump through intracellular Na 1 concentration which allowed us to quantitatively match temporal properties like burst duration and interburst interval along with the change of the maximal pump current. Since this description quantitatively describes changes of the burst duration, we conclude that this interaction controls the termination of the burst and interburst interval. Since, this mechanism produces shorter bursts and interburst intervals, it overrides the mechanisms described for LA bursting (A.A. Hill et al., 2001).