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

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 (Cymbalyuk et al., 2002). 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 (Imaxpump = 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 Imaxpump 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 Imaxpump 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⁺]_i became visibly higher (Fig. 1C,D2). These observations suggested the emergence of a new bursting regime as a result of the variation of Imaxpump and g¯NaP . Below, we call this regime the high-amplitude (HA) bursting.

We systematically varied dynamic-clamp parameters, Imaxpump 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⁺]_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⁺]_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 ± 2.2 Hz (Cymbalyuk et al., 2002).

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⁺]_i influx, and supported increased spike frequency (Fig. 2A). The strong positive correlation between voltage envelope and [Na⁺]_i oscillation amplitudes demonstrate that membrane potential depolarization and [Na⁺]_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⁺ concentration amplitude, and the spike frequency between the LA and HA bursting regimes. We used a two-threshold 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, Imaxpump ; 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 Imaxpump shows that g¯P 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 Imaxpump decreases HN burst duration and interburst interval

The Na⁺/K⁺ pump current is voltage-independent and can be activated during and between bursts. Because the Na⁺/K⁺ pump current can thus regulate neuronal activity throughout the burst cycle, it potentially affects both BD and IBI. The Na⁺/K⁺ pump current is activated by increased [Na⁺]_i. In both the LA and HA bursting regimes, [Na⁺]_i increases during bursts, but in the HA, the [Na⁺]_i reaches higher values because of the increased joint contribution of injected I_NaP and estimated I_NaF (Fig. 1). Then, the Na⁺/K⁺ 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⁺]_i-dependence, the Na⁺/K⁺ pump responds to an increase in intracellular sodium concentration ([Na⁺]_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⁺]_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 Imaxpump speeds up the period of HN bursting (Figs. 3, 4A,B). Since the parameter Imaxpump describes the maximal output of the I_pump but does not affect the pump’s sensitivity or the concentration of half-activation of I_pump, Imaxpump determines [Na⁺]_i and I_pump. Namely, an increase of Imaxpump from 0.3 to 0.7 nA decreased the median [Na⁺]_i by 18% [for Imaxpump = 0.3 nA, 90% confidence interval (CI): 12.59 13.35, range: 11.69–14.37, n = 13; for Imaxpump = 0.7 nA, 90% CI: 10.32–10.86, range: 9.86–11.38, n = 13] and the [Na⁺]_I amplitude by 27% (for Imaxpump = 0.3 nA, 90% CI: 1.75–2.44, range: 1.34–3.85, n = 13; for Imaxpump = 0.7 nA, 90% CI: 1.35–1.94, range: 1.02–2.82, n = 13; Figs. 3, 4C,D), and increased the amplitude of oscillations of I_pump by 27% (for Imaxpump = 0.3 nA, 90% CI: 0.07–0.10, range: 0.05–0.16, n = 13; for Imaxpump = 0.7 nA, 90% CI: 0.12–0.15, range: 0.08–0.19, n = 13; Fig. 4F) while not changing the median I_pump (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⁺]_i, [Na⁺]_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⁻⁶, n = 44), IBI (U = 156.0, p = 5.5 × 10⁻⁶, n = 44), Baseline [Na⁺]_i (U = 156.0, p = 5.5 × 10⁻⁶, n = 44), [Na⁺]_i amplitude of oscillations (U = 156.0, p = 5.5 × 10⁻⁶, n = 44), Baseline I_pump (U = 78.0, p = 1.0, n = 44), I_pump amplitude (U = 13.0, p = 0.00015, n = 44).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Increasing dynamic-clamp injected Imaxpump decreased burst duration, interburst interval, and [Na⁺]_i in the HA bursting regime. The average trend across experiments in the HA bursting regime shows that increasing dynamic-clamp injected Imaxpump decreases the following characteristics, burst duration (A) and interburst interval (B), the median [Na⁺]_i (C), and [Na⁺]_i amplitude (D). However, increasing Imaxpump does not affect the median I_pump (E) but increases the amplitude of the I_pump oscillations (F). Changes of different characteristics are shown as normalized pooled data from multiple experiments, represented by different colors consistent between subplots. The black line consists of the mean value across experiments for a given Imaxpump .

A simple HN neuron model explains how Imaxpump 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⁺ concentration, [Na⁺]_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 Imaxpump 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 steady-state 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⁺/K⁺ 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⁺]_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 Imaxpump . The simulated BD and IBI follow the same trend of the dependence on Imaxpump as observed experimentally and are within range of experimentally observed characteristics (relative error: BD = 0.10, IBI = 0.09; Fig. 5A1,A2).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

A two-dimensional HN model unravels the dynamics of the interaction of the persistent Na⁺ 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 Imaxpump . A, Envelope characteristics of the 2D model (cyan) match corresponding bursting characteristics obtained from a single dynamic clamp experiment: increasing dynamic-clamp injected Imaxpump decreases BD and IBI. These graphs show the effect of increasing Imaxpump 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⁺]_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 Imaxpump . B2, C2, The V-nullcline is the Z-shaped blue curve and the [Na⁺]_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).

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⁺]_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⁺]_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⁺]_i-nullcline marks the neuron’s states at which the inward Na⁺ fluxes are in precise balance with outward Na⁺ fluxes and the total Na⁺ flux is zero (Eq. 4). On the left side of the [Na⁺]_i-nullcline, the speed of [Na⁺]_i is positive and [Na⁺]_i increases, and on the right side of this nullcline, the speed reverses the sign and [Na⁺]_i decreases. The membrane potential V_m changes much faster than the intracellular Na⁺ concentration, since the factors determining time scales of the equations are more than three orders of magnitude different, Cm ≪ 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⁺ 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⁺]_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⁺]_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⁺]_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⁺]_i. The durations of the depolarized (BD) and hyperpolarized (IBI) phases are similarly determined by the [Na⁺]_i amplitude and distinctly determined by the speed of [Na⁺]_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 Imaxpump .

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Amplitude of [Na+]i in experimental data, and model nullclines in the phase plane. A1, Black curves are activation functions for three levels of Imaxpump : 0.3, 0.6, and 0.9 nA, showing I_pump as a function of [Na⁺]_i, color coded for each Imaxpump . Vertical projections to the left demark maximum and minimum I_pump values. Horizontal projections to the bottom demark maximum and minimum [Na⁺]_i values. A2, Zoomed inset of A1. These curves and ranges show the dependence of Na⁺/K⁺ pump activity as a function of [Na⁺]_i. Consistent with Figure 5, the range of [Na⁺]_i at Imaxpump 0.6 and 0.9 nA decreased relative to 0.3 nA, while the range of I_pump increased. B, C, Blue curve is the V-nullcline and magenta curve is the [Na⁺]_i nullcline. Green filled circles are 100 interpolated data points used to estimate the speed of [Na⁺]_i, d[Na+]idt , during the burst duration; orange circles are the 100 data points used to estimate the rate of change of [Na⁺]_i during the interburst interval. Vertical dashed lines project from the nullcline knee points and delimit the maxima and minima of [Na⁺]_i and thus its amplitude. B, Imaxpump : 0.2 nA. C, Imaxpump : 0.9 nA.

We investigated how the phase portrait evolved along with the changes of Imaxpump from 0.2 to 0.9 nA (Fig. 5B1–B3). The range of the oscillations of [Na⁺]_i shortens with the increase of Imaxpump as the nullclines shift to the left along the [Na⁺]_i axis and the distance between the knee points shrinks (Fig. 5B1–B3; Table 6). The parameter Imaxpump scales the sigmoidal activation curve of the Na⁺/K⁺ pump, so that with larger values of Imaxpump , the pump current reaches the same range at the smaller corresponding concentration values and more narrow range of [Na⁺]_i (Fig. 6A1,A2). As we increased Imaxpump from 0.3 to 0.9 nA, the I_pump activation curves became steeper while the median [Na⁺]_i and [Na⁺]_i amplitude decreased; Figure 6A1,A2 shows the activation curves for three levels of Imaxpump . The vertical and horizontal projections delimit the minimum and maximum I_pump and [Na⁺]_i, respectively, and facilitate visual comparison between Imaxpump 0.3, 0.6, and 0.9 nA. This comparison is consistent with the above analysis (Fig. 4C–F), where we determined that increasing Imaxpump decreases the [Na⁺]_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 Imaxpump produces a more significant I_pump, which compresses and shifts the [Na⁺]_i oscillations. This scaling effect explains the shift of the knee-points of the V_m-nullcline toward smaller intracellular Na⁺ concentration and the compression of the distance between them (Fig. 6B,C). This factor of the amplitude of [Na⁺]_i predicts simultaneous shortening of the BD and IBI with the same rate relative to the increase of Imaxpump .

The second controlling factor describes how the speed, d[Na+]idt , along the slow sections of the limit cycle changes with the increase of Imaxpump . These changes are caused by the change of the magnitude of the total Na⁺ flux (Eq. 4) and can be different for the depolarized and hyperpolarized phases of the cycle. More depolarized membrane potential produces higher Na⁺ influx through the fast Na⁺ and persistent Na⁺ currents and supports the motion along the cycle. With the increased Imaxpump,  the obtained Na⁺ concentration activates stronger pump current creating stronger Na⁺ 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⁺ current is fully deactivated and the persistent Na⁺ 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⁺ current is further diminished by deactivation, produces smaller Na⁺ influx and thus suggests a steeper dependence of the change of speed, d[Na+]idt , on Imaxpump 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+]idt . 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+]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 Imaxpump , 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 Imaxpump .

Approximation of the BD and IBI by the integration of the slow motion along with the slow branches of the V_m-nullcline roughly matches the values obtained from the simulated cycle with the relative errors of 18% (Table 7; Figs. 7, 8). With Imaxpump growing from 0.2 to 0.9 nA, the amplitude of oscillations of [Na⁺]_i shrinks by 57%, the average speed, d[Na+]idt , 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⁺ concentration amplitude and speed on Imaxpump , explain the mechanism underlying the decrease of BD and IBI with the increase of Imaxpump in the dynamic-clamp experiments. The amplitude of [Na⁺]_i is the dominating factor determining the burst duration, while both the [Na⁺]_i amplitude and speed determine the interburst interval.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

[Na⁺]_i amplitude decreases 57%, the [Na⁺]_i speed (dNadt) during the burst increases 11%, and the [Na⁺]_i speed during the interburst interval increases significantly 44% as a function of Imaxpump . These curves illustrate the relative change of the [Na⁺]_i amplitude and [Na⁺]_i speed as a function of Imaxpump . We collected the [Na⁺]_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⁺]_i amplitude or speed by its correspondent collected at the initial point (Imaxpump= 0.2 nA).

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

The estimated BD and IBI preserve the dependence on Imaxpump 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⁺]_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⁺]_i coordinates) to compute the currents and [Na⁺]_i speed. We divided the change in [Na⁺]_i amplitude by the [Na⁺]_i speed to get the estimated BD and IBI.

The Na⁺/K⁺ pump monitors cellular activity and provides a negative feedback

The Na⁺/K⁺ 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⁺]_i that it regulates considered as a measure of neuronal activity. By generating outward current in response to raised spiking activity, the Na⁺/K⁺ pump produces negative feedback which could stabilize functional regimes. Malfunctions in dynamics of [Na⁺]_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; Fröhlich and Bazhenov, 2006; 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⁺/K⁺ pump can play a crucial role in termination and in postictal depression (Cressman et al., 2009; 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⁺/K⁺ 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⁺-carrying currents are activated, and pump current is reduced or is challenged with Na⁺ 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⁺]_i as the measure of spiking activity in production of the negative feedback is also evident from the studies of the Na⁺-dependent K⁺ (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⁺ (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⁺ and K⁺ currents and the leak current (Doloc-Mihu and Calabrese, 2014, 2016). Further, [Na⁺]_i activates the Na⁺/K⁺ pump, and because the pump activity consumes ATP and oxygen, [Na⁺]_i may act as an yet undescribed homeostatic mechanism. The impact of these [Na⁺]_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⁺]_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⁺ 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⁺ concentration ([Na⁺]_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⁺/K⁺ pump current scaled by Imaxpump . The 2D model quantitively reproduced the dependence of the burst characteristics on Imaxpump 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⁺ concentration according to fast dynamics of the membrane potential and the slow dynamics of Na⁺ 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⁺ 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⁺/K⁺ pump activity, and the disbalance of the Na⁺ 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⁺ 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⁺ 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).