Six Types of Multistability in a Neuronal Model Based on Slow Calcium Current

Background Multistability of oscillatory and silent regimes is a ubiquitous phenomenon exhibited by excitable systems such as neurons and cardiac cells. Multistability can play functional roles in short-term memory and maintaining posture. It seems to pose an evolutionary advantage for neurons which are part of multifunctional Central Pattern Generators to possess multistability. The mechanisms supporting multistability of bursting regimes are not well understood or classified. Methodology/Principal Findings Our study is focused on determining the bio-physical mechanisms underlying different types of co-existence of the oscillatory and silent regimes observed in a neuronal model. We develop a low-dimensional model typifying the dynamics of a single leech heart interneuron. We carry out a bifurcation analysis of the model and show that it possesses six different types of multistability of dynamical regimes. These types are the co-existence of 1) bursting and silence, 2) tonic spiking and silence, 3) tonic spiking and subthreshold oscillations, 4) bursting and subthreshold oscillations, 5) bursting, subthreshold oscillations and silence, and 6) bursting and tonic spiking. These first five types of multistability occur due to the presence of a separating regime that is either a saddle periodic orbit or a saddle equilibrium. We found that the parameter range wherein multistability is observed is limited by the parameter values at which the separating regimes emerge and terminate. Conclusions We developed a neuronal model which exhibits a rich variety of different types of multistability. We described a novel mechanism supporting the bistability of bursting and silence. This neuronal model provides a unique opportunity to study the dynamics of networks with neurons possessing different types of multistability.


Introduction
Multistability is a fundamental attribute of the dynamics of neurons and neuronal networks [1][2][3][4]. As a feature it appears particularly advantageous for neurons which are part of multifunctional central pattern generators [5][6][7]. Some mechanisms underlying bistability, like the coexistence of tonic spiking and silence and the coexistence of tonic spiking and bursting, have been intensively studied and are well understood [8][9][10][11][12]. Surprisingly, there is a gap in our knowledge of the dynamical mechanisms supporting the bistability of bursting and hyperpolarized silence, bursting and subthreshold oscillations and multistability of bursting, subthreshold oscillations and silence. A classification of mechanisms supporting the multistability of oscillatory and silent regimes is yet incomplete, and remains a fundamental problem for both neuroscience and the theory of dynamical systems.
Bursting is a neuronal oscillatory activity consisting of groups of high-frequency spikes, ''bursts'', separated by intervals of quiescence. It is a basic regime of neuronal activity, which for functionally different neurons can signify either a normal or a pathological state. It is a commonly recorded functional regime of activity of the neurons in the central pattern generators (CPGs): oscillatory neuronal networks executing the motor control of rhythmic movements, like breathing in mammals and the heartbeat in invertebrates such as the medicinal leech [13][14][15][16].
Bursting activity is the result of an interplay of ionic currents which are voltage-gated on various timescales. We envisage a neuron as a slow-fast dynamical system. The coexistence of different attracting regimes of activity, i.e., multistability, is not uncommon for such systems. The qualitative theory of dynamical systems provides a rigorous description of scenarios producing multistability of regimes in the system's dynamics. Exemplary studies by Rinzel (1978) and by Guttman, Lewis and Rinzel (1980) formulated and answered a set of questions which describe a basic scenario of bistability of tonic spiking and silence. It is based on the presence of a repelling periodic orbit separating the basin of attraction of the tonic spiking periodic orbit from the equilibrium representing the silent regime. The scenario also describes the modulation of the neuron's dynamics in response to variations of a bifurcation parameter. According to this scenario, the unstable limit cycle emerges through a subcritical Andronov-Hopf bifurcation and disappears through a saddle-node bifurcation for periodic orbits; both bifurcations define the boundaries for the bistability region. This gives rise to hysteresis and catastrophe-like, fast and non-reversible transitions between silence and tonic spiking as the bifurcation parameter is varied. Guttman, Lewis and Rinzel (1980) showed that a switch between these regimes can be executed by a pulse of current in experiments on the squid giant axon in saline with low calcium concentration.
This study is precipitated by our keen interest in the dynamics of the leech oscillator heart interneurons, which constitute the core of the leech heartbeat timing network. They are found as pairs of mutually inhibitory neurons located in ganglia 3 and 4 [16]. This preparation provides a unique opportunity for studying cellular and network mechanisms of bursting. A leech heart interneuron can be decoupled from its network with bicuculline [17]; its endogenous and network activity is well described by the canonical model [18,19]. This model has been instrumental in predicting that these neurons have endogenous dynamics supporting bursting activity in a single cell, and has showed the high sensitivity of the bursting regime to variations of the leak current. This sensitivity explained why these interneurons show bursting activity while recorded extracellularly, and tonic spiking while recorded intracellularly. A rigorous analysis of the canonical model of a leech heart interneuron is a difficult subject. The model is a system of 14 ordinary differential equations with the variables operating on different time scales. We have numerically obtained a twoparameter bifurcation diagram, mapping oscillatory and stationary regimes on the plane of the leak current's parameters, conductance and the reversal potential of the leak current (g leak ,E leak ). It shows the borders from silence to bursting and from bursting to silence, bounding the zone where bursting and silence coexist [19]. To demonstrate the bistability experimentally in a leech heart interneuron, and to identify the implications of this regime for operation of the leech heartbeat CPG, first we develop its low dimensional model and carry out the analysis of this model.
The reduced model introduced here is based on the fast sodium current, I Na ; the slow, low-threshold calcium current, I CaS ; and the leak current, I leak . We have the following rationale for the choice of the currents. I CaS provides the slowest variable of the 14D canonical model of the leech heart interneuron. It has been suggested for the 14D model that I CaS underlies the bursting activity [19]. Also, we have demonstrated that I CaS introduced through dynamic clamp can reinstate the bursting activity of a tonically spiking leech heart interneuron having all Ca 2z -currents blocked and being isolated from other neurons by the application of saline containing Mn 2z [20]. The significance of this achievement is apparent in light of past experiments showing that these neurons are sensitive to the leak current's parameters, so that we could not intracellularly record bursting activity of a neuron pharmacologically singled out from the network. We showed that an intrinsic mechanism for regulating burst duration might be based on the kinetics of I CaS inactivation and we further corroborate that assertion in this study.
The model presented here makes a complementary study to our previous works on a simplified model since it employs the same fast subsystem supporting spiking as in Cymbalyuk and Calabrese, 2001, but differs in terms of its slow variable [12,[21][22][23][24][25]. That model represented a leech heart interneuron under a blockade of Ca 2z currents along with a partial block of outward currents. It included the fast sodium current, I Na , and non-inactivating slow potassium current, I K2 , so that it was described by a system of three differential equations. The inactivation of I Na and membrane potential constituted the fast subsystem; and the activation of I K2 was the slow variable. The present model allows us to focus the investigation on the potential roles of I CaS in the dynamics of a leech heart interneuron.
We construct a simplified model described by the set of ionic currents and their kinetics I CaS , I Na . We sweep the parameters determining the kinetics of the currents and choose one set which produces activity with temporal characteristics close to those recorded experimentally from leech heart interneurons. Fitting the experimental data is not the primary goal of the article; we bring the model into the ballpark with the experimental observables through simple sweeps over parameters' values. We describe mechanisms supporting multistability in the model under variation of parameters of leak current.

Materials and Methods
We present a simplified leech neuron model containing only the fast sodium (I Na ) and slow calcium (I CaS ) voltage-dependent currents and the leak current (I leak ). We refer to it as model {I CaS I Na }, according to the set of voltage-gated ionic currents which it contains. This model is described by the system of the following four equations: where maximum conductances and reversal potentials of I CaS and I Na are g CaS = 80 nS, g Na = 250 nS, E CaS = 0.135 V, and E Na = 0.045 V, correspondingly; conductance, g leak , and reversal potential, E leak , of the leak current are used as bifurcation parameters; C is the membrane capacitance, C = 0.5 nS. Function f ? (A,B,V ) is a steady-state activation (inactivation) function of a voltage-gated ionic current given by Here B is the half-activation (half-inactivation) membrane potential at which f ?~1 =2. In the model, the activation of I Na is considered to be instantaneous, so that m Na~m One can see that the inactivation of the calcium current, h CaS , is the slowest variable in the model. The bifurcation analysis was performed using the parameter continuation software CONTENT (Kuznetsov YA, Levitin VV, Skovoroda AR (1996) Continuation of stationary solutions to evolution problems in CONTENT. Report AMR9611. Amsterdam: Centrum/Voor Wiskunde en Informatica) freely available at http://www.staff.science.uu.nl/,kouzn101/CONTENT/. The solutions of this model were obtained using the Runge-Kutta method of the 4th order and a variable-order method based on numerical differentiation formulas, implemented as the ode15 s solver in Matlab (MathWorks, Inc.). Absolute and relative tolerances were set to 10 {8 and 10 {9 , respectively.
Because the model is based on a subset of inward currents and hence lacks all outward currents except for the leak current, the balance of inward and outward currents has to be restored. It could be achieved in the sense that we tune up the available parameters so that the model produces the activity with the temporal characteristics of the bursting measured experimentally [19]; the instance of such a model found here we call the canonical model. We took into account basic temporal characteristics of the bursting waveforms such as the spike frequency, burst duration, interburst interval, period, number of spikes and duty cycle. By spike frequency we mean the average frequency of spikes in a burst. The burst duration is measured as the time between the first spike and the last spike in a burst. The interburst interval is the time between the last spike of a burst and the first spike of the following burst. The period of bursting is measured as the time between the first spike of a burst and the first spike of the consecutive burst, which is the sum of the burst duration and interburst interval. The duty cycle is the fraction of the period occupied by the burst duration, i.e. the ration of the burst duration to the period. To explore the waveforms of the bursting we swept values of half-inactivation voltages B h and B hCaS of the two voltage-gated currents. Then, the parameters determining the leak current, g leak and E leak , were used as the bifurcation parameters following Cymbalyuk et al., 2002.

Canonical parameters
A considerable concern about the usability of this model was whether we could adjust it to produce bursting with characteristics close to experimentally recorded ones [19]. The targeted characteristics were: the burst duration was to be between 3.2 and 6.0 sec, the interburst interval was to be within 1.5 and 3.0 sec, the duty cycle was to be between 65.8 and 75.5%, the spike frequency was to be between 6.6 and 14.7 Hz. The model parameters supporting bursting with characteristics in the ballpark with these conditions were easily attained in three steps.
First, we varied the maximum conductances, and the inactivation kinetics for I Na and I CaS . Relative to the 14D model, we set larger values for g g Na~2 50 nS, and g g CaS~8 0 nS to make the corresponding currents more accentuated. Tuning the kinetic parameters B h and B hCaS was motivated by the notion of the window mode of a voltage-gated ionic current [21,26]. In this mode an inactivating current exhibits the properties of a persistent, non-inactivating current in the interval of membrane potentials where the steady-state activation and inactivation curves overlap. In the 4D model, the ''window'' mode of I Na could play a role similar to that of the sodium persistent current, which supports burst duration in the canonical 14D model [18,19]. A similar role would be played by the ''window'' mode of I CaS .
From the studies of the 14D model we can infer that in the 4D model the activation of I CaS would be responsible for the inception of a burst, while the inactivation of I CaS would control the burst termination. This mechanism of burst termination is similar to the one shown for the dynamics of a half-center oscillator assembled from two heart interneurons [18,20].
Second, we examined the activity of the model in response to variations of the half-inactivation voltage of I Na , B h , in the function h ? Na :f ? (500,B h ,V)~1=½1ze 500(V zB h ) . As B h is increased, the curve h ? Na shifts towards the more hyperpolarized values, closing the window voltage interval. This manipulation of B h directly affects the steady-state conductance g ? Na , which produces the window current supporting the burst phase of the bursting activity. The overall decrease of g ? Na decreases the period, the burst duration and the duty cycle of the model (Fig. 1). The bursting activity occurs within the range of B h [0.02888 0.03692] V ( Fig. 1). For some range of values of B h smaller than 0.02888 V, the model shows the tonic spiking activity. In the range of B h between 0.03692 V and 0.03790 V, the model exhibits the stable subthreshold oscillations; and it is silent for B h w0.03790 V. The parameters, g leak = 15.7 nS, E leak = 20.0505 V and B hCaS = 0.06 V fixed for this sweep, were chosen so that the model exhibits coexistence of bursting and silence. This sweep of B h allowed us to adjust the burst duration so that it falls within the experimentally measured range, marked by the grey rectangle in Fig. 1 A, while the interburst interval and spike frequency were not attained. The discontinuity of the graphs of the interburst interval ( Fig. 1 B) and the spike frequency ( Fig. 1 C) near B h = 0.36 V is due to the integer-number difference in the number of spikes.
Third, to adjust the interburst interval, we swept the halfinactivation voltage of I CaS , B hCaS . The increase of B hCaS prolongs the interburst interval of the bursting activity. The interburst interval grows monotonically from 0.53 sec to 3.77 sec as B hCaS is swept from 0.048 V to 0.06 V (Figs. 2 and 3). The graph also shows the effect of variation of B hCaS on the spike frequency within a burst (Fig. 3). Our parameter sweeps of the model were concluded with the following parameters B h = 0.031 V, B hCaS = 0.06 V. This adjusted model exhibits bursting activity with the burst duration of 4.5 sec, the interburst interval of 3.8 sec, the period 8.3 sec, and the duty cycle around 54.6%; the number of spikes per burst is 26, the spike frequency is 5.59 Hz (Fig. 2D). Leak current parameters are g leak = 15.7 nS and E leak = 20.0505 V. We further tuned the model by setting g leak = 15.2 nS, so that the temporal characteristics of bursting activity fit well to the experimental data (Fig. 4).

Results
Equilibria and oscillatory regimes: the (g leak E leak )-

parameter bifurcation diagram
Parameters of the leak current are remarkable targets for modulation of neuronal excitability. For small values of the leak conductance the neuron stays silent at the depolarized equilibrium. This regime of depolarized silence describes a neuron in depolarization block. On the other hand, for large values of the conductance the neuron stays silent at the hyperpolarized equilibrium. The regime of hyperpolarized silence describes a neuron at a hyperpolarized rest state, or in hyperpolarization block. For different intermediate values of the conductance, it exhibits a plethora of oscillatory regimes such as various bursting, tonic spiking, and subthreshold oscillations. These regimes represent different levels of excitability. To inventory this rich variety of the neuronal dynamics, we created a two-parameter (g leak ,E leak ) bifurcation diagram for equilibria and oscillatory regimes, which is shown in Fig. 5. On the diagram we map the areas of the hyperpolarized and depolarized silence (equilibria), the areas of tonic spiking (a stable periodic orbit), bursting activity, and various types of multistability. These areas are determined by these four types of codimension-one bifurcations: the Andronov-Hopf, the saddle-node and homoclinic bifurcations of equilibria, and the saddle-node bifurcations of periodic orbits.
The stable depolarized equilibrium loses stability through the supercritical Andronov-Hopf bifurcation. At the critical value, it gives rise to the stable periodic oscillations, see Figs. 5, 6, 7 and 8. This periodic orbit represents the tonic spiking activity of the neuron. The bifurcation occurs at the curve labeled AH1 in the (g leak ,E leak ) diagram in Fig. 5. At the bifurcation, the periodic orbit is born with zero magnitude and a non-zero frequency v, determined by the imaginary part of the characteristic exponents of the equilibrium state. The sign of first Lyapunov coefficient determines the stability of the new-born periodic orbit [27,28]. It is negative for the supercritical Andronov-Hopf bifurcation. The corresponding bifurcation curve corresponds to the transition from depolarized silence into tonic spiking activity.
A large area in the middle of the bifurcation diagram corresponds to the tonic spiking regime of the model. As the parameter g leak is increased from the bifurcation value, the magnitude of the stable periodic orbit increases. The orbit loses its stability through a period-doubling bifurcation on the curve PD in Fig. 5. With g leak increased further, the unstable orbit disappears through the homoclinic bifurcation. At this value the orbit becomes a homoclinic loop of the saddle equilibrium. The bifurcation occurs with critical values of g H leak and E H leak which are marked by the curve H1 on the diagram (Fig. 5). More precisely, to obtain this curve we exploited the fact that the period of the orbit grows as {lnjg leak {g H leak j near the homoclinic bifurcation, thus it can be arbitrarily large in the vicinity of the bifurcation [27,28]. The curve H1 marks the parameters' values corresponding to tonic spiking with a 25 second period. The area between the curves AH1 and PD supports the periodic tonic spiking in the model. Depending on the level, i.e., the value of the other parameter E leak , a further increase of g leak beyond the border H1 rightward leads to a transition from tonic spiking into either hyperpolarized silence or bursting (Fig. 5).
At a large value of g leak the neuron stays silent at the hyperpolarized equilibrium, which is a stable focus. Lowering g leak makes it unstable through another Andronov-Hopf bifurcation defining the curve AH2 in the parameter plane in Fig. 5. For the most part, the bifurcation is a subcritical one which means that it gives rise to an unstable, subthreshold periodic orbit of a saddle type. However, the segment on the curve AH2 bounded by the points labeled AHG1 and AHG2 corresponds to a supercritical Andronov-Hopf bifurcation, giving rise to the stable subthreshold oscillations. At these two points the first Lyapunov coefficient is zero. Each point locates a Bautin bifurcation, which locates the birth of the saddle-node periodic orbit with zero amplitude and non-zero frequency and has codimension 2. The feature of such a bifurcation is that its unfolding includes a bifurcation curve of saddle-node periodic orbits, that originates from the codimension 2 point. These two points lay on the intersection of the Andronov-Hopf bifurcation AH2 and the saddle-node bifurcation for periodic orbits SN o1 (Fig. 5 and Fig. 6). On the curve SN o1 , two   To the left of AH1, the equilibrium is stable. To the right, it becomes unstable, giving rise to stable tonic spiking. The orbit of spiking loses stability at the period doubling bifurcation, the green curve PD. Followed further, the tonic spiking periodic orbit disappears through a homoclinic bifurcation of an equilibrium marked by the red curve H1. The A-H bifurcation curve (AH2) for hyperpolarized equilibrium is shown in light blue. On this curve, the two points AHG1 and AHG2 mark the Bautin bifurcations ('^'). They bound the section of the curve where the A-H bifurcation is supercritical and gives rise to the stable subthreshold oscillations. The outer sections, above AHG1 and below AHG2, mark the subcritical A-H bifurcation, giving rise to the saddle orbit. The range where the saddle orbit exists is bounded by the homoclinic bifurcation of the saddle equilibrium, the light brown curve H2. Passage through the supercritical section leads to the onset of stable subthreshold oscillations. These oscillations vanish through a saddle-node bifurcation of periodic orbits on the dashed black curve SN o2 (Fig. 6). The area supporting bursting is obtained numerically (mapped in orange, light blue, and partially in pink) and its border is marked by '+'s. This border is bounded by the curves PD and H2. In the pink zone there coexist bursting and stable subthreshold oscillations (Fig. 8). The bright blue patch corresponds to the bistability of bursting and silence; it is bounded by the curves AH2, H2 and H1. The yellow area between the curves AH2 and PD corresponds to the coexistence of tonic spiking and the hyperpolarized silent regime. The dotted lines indicate the four levels of E leak used in the diagrams: 20.048 V (1, Fig. 7-1), 20.04938 V (2, Fig. 7-2), 20.04958 V (3, Fig. 8-3), and 20.0505 V (4, Fig. 8-4). The dark brown curve, SN e , corresponds to the saddle-node bifurcation at which the hyperpolarized equilibria disappear (Figs. subthreshold periodic orbits, stable and saddle, coalesce and vanish. Outside the interval between AHG1 and AHG2, the periodic orbit emerges unstable through the subcritical Andronov-Hopf bifurcation at AH2. With g leak increased, it terminates at the homoclinic orbit of the saddle equilibrium. This event defines the curve H2 (Fig. 5, Fig. 6).
The bursting activity is observed in the model for quite a large range of the leak current parameter values. This range is bounded by '+' in Fig. 5 and the boundaries of the area are special interest here. The onset of bursting appears to be in association with a rapid period doubling cascade leading to chaos in the model [29]. The final event of the scenario involves a homoclinic bifurcation of an unstable periodic orbit becoming the homoclinic loop of the saddle equilibrium. The reduced 4D model also demonstrates the coexistence of long-period, irregular bursting with chaotic tonic spiking oscillations. This bistability is observed in the narrow stripe bounded by the curves PD and H1. This type of transition has been described for the Hindmarsh-Rose model [30].
To expand: (1) The basin of tonic spiking and the hyperpolarized silent state are separated by the stable manifold of the saddle equilibrium.
In the diagram the zone of this bistability is bounded by the subcritical Andronov-Hopf bifurcation curve AH2, and the period-doubling bifurcation curve PD (Fig. 5, Fig. 7.1 for E leak = 20.048 V). (2) The tonic spiking and subthreshold oscillations are separated by the stable manifold of the saddle equilibrium. This bistability area is determined by the range of parameter values supporting the stable subthreshold oscillations, i.e. it is limited by the saddle-node bifurcation curve for the periodic orbits,SN o2 and the supercritical Andronov-Hopf bifurcation curve AH2 (Fig. 6, Fig. 7.2 for E leak = 20.04938 V). (3) The tonic spiking and bursting coexist near the transition border between them (Fig. 5, Fig. 9 C-D). The mechanism underlying this type of multistability has been previously investigated in the Hindmarsh-Rose model [30].
(4) The bursting and subthreshold oscillations are separated by the stable manifold of the unstable subthreshold oscillations. Similarly to (2), the area of this bistability is determined by the range of parameter values supporting the stable subthreshold oscillations, and hence is bounded by the saddle-node bifurcations' curves for the periodic orbits,SN o2 and SN o1 (Fig. 6, Fig. 8.3 for E leak = 20.04958 V, Fig. 9 A-B). (5) The bursting and the hyperpolarized silent regime are separated by the stable manifold of the unstable subthreshold periodic orbit similar to (4). The corresponding zone is bounded by the range of parameters' values supporting the unstable subthreshold oscillations, i.e. the area is bounded by the subcritical Andronov-Hopf bifurcation curve AH2, and the homoclinic bifurcation curve H2 (Fig. 6, Fig. 8.4 for E leak = 20.0505 V, Fig. 11).  (6) The bursting, subthreshold oscillations and silence are separated by the stable manifolds of the two saddle orbits.
One saddle orbit appears through the subcritical Andronov-Hopf bifurcation and disappears through a saddle-node bifurcation for orbits SN o1 . At SN o1 the stable subthreshold oscillations appear (the stable periodic orbit). In turn they disappear at a saddle node bifurcation SN o2 , for a smaller value of the bifurcation parameter, where the second saddle orbit appears, which creates the barrier between the bursting and the stable subthreshold oscillations. The second unstable orbit disappears at the homoclinic bifurcation. Thus the tristability is bounded by either the Andronov-Hopf bifurcation or the second saddle-node bifurcation SN o2 on one side and by the first saddle-node bifurcation on the other side SN o1 . We illustrated the three regimes by switching neuron's activity from the subthreshold oscillations in either bursting or silent regimes (Fig. 10).

Switching between bursting and silent regimes by a pulse of current
In the previous section, we applied the bifurcation analysis to find the range of parameters' values where bursting and silent regimes co-exist. The analysis predicts that a model with parameters taken from this range starting with different initial conditions would demonstrate one regime or the other. These results alone leave a concern as to whether these regimes would be observable. It might be difficult to demonstrate both regimes if the basin of attraction of one of them is too small. For example, for E leak = 20.0505 V, the coexistence is determined for g leak between 15.466 nS and 15.776 nS. If g leak is picked in the vicinity of the Andronov-Hopf bifurcation, the basin of attraction to the equilibrium is small, and it is hard to achieve the switch from bursting into the silent regime. If g leak is picked close to the homoclinic bifurcation, then the situation reverses and it is harder to switch from the silent regime into bursting. With the following Figure 8. The one-parameter (g leak ) bifurcation diagrams corresponding to levels 3 and 4 in Figure 5. The diagrams show the evolution of equilibria and oscillatory regimes for two values of the leak reversal potential: E leak~{ 0:04958V (3) and E leak~{ 0:0505V (4) plotted against the bifurcation parameter g leak . Labeling is the same as in Fig. 7. Numbers 3 and 4 correspond to the dashed lines 3 and 4 in Fig. 5. Here, the blue rectangles determine the range of bistability of the bursting and hyperpolarized equilibrium. In (3) the pink rectangle marks the range of the coexistence of bursting and the stable subthreshold oscillations. The critical values for the Andronov-Hopf bifurcation (AH2) and the saddle-node bifurcation SN o1 are very close to each other and marked by the single vertical line AH2jSN o1 . doi:10.1371/journal.pone.0021782.g008  The question is, how to test for multistability? Although we address this question using a rather abstract model, we have a preference for those perturbations which can be easily implemented in an experimental set up, so that the predictions made could be tested in living neurons. A procedure of this sort has been utilized by Guttman, Lewis & Rinzel (1980). Following their approach, we used perturbations of the model made by ''injection'' of a square pulse of current into the neuron. We tested whether such perturbation could be used to switch from one regime to the other. The duration of the pulses was set to 0.03 sec, which is close to the effective width of a synaptic current pulse in the leech heart interneuron. We tried different amplitudes of the pulse. To switch from the silent regime into bursting, one needs to provide a pulse of current with a sufficiently large amplitude. The pulse can be either depolarizing or hyperpolarizing. For each polarity the critical value of the pulse which is just sufficient for the switch is different. Thus, there are two thresholds for the switch from silence into bursting. The description is more complicated for the switch from bursting into silence; and one more parameter of stimulation has to be taken into account. It is a phase within the period of bursting. Interestingly enough, for the parameters chosen one could switch from bursting into silence again by either a hyperpolarizing or depolarizing pulse of current (Fig. 11).
An important question for experimental testing is whether the unstable oscillations responsible for the bistability of bursting and silent regimes could be recorded. The model study predicts that it might be possible. A projection of a trajectory onto (m CaS , h CaS )plane before and after a pulse was injected into the neuron showed that the phase point oscillated in the vicinity of the unstable subthreshold orbit sufficiently long to trace a few cycles of the unstable oscillations (Fig. 12). Here, the green dot represents the stable hyperpolarized equilibrium, whose basin of attraction is bounded by the stable manifold of the unstable periodic orbit. The orbit is shown in red. Application of a hyperpolarizing I inj = 20.029 nA disturbs the neuron but the state is still in the basin of attraction of the equilibrium. However, application of a pulse of a larger magnitude I inj = 20.03 nA moves the model outside of the basin of the silent attractor away into bursting (Fig. 12). These two values specify a threshold in terms of the critical amplitude of the hyperpolarizing pulse of current.

Discussion
Either a single neuron or neuronal networks can exhibit bistability as demonstrated in a number of theoretical and  experimental studies. The coexistence of neuronal activity regimes -silence, subthreshold oscillations, tonic spiking and burstingwith each other has been observed in experimental studies [1,9,13,31,32]. Multistability has clear implications for dynamical memory and information processing in a neuron [3,4,31,[33][34][35]. In the area of motor control it could be a major mechanism of operation of multifunctional central pattern generators [5][6][7]36]. Multistability can be classified according to regimes which coexist: coexistence of two silent regimes [1], coexistence of tonic spiking and silence [1,2,4,[8][9][10]37], coexistence of bursting and silence [19], coexistence of bursting and subthreshold oscillations [38], coexistence of different tonic spiking regimes [22], coexistence of different bursting regimes [19,39,40], and coexistence of bursting activity and tonic spiking [11,12,19,32,39]. One of the most studied types of bistability is the coexistence of tonic spiking and silence [8,9]. In contrast, the coexistence of bursting and silence has not been the focus of any theoretical or experimental study of dynamics of a single neuron. Here we filled this gap in part by describing multiple scenarios of such coexistence in the dynamics of a low-dimensional model.
One of the first and most extensively studied types of bistability is the coexistence of tonic spiking and silence. The work by Rinzel (1978) set the standard for the investigation of bistability. It found the bistable properties of the Hodgkin-Huxley model of a squid giant axon. In this work, Rinzel formulated a set of questions which have to be answered in order to describe the mechanisms of bistability. What is the unstable, hyperbolic regime which separates the two observable, attracting regimes? How does this unstable regime appear and disappear as controlling bifurcation parameters are changed? Rinzel demonstrated that the unstable subthreshold oscillations are the regime which separates the tonic spiking periodic orbit from the equilibrium representing silence. The unstable oscillations are born through the Andronov-Hopf bifurcation and terminate at the saddle-node bifurcation for periodic orbits, where the stable orbit is the one corresponding to the tonic spiking. These two bifurcations define the range of the bistability in terms of the controlling parameter. Guttman et al. (1980) confirmed experimentally that the squid giant axon exhibits this type of bistability under low Ca 2z bath concentration. Repetitive firing and silence regimes can co-exist; and perturbation by a pulse of current can be used to switch between two regimes (Guttman et al., 1980). Bistability is of particular interest for motor control, since the coexistence of tonic spiking and silence has been shown in cerebellar Purkinje cells [4,37] and has been induced by application of serotonin in spinal motoneurons from different species including the cat, rat, and mouse [1,2].
In our previous studies, a 3D reduced model of the leech heart interneuron exhibited four types of bistability under different parameter regimes: (1) tonic spiking and bursting; (2) two different tonic spiking regimes; (3) tonic spiking and hyperpolarized silence; (4) bursting and depolarized silence [12,22,24]. In all four types, the regimes of activity are separated in the 3D phase space by the 2D stable manifold of the saddle type regime, either periodic orbit or equilibrium state. The topology of the manifold determines the threshold between the co-existing regimes.
Here, we developed a low dimensional, 4D model, which exhibits temporal characteristics close to those recorded from leech heart interneurons. Although the model is simple, it maintains a number of biophysical correspondences to the living counterpart. With this model we explored and accentuated the role of the low threshold slowly inactivating calcium current, as one sufficient to support a plethora of different mechanisms supporting the coexistence of different regimes of activity. It shows six types of multistability: (1) tonic spiking and the hyperpolarized silent regime; (2) tonic spiking and subthreshold oscillations; (3) tonic spiking and bursting; (4) bursting and subthreshold oscillations; (5) bursting and the hyperpolarized silent regime; and (6) bursting, subthreshold oscillations, and silence. We illustrated some of these mechanisms by a series of numerical experiments to show that they are sufficiently robust to be observed. We showed that switching between bursting activity and silence can be controlled by a pulse of current.
If we compare the lists of the mechanisms described for the two models, 3D and 4D, we observe one mechanism in common, the one underlying the co-existence of tonic spiking and hyperpolarized silent regimes which is based on a saddle equilibrium. We hypothesize that under different parameter regimes the two lists would record more mechanisms in common. The mechanisms of multistability described for these two models are generic and could be found in a variety of models under different parameter regimes.
Rhythmic motions of animals like swimming, crawling, walking, scratching, heart beating, and ingestion or rejection of food are controlled by specialized oscillatory neural networks, CPGs. It is a fundamental question whether each of these types of motion is controlled by a separate pattern generator. For a number of circuits it has been shown that some interneurons participate in multiple tasks, being part of a multifunctional central pattern generator [5][6][7]41]. This solution appeals as much more parsimonious, and thus evolutionarily advantageous, compared to the one recruiting specialized circuits for each type of behavior. On the other hand, it is also appealing for mathematical modeling from the perspective of the theory of dynamical systems since the behavior of oscillatory neural networks, especially with heterogenous membrane and synaptic properties, is not well understood [24,36]. Bistability as a membrane property of single neurons appears to be a natural feature of building blocks for such CPGs. In the future we plan to use this new model to explore the advantages and disadvantages of different types of multistability in small cental pattern generator circuits.