Propensity for Bistability of Bursting and Silence in the Leech Heart Interneuron.

The coexistence of neuronal activity regimes has been reported under normal and pathological conditions. Such multistability could enhance the flexibility of the nervous system and has many implications for motor control, memory, and decision making. Multistability is commonly promoted by neuromodulation targeting specific membrane ionic currents. Here, we investigated how modulation of different ionic currents could affect the neuronal propensity for bistability. We considered a leech heart interneuron model. It exhibits bistability of bursting and silence in a narrow range of the leak current parameters, conductance (gleak ) and reversal potential (Eleak ). We assessed the propensity for bistability of the model by using bifurcation diagrams. On the diagram (gleak , Eleak ), we mapped bursting and silent regimes. For the canonical value of Eleak we determined the range of gleak which supported the bistability. We use this range as an index of propensity for bistability. We investigated how this index was affected by alterations of ionic currents. We systematically changed their conductances, one at a time, and built corresponding bifurcation diagrams in parameter planes of the maximal conductance of a given current and the leak conductance. We found that conductance of only one current substantially affected the index of propensity; the increase of the maximal conductance of the hyperpolarization-activated cationic current increased the propensity index. The second conductance with the strongest effect was the conductance of the low-threshold fast Ca2+ current; its reduction increased the propensity index although the effect was about two times smaller in magnitude. Analyzing the model with both changes applied simultaneously, we found that the diagram (gleak , Eleak ) showed a progressively expanded area of bistability of bursting and silence.

Introduction

Both single neurons and neuronal networks ubiquitously exhibit multistability of activity regimes. As a result of this coexistence of regimes, a neuron or a neuronal network could be switched from one regime to another by noise or by a short transient perturbation such as a pulse of synaptic current (Guttman et al., 1980; Paydarfar et al., 2006). Multistability has been implicated in an assortment of brain functions. In information processing and maintaining short-term memory, multistable neurons can play roles of the memory units and the toggle switches (Marder et al., 1996; Egorov et al., 2002; Williams et al., 2002; Loewenstein et al., 2005; Durstewitz and Seamans, 2006; Tahvildari et al., 2007). In motor control, the multistability of interneurons may pertain to the robust and flexible operation of multifunctional central pattern generators (CPGs), the oscillatory neuronal networks controlling different rhythmic motor patterns, like crawling and swimming in the medicinal leech (Getting, 1989; Marder, 1994; Jing and Weiss, 2001; Sutton et al., 2004; Venugopal et al., 2007; Briggman and Kristan, 2008; Schwabedal et al., 2014; Lyttle et al., 2017). Coexistence of silent and tonic spiking regimes exhibited by spinal motoneurons plays an important role in posture maintenance (Hounsgaard et al., 1984; Conway et al., 1988; Eken and Kiehn, 1989; Hounsgaard and Kiehn, 1989; Kiehn and Eken, 1998; Lee and Heckman, 1998b; Carlin et al., 2000; Perrier and Hounsgaard, 2000). Besides its functional roles, multistability of neuronal networks can accompany pathological conditions such as epilepsy (Foss and Milton, 2000; Hahn and Durand, 2001; Fuentealba et al., 2005; Fröhlich and Bazhenov, 2006; Takeshita et al., 2007; Krishnan et al., 2015). In most cases, the mechanisms underlying multistability and the factors leading to multistability are not well understood.

Invertebrate nervous systems provide the prime advantage of gaining insights into the dynamics of multistable systems on the cellular level since their neurons can be uniquely identified by position, morphology, and electrical activity. The phenomenon of bistability in neuronal dynamics was predicted theoretically and confirmed experimentally in various identified invertebrate neurons. Classic examples are the theoretical studies that revealed the bistability of tonic spiking and silence in the squid giant axon (Rinzel, 1978a; Best, 1979; Guttman et al., 1980; Paydarfar et al., 2006) and multistability of tonic spiking and multiple bursting regimes in the R15 neuron in Aplysia (Canavier et al., 1994; Lechner et al., 1996; Butera, 1998; Newman and Butera, 2010). Here, we focused on bistability of bursting and silence in a well-studied, leech generic neuron, an oscillator leech heart interneuron (HN) (Stent et al., 1979; Calabrese et al., 1995, 2016; Hill et al., 2001; Cymbalyuk et al., 2002; Kueh et al., 2016). In our previous work, we predicted that the leech heart interneuron could exhibit bistability of bursting and silence under conditions with elevated leak conductance, although the range of parameters exhibiting bistability is very narrow (Cymbalyuk et al., 2002; Malashchenko et al., 2011a; Marin et al., 2013). We hypothesize that to reveal this bistability experimentally we need to identify neuromodulatory conditions expanding this range.

Multistability of neurons is usually reported as a result of neuromodulation of membrane currents. For example, in motoneurons bistability results from cellular membrane dynamics, and synaptic input plays the role of a toggle-switch between two regimes (Kiehn and Eken, 1998). Such cellular bistability is an endogenous property of motoneurons that can be recruited when necessary by neuromodulators targeting the dynamics of intrinsic ionic currents (Hounsgaard and Kiehn, 1985; Conway et al., 1988; Lee and Heckman, 1998a, 1999; Shapiro and Lee, 2007). Numerous studies concentrate on assessing the partaking of different ionic currents in the dynamics producing multistability (Hounsgaard and Kiehn, 1989; Hsiao et al., 1997; Lee and Heckman, 1998a; Hughes et al., 1999; Williams et al., 2002). The open question which we address here is: Can different currents be evaluated with respect to the propensity of a single neuron for multistability?

The bifurcation theory provides powerful tools for studies of multistability on the cellular level. It allows one to determine specific biophysical parameters and the range of their variation which would support coexistence of particular regimes of activity. For example, the coexistence of tonic spiking and silence in the dynamics of a single neuron has been extensively studied in models using bifurcation analysis (Rinzel, 1978a,b; Best, 1979; Lee and Heckman, 1998a, 1999; Hahn and Durand, 2001; Williams et al., 2002; Guckenheimer et al., 2005; Tabak et al., 2006; Fernandez et al., 2007; Gutkin et al., 2009; Terman and Ermentrout, 2010; Dovzhenok and Kuznetsov, 2012; Yu et al., 2012; Krishnan et al., 2015). Each observable regime of neuronal activity can be described as an attractor and can be mapped in parameter space by using bifurcation diagrams. These diagrams allow one to determine the borders of parameters where attractors are observed and associate bifurcations of stationary and oscillatory regimes with these borders. The mechanisms underlying multistability specify the unstable regimes that form barriers between the coexisting attractors. Analysis of the bifurcations at which the unstable regimes appear and disappear on the diagrams is a key part of the description of the mechanisms supporting multistability. For example, in cardiac pacemakers and the squid giant axon the threshold separating silent and spiking regimes is the stable manifold of a saddle orbit, which represents unstable subthreshold oscillations (Rinzel, 1978a; Landau et al., 1990). The saddle orbit appears at the subcritical Andronov-Hopf bifurcation and disappears at the saddle-node bifurcation for periodic orbits when injected current is varied. These excitable cells, which were initially exhibiting a single observable regime, could turn bistable after the modulation of membrane ionic currents. For instance, the bifurcation analysis of the Hodgkin-Huxley model predicted that bistability could occur if a steady depolarizing current is applied or if the external K+ concentration is elevated (Rinzel, 1978a; Hahn and Durand, 2001). In a cardiac cell, application of the depolarizing current also leads to bistability (Jalife and Antzelevitch, 1980).

In this article, we investigate the bistability of bursting and silent regimes of a leech heart interneuron. This bistability was reported in our previous studies of a neuronal model (Cymbalyuk et al., 2002; Malashchenko et al., 2011a; Marin et al., 2013). These studies lead us to the problem of identifying biologically relevant conditions under which this bistability can be revealed experimentally. We considered the well-studied model of the leech heart interneuron. It contains eight voltage-gated currents and a leak current (Hill et al., 2001). All currents were measured except for the fast sodium current. Their dynamics were described based on Hodgkin-Huxley formalism and assembled into a canonical model (Hill et al., 2001). To investigate the model behavior, the reversal potential and conductance of the leak current were used as the controlling parameters. Depending on their values, the model exhibited tonic spiking, bursting activity, subthreshold oscillations, or hyperpolarized and depolarized silent regimes. For the narrow range of the leak current parameters, the model exhibited bistability of bursting activity and silence (Cymbalyuk et al., 2002; Malashchenko et al., 2011a). We suggest that the width of this range can be used as an index evaluating a neuronal propensity for bistability. We assume that the larger the area of bistability bounded by the leak current parameters, the more probable it is to find these neurons bistable experimentally.

Since bistability was observed within a small range of the leak conductance, we looked for the conditions that could expand this range. In our previous work we determined the mechanism supporting the coexistence of bursting and silence. Similar to the mechanism of bistability of tonic spiking and silence detected in cardiac cells and the squid giant axon, it was based on the saddle periodic orbit. By varying a controlling parameter, the leak conductance, we found that the saddle orbit originates at a sub-critical Andronov-Hopf bifurcation. The increase of the parameter led to the termination of the saddle orbit through a homoclinic bifurcation. These two bifurcations bounded the range of the leak current parameters where the saddle orbit existed (Malashchenko et al., 2011a) and thus limited the range of coexistence of bursting and silence.

In the present work, we investigated the leech heart interneuron model to ascertain the conditions that increase the propensity of these neurons for bistability. We applied a general method that allows one to identify how up- and down- regulations of currents affect neuronal bistability (Barnett et al., 2013). We used the knowledge about the mechanism supporting the bistability and investigated the effects of modulation of every voltage-gated ionic conductance on the propensity of a neuron for bistability. The results of this article were included into the dissertation of Tatiana Dashevskiy (Malaschenko) (Malaschenko, 2011).

Materials and methods

The biophysical model of the leech heart interneuron is a single compartment model based on Hodgkin-Huxley formalism (Hill et al., 2001). It contains eight voltage-gated conductances describing two types of sodium currents, fast INa and persistent IP, two low-threshold calcium currents, slow inactivating ICaS and fast ICaF, three potassium currents that are delayed rectifier IK1, non-inactivating IK2 and transient IKA, and hyperpolarization-activated cationic current Ih. The dynamics of currents through the membrane is described by a system of 14 differential equations. Parameters of the models had been chosen so that the activity produced by the model is close to the activity recorded experimentally. The bistability of bursting and silence in this model could be induced by up-regulation of the leak conductance. In this paper, it was elevated to 10.7 nS. To obtain trajectories in Figures 1, 2, the integration of models was performed by a variable-step, variable-order methods based on the numerical differentiation formulae by using a Matlab ODES solver ode15s, Mathworks Inc. (Shampine and Reichelt, 1997). Absolute and relative tolerances were 10−8 and 10−9, correspondingly. Initial and maximal step sizes were 0.01 and 1, correspondingly. To demonstrate the bistability of bursting and silence in a single neuron, a negative square pulse of current was applied while neuron was initially in a bursting regime. Pulse was delivered near the first spike of a burst. The duration of the pulse was set to 0.03 s.

Figure 1

Perturbation of the bursting activity of a bistable neuron by a pulse of current can switch the regime into the silent regime. Bursting is shown in dark blue. Switch from bursting into silence can be performed by either depolarizing (A) or hyperpolarizing (B) pulse. Parameters of the model were kept the same as in Cymbalyuk et al. (2002), except for g = 10.7 nS.

Figure 2

Spontaneous transition from bursting to silence at the transition border from bursting to silence. The model exhibited bursting activity for about 500 s and then the activity spontaneously ceased. The leak current parameters were g = 10.8437 nS, Eleak = −0.0635 V. At this value of g, the model was located in parameter space near the homoclinic bifurcation for the saddle periodic orbit. The homoclinic bifurcation appeared at g = 10.87 nS. Initial conditions: V = −0.05485488 V, mCaF = 0.007453999, hCaF = 0.3851188, mCaS = 0.0386471, hCaS = 0.3093507, mK1 = 0.007837126, hK1 = 0.9157689, mK2 = 0.05334662, mKA = 0.1961155, hKA = 0.209315, mh = 0.3366125, mP = 0.1307736, mNa = 0.02026809, hNa = 0.999996, and parameters gleak = 10.8437 nS.

A previous study Malashchenko et al. (2011a) showed that the range of Eleak and gleak where bistability can be observed is limited by the sub-critical Andronov-Hopf and homoclinic bifurcations. These are co-dimension 1 bifurcations. Within this range the unstable periodic orbit exists. Bifurcation diagrams in Figures 3, 4, 7 were computed using CONTENT software which is freely available at http://www.staff.science.uu.nl/~kouzn101/CONTENT/. We integrated the system by using the Runge-Kutta method of the 4-th order. The tolerance of integration was fixed to 10−9 and the minimal step size of integration was set 10−14. The software calculates the Andronov-Hopf bifurcation curves and curves of limit cycles with a given period. In Figures 3, 4 the diagrams were calculated for the pairs of parameters: (g, ḡ), (g, ḡ), (g, ḡ), (g, ḡ), (g, ḡ), (g, ḡ), (g, ḡ), and (g, ḡ). The maximal value of voltage-gated conductance in the diagrams was constrained by the doubled value of the canonical value of the parameter (Hill et al., 2001). If bursting activity disappears for a smaller value than the doubled value of the canonical value, than the smaller value defines the maximal. The minimal value of a conductance in the diagram was 0 nS, unless the bursting activity disappears for a larger value of the conductance. The homoclinic bifurcation curves in Figures 3, 4 were computed using the procedure which includes 1D and 2D bifurcation analysis. First, we calculated a curve of the equilibria. The main bifurcation parameter was g. The maximum step size of the change of g was set to 0.001 nS. We detected the sub-critical Andronov-Hopf bifurcation point where the stable stationary state, representing the silent regime, lost stability. At this point, an unstable sub-threshold periodic orbit collapsed onto the stationary state. Starting from the sub-critical Andronov-Hopf point, we traced the evolution of the unstable periodic orbit. For the orbit continuation, the number of traced points on a periodic orbit was set to 120. The amplitude in the CONTENT window “switch data” was 10−5. As we systematically increased g, the period of the orbit grew logarithmically toward infinity. We continued this procedure until the orbit had period of 100 s. The obtained value of g approximates a homoclinic bifurcation value (Malashchenko et al., 2011a). We used two parameter diagrams to calculate the equi-periodic curve where the model exhibited the unstable periodic orbit with the same period of 100 s. One bifurcation parameter was g and the other parameter was the maximum conductance of a voltage-gated current under investigation. To calculate the curve of this constant period we used the tolerance of 10−9. We denote this curve as an approximation of the homoclinic bifurcation curve.

Figure 3

The two parameter bifurcation diagrams mapping the bursting and silent regimes with the variations of the maximal conductances of IP, INa, ICaF, ICaS vs. gleak. Up- and down-regulation of the conductances of inward currents change the index of propensity which is defined as the width of the range of gleak supporting bistability observed for fixed value of the current conductance ḡ. For each value of ḡ index of propensity was calculated. The blue curve marks the sub-critical Andronov Hopf bifurcation, the red curve represents a homoclinic bifurcation. The saddle periodic orbit exists between these two bifurcations. The yellow dashed curve located near homoclinic curve depicts the value of parameters where transition from bursting to silence occurs. The green curve marks the transition from tonic spiking to bursting. The horizontal dashed black line plotted across the diagram corresponds to the canonical parameter of the voltage-gated current conductance under analysis. (A,B,D) 2D diagrams are plotted in parameters space of (ḡ, gleak), (ḡ, gleak), and (ḡ, gleak). A slim effect on index of propensity was detected. (C) 2D diagram in parameter space of (ḡ, gleak) shows that the decrease of ḡ increases the index of propensity. Eleak is fixed to −0.0635 V.

Figure 4

The two parameter bifurcation diagrams mapping the bursting and silent regimes with the variations of the maximal conductances of IK1, IK2, IKA, Ih vs. gleak. The color and style of the curves are coded as in the Figure 3. (A–C) 2D diagrams are plotted in parameters space of (ḡ, gleak), (ḡ, gleak), and (ḡ, gleak) and no substantial change of the index of propensity is detected. (D) 2D diagram in parameter space of (ḡ, gleak) shows that the increase of ḡ increases the index of propensity. Eleak is fixed to −0.0635 V.

The index of the propensity for bistability was defined as the width of the range of g supporting the bistability of bursting and silence. The index of propensity is bounded by two critical values of g, one determining the subcritical Andronov-Hopf bifurcation and the other corresponding to the transition from bursting to silence. The last is calculated empirically by integrating the system. We varied each conductance one at the time and detected the index of propensity. The slope of the index of propensity vs. g was used to analyze the effect of each current on the propensity for bistability.

Model of leech heart interneuron

A leech heart interneuron is modeled as a single isopotential compartment with the membrane conductances represented in terms of the Hodgkin and Huxley formalism (Hill et al., 2001). The dynamics of membrane potential (V) of a neuron are governed by

where C is the total membrane capacitance (0.5 nF), I is an intrinsic voltage-gated current, I is the injected current and I is the leak current. The voltage-gated currents are given by

where ḡ is the maximal conductance, E is the reversal potential of a voltage-gated current, f(t) is 1 during the pulse, otherwise f(t) is 0, m and h are the activation and inactivation gating variables, respectively. These variables are governed by the following equations:

where the steady-state activation and inactivation functions are given by the following function:

except for the steady-state activation of I which is given by

All time constants are described by the function , except for the inactivation time constants for I and I

The canonical values of the reversal potentials are ENa = −0.045 V, ECaS = 0.135 V, EK = −0.07 V, Eh = −0.021 V, and the maximal conductances are ḡ = 200 nS, ḡ = 7nS, ḡ = 5nS, ḡ = 3.2nS, ḡnS, ḡ = 80nS, ḡ = 80nS, ḡ = 4nS, g = 9.9nS.

Results

Switch of activity from bursting to silence

In our previous work, we investigated the mechanisms supporting bistability of bursting and silence in the canonical model of the leech heart interneuron (Malashchenko et al., 2011a). By performing bifurcation analysis, we showed that for a specific set of parameters of conductance and reversal potential of the leak current (g, E) a silent regime coexisted with bursting activity, and the result was mapped on a two parameter diagram of stationary and oscillatory states. This bistability is supported by a saddle periodic orbit. The stable manifold of this orbit separates basins of attractions of silent and bursting regimes and sets a threshold between them. Figure 1 illustrates the switch from the bursting regime into silence. If the perturbation moved the phase point across the threshold into the basin of attraction of a stable rest state, the bursting activity ceased (Figure 1). If during the perturbation, the phase point did not cross the threshold, the neuron resumed bursting activity. Either a hyperpolarizing or depolarizing pulse of current could produce the switch. The switch could happen if the pulse is appropriately adjusted in terms of amplitude and timing. The saddle periodic orbit appeared at the sub-critical Andronov-Hopf bifurcation when the leak conductance, g, progressively increased, and it terminated at the homoclinic bifurcation (Malashchenko et al., 2011a). For the set of parameters near the Andronov-Hopf bifurcation the amplitude of the saddle orbit was very small and the switch from bursting into silence by perturbation was hard to achieve. The amplitude of the saddle orbit grew as it approached a homoclinic bifurcation, and for the parameter values near the homoclinic bifurcation the switch into the silent regime could be achieved easily. The transition from bursting to silence is characterized by long-lasting transient bursting activity (Figure 2). If the parameters are set close to critical values of the transition, the model of the leech heart interneuron can exhibit bursting activity for several hundreds of seconds and then suddenly switch from bursting into silence. In Figure 2, the model exhibits bursting activity for 500 s and after 500 s the bursting ceases. The activity can be temporally restored by applying a pulse of injected current.

In order to locate the area of bistability in a two parameter bifurcation diagram, we mapped two borders demarcating the transition from silence to bursting and the transition from bursting to silence (Figures 3, 4, 7). The transition from silence to bursting was located by the critical parameters of the Andronov-Hopf bifurcation. The border of the transition from bursting to silence was detected by numerical integration of bursting trajectories for different values of g and E parameters. In order to calculate the parameter values corresponding to the transition from bursting to silence, we had to be cautious about possible co-existence of the bursting and silent regimes; since the neuron's activity depends on the choice of initial conditions. The initial conditions define the state of the neuron as the set of the state variables: the membrane potential, and the activation and inactivation gating variables. In (g, E) bifurcation diagram, each point indicating the transition from bursting to silence was found by numerical integration of the model. First, the initial conditions and g value corresponding to bursting activity were picked. In order to confirm that the bursting was an attractor, we verified that the neuron exhibited bursting for at least 2,000 s. If bursting persisted, we used the set of the state variables at the last phase point of the trajectory of bursting as the initial conditions for the next trial in which we slightly increased g and again performed integration for 2,000 s. The procedure was repeated until the bursting was not observed anymore. The detected critical g was mapped in the 2D diagrams as a transition point from bursting to silence. We assigned E to different values and, by using the described procedure, calculated the corresponding g values of the transition. The collected (g, E) parameters were plotted as the right-side border of the area of bistability. The border is located near the homoclinic bifurcation curve that limits the range of parameter values of the existence of the unstable regime in (g, E) parameter plane. For example, for E = −0.0635 V the bistability was recorded for the g values between 10.67 to 10.84 nS. At the value of g = 10.67 nS the neuron underwent the subcritical Andronov-Hopf bifurcation and the transition from silence to the bursting regime occurs. For values of g larger than 10.84 nS the neuron exhibits only the hyperpolarized silence, the value 10.84 nS indicated the transition point from bursting to the silent regime.

Effects of modulation of ionic conductances on the propensity for bistability

The bifurcation analysis of the bistable system demonstrated the existence of a saddle orbit, and that the range of bistability is limited by values of g corresponding to the Andronov-Hopf and homoclinic bifurcations at which this saddle orbit originates and terminates (Malashchenko et al., 2011a). The transition from bursting to silence was always located between these bifurcations, usually in the vicinity of the homoclinic bifurcation. We traced the Andronov-Hopf and homoclinic bifurcations and mapped them on 2D diagrams. The leak conductance g and a maximal conductance of each ionic current were used as bifurcation parameters. The purpose was to investigate the relative positions of the sub-critical Andronov-Hopf (AH) and homoclinic (Hom) bifurcation curves and the curve (BSi), defining the transition from bursting into silence, while changing g and the maximal conductance of the investigated ionic current, ḡ, where the subscript i identifies the current. This way, we evaluated how up- and down-regulation of a maximal ionic conductance affected the range of the leak conductance where the bursting and silence coexist and where the saddle orbit exists. The difference between the values of g corresponding to the transition from bursting to silence (BSi) and g corresponding to the sub-critical Andronov-Hopf bifurcation (AH) for the same ionic reversal potential value gives the width of the g range where bistability of bursting and silence can be observed. We suggest that this width of the g range can be used as an index evaluating the neuronal propensity for bistability and call it the propensity index H. To find the whole range of parameters where bursting existed, the curve of the transition from bursting to tonic spiking activity (BTs) was also calculated.

We calculated the 2D bifurcation diagrams using the bifurcation software CONTENT for the pairs of parameters, the leak conductance and the maximal conductance of an ionic current. For the inward currents, the 2D diagrams, (g, ḡ), (g, ḡ), (g, ḡ) and (g, ḡ) were calculated, where ḡ, ḡ, ḡ, ḡ are the maximal conductances of the persistent Na+, fast Na+, fast Ca2+ and slow Ca2+ currents (Figures 3A–D). For the bifurcation diagrams we set the upper and lower boundaries of each maximal ionic conductance, ḡ so that the upper boundary of ḡ did not exceed its doubled canonical value and the lower boundary was set to zero nS. If the bursting regime disappeared at the value smaller (larger) than the established upper (lower) boundary, this value was used as the upper (lower) boundary. Figure 3A shows the sub-critical Andronov-Hopf and homoclinic bifurcation curves in (g, ḡ) parameter space; the saddle orbit exists within the range of parameters determined by these curves. The corresponding g width (horizontal slice) for the given ḡ value was very narrow and stayed relatively constant as ḡ was up- or down-regulated from its canonical value of 7 nS. This analysis predicts that the propensity index is not noticeably affected by the change of ḡ. To locate the range of bistability, we calculated the transition from bursting to silence (BSi) that defined the right-side border of the bistable area. The transition BSi curve passes very close to the homoclinic bifurcation curve. The left border of the bistable range is defined by the sub-critical Andronov-Hopf bifurcation curve. The area of bistability in the 2D diagram was approximately located between the Andronov-Hopf and homoclinic curves, and more precisely it was determined by the Andronov-Hopf and BSi curves. The change of ḡ affects the bursting area only slightly that is observed in the width of the range of parameters between BTs and BSi curves. The decrease of g beyond the BTs curve showed that the model can exhibit only tonic spiking activity.

Let us now consider the influence of modulation of the fast Na+ current on the propensity index. The increase of ḡ from the canonical value, 200 nS, up to 400 nS only slightly decreased the width of g range between the sub-critical Adronov-Hopf and homoclinic bifurcation curves (Figure 3B). Similar to the analysis performed for ḡ in Figure 3A we calculated the BTs (transition from bursting to tonic spiking), AH, BSi, and Hom curves. According to the diagram in Figure 3B the bursting area and the area of bistability were little altered when ḡ was changed. The effect of up- or down-regulation of Ca2+ currents on the AH and Hom bifurcation is shown in Figures 3C,D. The down-regulation of ḡ increased the g interval between the two bifurcation curves, whereas the down-regulation of ḡ produced a smaller but opposite effect. Considering the g interval between AH and BSi curve, the propensity index H increased as ḡ decreased, and in contrast H decreased as ḡ decreased.

Further, we investigated the effect of four other currents, one predominantly inward and three outward. The hyperpolarization-activated current, Ih, has the reversal potential of −21 mV and mostly acts as an inward current. Figure 4 shows 2D bifurcation diagrams with different sets of bifurcation parameters: (g, ḡ), (g, g), (g, ḡ) and (g, ḡ), where ḡ, ḡ, ḡ and ḡ are the maximal conductances of the delayed rectifying K+, fast K+, transient K+, and hyperpolarization-activated ionic currents, respectively. The increase of ḡ from 0 to 200 nS had almost no effect on the interval of g values between AH and Hom bifurcation curves (Figure 4A), keeping the index of propensity H almost constant. On the other hand, the range of (g, ḡ) parameters supporting the bursting regime in the model increased significantly as ḡ decreased from the canonical value of 100 to 0 nS. This range is located between BTs and BSi curves. Figures 4B,C show that the changes of ḡ and ḡ conductances had no significant effect on the g interval between AH and Hom bifurcation curves. In the same way, we investigated the effects of modulation of ḡ. The elevation of ḡ strongly affected both the area between AH and BSi and the g interval between the AH and Hom bifurcation curves. For instance, the increase of ḡ from the canonical value of 4 nS to 8 nS expands the g interval between bifurcations by a factor of 3.2 (Figure 4D). Moreover, the range of the parameters where the bursting is observed, i.e., the area between BTs and BSi curves also significantly increased with the elevation of ḡ.

Let us further explore the data presented in Figures 3, 4. The bifurcation diagrams in Figures 3, 4 are difficult to assess visually in terms of the change of the propensity index H. Therefore the data from these diagrams were used to construct graphs showing the dependence of the propensity index on the maximal ionic conductances (Figures 5, 6). The value of the propensity index for the canonical set of parameters H corresponds to 0.17 nS. A maximal ionic conductance ḡ increases the range of bistability of bursting and silence if its up- or down- regulation produced a value of H larger than the canonical value H.

Figure 5

Effect of the variations of the maximal conductances of IP, INa, ICaF, ICaS on the propensity index for bistability of bursting and silence. The propensity index is defined as the range of gleak where bistability of bursting and silence is observed. The dashed yellow vertical line corresponds to the canonical value of a parameter in the model. (A) Manipulation of ḡ almost does not affect the index of propensity. (B) The decrease of ḡ to 100 nS slightly increases the index of propensity from 0.17 to 0.2. (C) The complete reduction of ḡ increased the index from 0.17 to 0.23. (D) The increase of ḡ almost does not affect index of propensity.

Figure 6

Effect of the variations of the maximal conductances of IK1, IK2, IKA, Ih currents on the index of propensity for bistability of bursting and silence. The dashed yellow vertical line corresponds to the canonical value of a parameter in the model. (A) Increase or decrease of ḡ mostly decreases the index of propensity. (B,C) The changes of ḡ and ḡ do not lead to a significant increase of index of propensity. (D) The increase of ḡ significantly increases the index of propensity. The change of ḡ from 4 to 8 nS increases the index of propensity from 0.17 to 0.31.

We investigated dependence of the propensity index H on ḡ (Figure 5A). The value of H stays relatively constant near its canonical value and decreases if ḡ is further up- or down-regulated. The similar analysis established that the increase of ḡ from 100 to 400 nS decreased the propensity index monotonously (Figure 5B). In contrast, when the dependence of the propensity index on ḡ and ḡ was evaluated, we found that the calcium currents affected H differently. H was increased by up-regulation of ḡ and reached a maximum near a ḡ value of 5 nS. In contrast, it was down-regulation of ḡ that elevated H (Figures 5C,D).

A similar procedure was used to summarize the effects of regulation of maximal conductances ḡ, ḡ, ḡ and ḡ on the propensity index. The index varied in a jagged manner with the increase of ḡ from 0 to 200 nS (Figure 6A) and had a maximum near the canonical value of ḡ. The increase of ḡ slightly decreased the propensity index, whereas the increase of ḡ almost did not affect it (Figures 6B,C). Interestingly, the increase of ḡ from 2 to 8 nS significantly elevates the propensity index (Figure 6D). Notice, among all maximal conductances the increase of ḡ affects the most the relative position between AH and Hom bifurcation curves. The data presented in Figures 5, 6 are summarized in Table 1. It contains the maximal and minimal values of the propensity index (HMax and HMin) for the regulation of a given maximal ionic conductance and the difference between these values (ΔH), ΔH = HMax – HMin.

Table 1

Minimal and maximal values of propensity index for the different voltage-gated conductances.

	ḡNa, nS (100–400)	ḡP, nS (4–14)	ḡK1, nS (0–200)	ḡK2, nS (30–140)	ḡKA, nS (0–160)	ḡCaF, nS (0-10)	ḡCaS, nS (2–6.5)	ḡh, nS (2–8)
Hmin	0.1366	0.1366	0.0997	0.1054	0.1571	0.1347	0.1223	0.0751
ḡ(Hmin)	400	14	200	120	160	10	2	2
Hmax	0.1993	0.17	0.1834	0.1832	0.1803	0.2337	0.1929	0.3200
ḡ(Hmax)	100	7	85	40	0	0	5	8
ΔH	0.0627	0.0334	0.0837	0.0778	0.0232	0.099	0.0706	0.2449

Minimal and maximal values are marked in the table by H.

Summing up, we sorted all currents into three groups based on the achievable gain of the propensity index, ΔH: small, intermediate and significant. Comparing the data presented in Table 1, we conclude that the increase of Ip and IKA has the smallest effect on the index of propensity. Manipulation of maximal conductances of currents ICaS, ICaF, IK1, IK2, INa produces an intermediate effect on the propensity index. The up-regulation of Ih leads to the most significant increase of the propensity index for bistability of bursting and silence. This still leaves the question of how the ionic conductance manipulation will be reflected on (g, E) diagram. Here, we assume that the larger the area of bistability in the diagram, the easier it would be to find such neurons bistable.

Changes in the 2d bifurcation diagram of stationary and oscillatory states when Ih and ICaF are modified

The critical values of parameters (g, E) defining the Andronov-Hopf and homoclinic bifurcation curves limit the range of parameters (g, E) where the bistability of bursting and silence can be observed (Malashchenko et al., 2011a). The area of bistability defined in (g, E) bifurcation diagram for the canonical set of ionic current parameters was used as a control. We compare the ranges of bistability under different ionic current modifications. We assume that if up- or down-regulation of ionic currents modifies the area of bistability such that it becomes larger than the control, then the more feasible it will be to record the bistability in an experimental setting.

To test our hypothesis that the up-regulation of ḡ significantly increases the range of the bistability, we computed the bifurcation diagrams using the leak current parameters for the value of ḡ fixed to 8 nS (Figure 7B) and compared it with the control diagram (Figure 7A). The diagrams depict two areas of parameters values, one supporting the bursting activity and the other supporting the bistability of bursting and silence. These diagrams also depict the sub-critical AH and Hom bifurcation curves. The diagram in Figure 7B indicates that the elevation of ḡ significantly increases the area of bistability of bursting and silence in the (g, E) diagram.

Figure 7

Bifurcation diagrams of the oscillatory and rest regimes in parameter space (gleak, Eleak) for different ḡ and ḡ. The sub-critical Andronov-Hopf bifurcation of the hyperpolarized rest state (silent regime) is shown by the blue curve and marks the boundary where silent regime loses stability. The green curve corresponds to the homoclinic bifurcation of the unstable subthreshold oscillations at which the oscillations disappear. The red dots are numerically found points where the transition from bursting activity into silence occurs. The area between the blue curve and the section of the red curve with the dots marks the parameter regime where bursting coexists with rest state. The left section of the red curve locates period-doubling bifurcation of the large amplitude periodic spiking and identifies the transition from tonic spiking to bursting. The diagrams are computed for ḡ = 4 nS (A), ḡ = 8 nS (B), and ḡ = 8 nS and ḡ = 0 nS (C). The increase of ḡ increases index of propensity and the area of bistability. The width of bi-stable area is increased by factor 1.8 for canonical value of Eleak = −0.0635V. In (C) ḡ and ḡ are changed to simultaneously increase the index of propensity. The elimination of fast Ca2+ current and the increase of Ih current remarkably expanded the width of bistable area.

The data presented in Table 1 show that modulations of five currents can be classified as having an intermediately strong effect on the propensity index: ICaS, ICaF, IK1, IK2, INa. Among them, the strongest effect is observed with the down regulation of ḡ. Therefore, we investigated whether the superposition of two modulatory conditions could increase the range of bistability even more. We calculated the (g, E) bifurcation diagram with ḡ up-regulated to 8 nS, ḡ set to 0 nS. Figure 7 shows that the range of bistability is increased in the Figure 7C compared to the Figures 7A,B.

Discussion

Multistability of oscillatory and silent regimes is a common phenomenon exhibited by excitable systems such as neurons and cardiac pacemakers (Jalife and Antzelevitch, 1980; Hounsgaard et al., 1984; Conway et al., 1988; Eken and Kiehn, 1989; Hounsgaard and Kiehn, 1989; Landau et al., 1990; Hsiao et al., 1997; Kiehn and Eken, 1998; Lee and Heckman, 1998a, 1999; Carlin et al., 2000; Perrier and Hounsgaard, 2000). Studying neuronal dynamics on the cellular level allows thorough investigation of the mechanisms supporting bistability, defined as coexistence of two observable regimes, i.e., attractors, in the phase state space. It has been the focus of numerous experimental and theoretical studies (Rinzel, 1978a; Guttman et al., 1980; Canavier et al., 1994; Lechner et al., 1996; Hsiao et al., 1997; Butera, 1998; Lee and Heckman, 1998a, 1999; Crunelli et al., 2005; Cymbalyuk and Shilnikov, 2005; Cymbalyuk et al., 2005; Shilnikov et al., 2005; Fröhlich and Bazhenov, 2006; Newman and Butera, 2010; Malashchenko et al., 2011a,b; Dovzhenok and Kuznetsov, 2012; Barnett et al., 2013; Marin et al., 2013; Krishnan et al., 2015). A classic example of a bistable system is the squid giant axon exhibiting co-existence of tonic spiking and silence under conditions of low Ca2+ bath concentration (Guttman et al., 1980). Spinal motoneurons from the cat, turtle and frog show bistability of tonic spiking and silent regimes (Hounsgaard and Kiehn, 1989; Hsiao et al., 1997). Bistability of tonic spiking and silence was observed in neurons from layer V of the Entorhinal cortex and in the Purkinje neurons in rats (Egorov et al., 2002; Williams et al., 2002; Loewenstein et al., 2005; Tahvildari et al., 2007).

Theoretical analysis of biophysically realistic models can accurately predict whether particular neurons are capable of producing coexisting regimes of activity; and bifurcation analysis allows one to investigate the dynamical mechanisms supporting multistability. Because bistability requires an unstable regime, the mechanisms supporting it can be classified according to the types of separating regimes involved. Typically, a separating regime is either a saddle equilibrium or a saddle orbit. For example, the stable manifold of the saddle equilibrium separates tonic spiking and silence observed in the simplified model of the cerebellar Purkinje neurons (Loewenstein et al., 2005; Fernandez et al., 2007). On the other hand, the stable manifold of the saddle periodic orbit is also a common cause of bistability (Rinzel, 1978a,b; Hahn and Durand, 2001; Shilnikov et al., 2005; Fröhlich and Bazhenov, 2006; Malashchenko et al., 2011a,b; Dovzhenok and Kuznetsov, 2012). Up until recently, most studies have been focused on the bistability of spiking and silence. We are interested in mechanisms of the bistability of bursting and silence that recently has been investigated in a leech heart interneuron model (Malashchenko et al., 2011a,b; Marin et al., 2013).

Here, we show that by using the knowledge of the mechanism, we can elucidate the propensity of the neuron for bistability and the effects of modulations of different ionic currents on it. Since bistability requires presence of an unstable regime, one could asses the range of bistability by determining the bounds of a biophysical controlling parameter where the unstable regime exists. The knowledge of the range of parameters supporting bistability allows us to introduce a measure of the propensity for bistability. We define the propensity index as the width of the range of the controlling parameter where bistability can be observed. As the controlling parameter, one can choose among many biophysical parameters, e.g. the conductance or reversal potential of any ionic current of interest. Because of its special role in the dynamics of excitable cells, following Barnett et al. (2013) in this work, we used the conductance of the leak current as the controlling parameter defining the propensity index. We suggest that the larger the range of values of the leak conductance supporting bistability, the more probable it is to find neurons in the bistable state.

The leak conductance was singled out as such due to the following rationale. The leak current plays an important physiological role in controlling the excitability of neurons. In contrast to other conductances, the leak conductance does not depend on the membrane potential, yet it is a target of modulation. Biophysically, the leak current is carried through multiple channels. Channels with mixed ionic permeability to Na+, Cl−, K+, and Ca2+ ions, channels permeable to K+ such as TASK channels or channels with primarily Na+ permeability such as NALCN are examples of leak channels that control neuronal resting membrane potential (Bayliss et al., 2003; Lu et al., 2007; Koizumi et al., 2010). Depending on the composition of these channels, neuromodulators can reduce or augment the leak conductance. For example, the TASK leak channels are pH-sensitive, and an increase of the extracellular pH up-regulates the leak conductance (Washburn et al., 2002; Larkman and Perkins, 2005). Neurotransmitters including serotonin and noradrenalin close TASK channels, leading to the decrease of the total leak conductance (Sirois et al., 2002; Perrier et al., 2003).

In biophysical models, the activity of leech heart interneurons is highly sensitive to modifications of the leak current (Cymbalyuk et al., 2002; Marin et al., 2013). The notion of such sensitivity is supported by the experiments showing that the activity of leech heart interneurons depends on the method of recording. A leech heart interneuron pharmacologically isolated with bicuculline produces tonic spiking activity if the recording is made intracellularly, whereas it exhibits bursting activity when recording is performed extracellularly (Cymbalyuk et al., 2002). It has been suggested that the intracellular recording with a sharp microelectrode introduces an additional shunting component to the leak current. Similarly, sensitivity to the method of recording has also been reported in Xenopus neurons (Aiken et al., 2003). The intracellular recording with a sharp microelectrode shows that Xenopus neurons fire single short burst or single spikes in response to a depolarizing pulse of current, whereas the whole-cell recording shows tonic spiking activity in response to the same stimulation.

The universality and special role of the leak current parameters has been exploited in analysis of neuronal dynamics in a variety of models (Guckenheimer et al., 2005). Typically low dimensional systems that are based on Hodgkin-Huxley formalism can be dissected in different time scales, slow and fast (Rinzel, 1985; Bertram et al., 1995; Guckenheimer et al., 2005). Such slow-fast decomposition allows consideration of the slow variable as constant within the fast time scale and using it in turn as a controlling parameter. In neuronal models containing slow currents the leak current could be used as a tool to classify the dynamic regime of a model. Slow currents could be approximately considered as constant on a fast time scale, and the slow currents can be lumped together as a component of the leak current (Guckenheimer et al., 2005). Guckenheimer et al. (2005) applied this idea to describe the neuronal dynamics in several models for which the Na+ current is the dominant fast current. The authors considered models of the Aplysia neuron, the thalamic relay neuron, a simplified model of the leech heart interneuron, and a neuron from the respiratory center. They detected bifurcations of a fast subsystem that describe the transitions between spiking and silent phases. The analysis was presented in the form of 2D bifurcation diagrams depicting the regions of oscillatory and stationary states.

In this work, we measured the propensity of a neuron for bistability in terms of the leak conductance. The mechanism has been thoroughly studied in our previous works (Malashchenko et al., 2011a,b). Similarly to the Rinzel mechanism (Rinzel, 1978a), the mechanism supporting bistability of bursting and silence requires the presence of a saddle periodic orbit allowing separation of the stable stationary state (silence) and bursting regime (Malashchenko et al., 2011a,b). The range of controlling parameters supporting the bistability is bounded by the two bifurcations at which the saddle orbit appears and disappears. This orbit exists in a range of leak current parameters that is bounded by the sub-critical Andronov-Hopf and homoclinic bifurcations. In this study, we investigated the effect of up- or down-regulation of each current on the bounds limiting the separating regime and on the propensity index. The impact of every current was classified as small, intermediate or significant (Table 1). Only one current, the hyperpolarization-activated current, was shown to significantly increase the propensity index. The second most influential current, the fast Ca2+ current, ICaF, leads to a moderate increase of the propensity index. Based on our modeling study, we predict that by manipulating these currents experimentally one could reveal the bistability of bursting and silence in a single leech heart interneuron experimentally.

By analyzing how up- or down-regulations of ionic currents affect propensity index, we could assess whether certain neuromodulators promote multistability. In some neurons, it has been shown that neuromodulators could lead to multistability. For example, serotonin induces bistability of spinal motoneurons in different animals (Hounsgaard and Kiehn, 1989; Perrier and Hounsgaard, 2000). Bistability of tonic spiking and silent regimes was found in the presence of serotonin in an in vitro preparation of the spinal cord of the turtle. Serotonin promotes the depolarized silent regime and supports Na+-dependent tonic spiking activity. By blocking the Na+ current with TTX, the depolarized rest state was shown to be maintained by an L-type calcium current and a Ca2+-dependent mechanism. A perturbation of the neuron by a brief pulse of current can switch the regime of activity between the hyperpolarized and depolarized rest states (Perrier and Hounsgaard, 2000; Perrier et al., 2000). In the leech heart interneuron, there are two identified Ca2+ currents, slow and fast. Our analysis of the model indicated that these currents have different influence on the propensity for bistability: the increase of slow Ca2+ conductance weakly increases the propensity for bistability of bursting activity and silence, whereas the increase of fast Ca2+ conductance reduces the propensity.

Through the augmentation of ionic currents, neuromodulators play an important role in controlling motor behavior. In the medicinal leech, the endogenous peptide myomodulin has been shown as a modulator of bursting activity of the leech heart interneurons. The myomodulin decreases the period of bursting activity by increasing the conductance of h-current and decreasing the Na+/K+ pump current (Tobin and Calabrese, 2005). Based on our prediction that h-current leads to a stronger propensity for bistability, we hypothesize that myomodulin will increase the propensity of a heart interneuron for bistability. In other biological systems, the role of the h-current in bistability could be different, e.g., the bistability of tonic spiking and silence of cerebellar Purkinje neurons observed in vitro. Some studies show the blockade of this current induces bistability (Williams et al., 2002) whereas another reported bistability of Purkinje neurons in the presence of h-current (Loewenstein et al., 2005). The ionic mechanism of the Purkinje cell bistability has been suggested to rely on non-inactivating inward currents such as the persistent Na+ current that supports the depolarized tonic spiking regime (Williams et al., 2002). In contrast, although not directly comparable, our results suggest that persistent Na+ current plays an insignificant role in affecting the propensity for bistability of bursting and silence. The apparent discrepancy might be related to the difference in the separating regimes supporting bistability in the Purkinje cell (a saddle rest state) and the leech heart interneuron (a saddle orbit).

Several studies have indicated that multistability may play an important role in organizing multifunctional Central Pattern Generators (CPGs). These CPGs can generate different motor programs such that the switch between different regimes can be achieved either by external stimuli or by neuromodulators (Venugopal et al., 2007; Briggman and Kristan, 2008; Bondy et al., 2015). For example, Briggman and Kristan (2008) suggest that the crawling and swimming CPGs of the leech could be an example of a multifunctional network. The authors show that the vast majority of neurons are involved in both behaviors, swimming or crawling. The modeling of multistable systems can open a new vista onto the understanding of dynamics of currents and control of multistability. The method described in this work can be applied to other neuronal models if the mechanism and the separating regime supporting bistability are known. Concerning our model, we found that only two currents substantially affected the index of the propensity for bistability. These results suggest that under the conditions reported here the coexistence of bursting and silence in an isolated leech heart interneuron could be revealed experimentally.

Author contributions

TD and GC: provided substantial contribution to conception and design, acquisition of data, analysis and interpretation of data, and writing of the paper. Both authors approved the final version for publication.

Conflict of interest statement

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