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).

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 I_Na and persistent I_P, two low-threshold calcium currents, slow inactivating I_CaS and fast I_CaF, three potassium currents that are delayed rectifier I_K1, non-inactivating I_K2 and transient I_KA, and hyperpolarization-activated cationic current I_h. 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⁻⁸ and 10⁻⁹, 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.

A previous study Malashchenko et al. (2011a) showed that the range of E_leak and g_leak 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⁻⁹ and the minimal step size of integration was set 10⁻¹⁴. 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_leak, ḡ_P), (g_leak, ḡ_Na), (g_leak, ḡ_CaF), (g_leak, ḡ_CaS), (g_leak, ḡ_K1), (g_leak, ḡ_K2), (g_leak, ḡ_KA), and (g_leak, ḡ_h). 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_leak. The maximum step size of the change of g_leak 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⁻⁵. As we systematically increased g_leak, 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_leak 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_leak 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⁻⁹. We denote this curve as an approximation of the homoclinic bifurcation curve.

The index of the propensity for bistability was defined as the width of the range of g_leak supporting the bistability of bursting and silence. The index of propensity is bounded by two critical values of g_leak, 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_leak was used to analyze the effect of each current on the propensity for bistability.

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 Ca²⁺ 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 Ca²⁺ 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 Ca²⁺ current, I_CaF, 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 Ca²⁺-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 Ca²⁺ 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 Ca²⁺ conductance weakly increases the propensity for bistability of bursting activity and silence, whereas the increase of fast Ca²⁺ 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.

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.