Abstract
Background
Development of effective and plausible numerical tools is an imperative task for thorough studies of nonlinear dynamics in life science applications.
Results
We have developed a complementary suite of computational tools for twoparameter screening of dynamics in neuronal models. We test a ‘bruteforce’ effectiveness of neuroscience plausible techniques specifically tailored for the examination of temporal characteristics, such duty cycle of bursting, interspike interval, spike number deviation in the phenomenological HindmarshRose model of a bursting neuron and compare the results obtained by calculusbased tools for evaluations of an entire spectrum of Lyapunov exponents broadly employed in studies of nonlinear systems.
Conclusions
We have found that the results obtained either way agree exceptionally well, and can identify and differentiate between various fine structures of complex dynamics and underlying global bifurcations in this exemplary model. Our future planes are to enhance the applicability of this computational suite for understanding of polyrhythmic bursting patterns and their functional transformations in small networks.
1 Introduction
Individual and networked neurons can generate various complex oscillations known as bursting, formed by alternating fast repetitive spiking and quiescent or subthreshold oscillatory phases. Bursting is a manifestation of composite, multiple time scale dynamics observed in various fields of science as diverse as food chain ecosystems, nonlinear optics, medical studies of the human immune system, and neuroscience. The role of bursting is especially important for rhythmic movements determined by Central Pattern Generators (CPG). Many CPGs can be multifunctional and produce polyrhythmic bursting patterns on distinct time scales, like fast swimming and slow crawling in leeches [1]. Such CPGs are able to switch between different rhythms when perturbed [2,3].
In mathematical neuroscience a deterministic description of endogenously oscillatory activities, like twotime scale bursting, is done by revealing generic properties of mathematical and realistic models of neurons; the latter are derived through the HodgkinHuxley formalism for gating variables. Either bursting model falls into a class of dynamical systems with at least two time scales, known as slowfast systems.
Configurations and classification schemes for bursting activities in neuronal models first proposed in [4] and extended in [5,6] are based on geometrically transparent mechanisms that initiate and terminate the socalled slow motion manifolds composed of the limiting solutions, such as equilibria and limit cycles, of the fast subsystem of a model [711]. These manifolds constitute the backbones of bursting patterns in a neuronal model. A typical HodgkinHuxley model possesses a pair of such manifolds [4]: quiescent and tonic spiking. The existing classifications of bursting are based on codimensionone bifurcations that initiate or terminate the fast trajectory transitions between such onedimensional [1D] and twodimensional [2D] slow motion manifolds in the phase space of a model. These classifications single out the classes of bursting by subdividing neuronal models into the following types: elliptic or Hopffold; squarewave burster, or foldhomoclinic; parabolic, or circlecircle class describing tophat models. These terms are either due to specific shapes of voltages traces in time, or after the static underlying bifurcations that occur in the fast subsystem of the given neuron model.
The types of the static bursting configurations in the HindmarshRose model shown in Figures 1 and 2 are also named fold/homoclinic and fold/Hopf, as this would indicate that the terminal phases of the fast spiking and slow quiescent periods are determined, respectively, by a homoclinic bifurcation of a saddle equilibrium state, or a supercritical AndronovHopf bifurcation, along with a saddlenode bifurcation of equilibria, respectively, which all occur in the fast subsystem of the model. In the next section we will examine the transition bifurcation patterns between these types of bursting.
Fig. 1. (A) Squarewave and (B) plateaulike bursting traces in the HindmarshRose model at , respectively. Transformations of bursting can be detected quantitatively by a sudden change of the number of spikes per burst.
Fig. 2. Top: a plateaulike or fold/Hopf bursting starts after the spiking manifold becomes tangent to the middle, saddle branch of and terminates further through the reverse supercritical AndronovHopf bifurcation on the upper depolarized branch of . Bottom: the primary feature of the squarewave bursting activity in the HR model also referred to as of fold/homoclinic type is the termination of the spiking manifold by the homoclinic bifurcation in the phase space of the fast subsystem. In both cases: fold stands for a saddlenode bifurcation at the turning point (SN) on the lower, hyperpolarized branch of .
These manifolds, especially their stable branches, can be easily traced and visualized in the phase space by utilizing the slow variable as a sweeping parameter in the decoupled fast subsystem. Far from bifurcations, this slowfast dissection approach allows for exhaustive simplifications that let one treat the dynamics of the full model as on overlay the uncorrelated dynamics of its fast subsystem mediated by repetitive passages of the slow variable.
The slowfast dissection has been proven to work very well for a loworder model of a bursting neuron as long as it stays away from a bifurcation that is due to reciprocal interactions of the dynamics of both subsystems. Such a bifurcation, underlying a bursting transition(s), gives rise to the emergence of dynamical phenomena that can only occur in the full system. For example, this occurs when the dynamics of the fast subsystem falls to the time scale of the slow subsystem, particularly near saddlenode and homoclinic bifurcations. A classic example is the onset of chaotic dynamics of finite subshift shown by D. Terman [12] at the transition between tonic spiking and bursting that happens to be generic for a squarewave burster like the HindmarshRose model [1315]. In addition [12] gives an explanation of common spike adding cascade in classical squarewave bursters which is due to a slow passage of the phase point near the saddle in the fast subsystem. Note that the nature of spike adding cascade could be bifurcationally different, as for example in the leech interneuron model due to the blue sky catastrophe [16] prolonging the burst duration phase, or because of homoclinics of a saddle periodic orbit [17,18] playing the role of the chaotic potential barrier, loosely speaking, that bursting should overcome to get an extra spike.
Complex dynamics, including quick period doubling cascades in square wave bursters [19,20] can also be explained in terms of codimensiontwo homoclinic bifurcations, including inclinationswitch and orbitflips, that occur at the transition [21,22] from tonic spiking to bursting. Recent breakthrough examples of novel transitions to and from bursting due to reciprocal interactions of slowfast dynamics include various homoclinic of saddles and saddlefoci, the blue sky catastrophe, bistability due nontransverse homoclinics to saddlenode periodic orbits, canardtori [13,1618,2326]. The range of bifurcations and dynamical phenomena giving rise to bursting transcends the existing static classification schemes based solely on slowfast dissection.
Indepth understanding of the generic mechanisms combined into a broad global picture on the transition patterns between the activity types in typical models of individual neurons presents a fundamental challenge for the theory of applied dynamical systems. In response to variations of intrinsic parameters, or an external applied current, like I in (1), a neuron model should demonstrate, migrate, and switch flexibly between various types of activities such as quiescence, tonic spiking and bursting. In addition, nonlinearity of the model can often imply bi or multistability of several coexisting activities at the same parameter values. Bistability of coexisting oscillatory patterns originate near global bifurcations taking place in the model. Multistability is well noticeable when a targeted activity can be robustly selected by choosing other initial conditions or by temporal perturbations, like applied external current. Having ascertained such a global picture we can make consistent predictions for determining basic principles of the functioning of coupled neurons on networks where they receive mixed, inhibitory and excitatory perturbations from other neurons and synergetically reciprocate.
The applied dynamical systems community has developed a universal and versatile suite of computational tools and techniques for comprehensive examination of diverse nonlinear systems of various origins including life sciences. These tools allow a researcher to mono and biparametrically scan the model in question to search for specific transformations corresponding to local and global bifurcations that often hard to detect by standard means. One approach to an initial examination of an unknown model, often called ‘bruteforce’ approach, is to evaluate the largest Lyapunov exponent. This approach has long been widely utilized in nonlinear dynamics for initial detection of bifurcations of steady state and oscillatory solutions. The bruteforce approach contrasts drastically from examination of fine bifurcation structure of limit solutions of the system. Nevertheless, if performed extensively the brute force approach reveals adequately the backbone of the bifurcation structure of the model, which can be further enhanced and complemented with detailed bifurcation analysis that would provide the finishing touches, in the form of bifurcations curves, to the initial brute force diagram. We point out that transformations reshaping bursting waveforms en route toward the tonic spiking activity are caused by atypical bifurcations due to the presence of two, or more time scales in bursting. Because of that bifurcations of stiff bursting solutions, especially irregular ones, are hard to trace down by parameter continuation software packages such as CONTENT and AUTObased packages, which are specifically designed mainly for explorations of equilibria and ‘typical’ periodic orbits.
The primary goal of this paper is to demonstrate that straightforward methods used in neuroscience experimental studies can be as effective as conventional tools, based on the bifurcation and Lyapunov exponent theory, employed in nonlinear dynamics studies. In this paper, we revisit and examine transformations of various oscillatory activity types in the phenomenological HindmarshRose model of bursting neurons, viewed, so to speak, through the prism of neuroscience plausible methods. Next we compare our findings with the results obtained using the evaluation of a maximal Lyapunov exponent that was presented in detail in [27,28] which we consider as an etalon. More specifically, as a part of the comparison test, we place next to each other the bifurcation diagrams found using calculusbased computational tools yielding the whole spectrum of the Lyapunov exponents for complete solutions of the model and those obtained through examinations of 1D voltage traces, which are typically available in experimental studies. We then extract various qualitative temporal characteristic of neuronal activity from nontransient fragments of such traces, including the number of spikes per regular burst, deviations of the spike numbers in case of chaotic bursting, interspike intervals, burst duration and period, and the duty cycle which is the ‘spiking’ fraction of the bursting period. By varying two control parameters of the model, we basically perform biparametric sweeps of its dynamics that are aimed to detect in a very straightforward manner various global bifurcations including • transitions between quiescence, tonic spiking and bursting activities including ones through various homoclinic bifurcations; • identify regular and chaotic transformation of bursting, including a change of the bursting topology accompanying square wave to plateaulike transitioning, as well as forward and backward sequences of spikeaddition and deletion, and so forth.
2 Materials and methods
The phenomenological system of ODEs proposed by Hindmarsh and Rose [29,30] for modeling bursting and spiking oscillatory activities in isolated neurons is given by:
here, x is treated as the membrane potential, while y and z describe some fast and slow gating variables for ionic currents, respectively. Slow ‘activation’ of z is due to the small parameter . The parameters in (1) are typically set as follows and , so that regular bursting oscillations in the model at an ‘applied current’ , which belongs to the squarewave type at , and transforms to a plateaulike bursting at , see Figure 1. Along with ‘intrinsic,’ b, and ‘external,’ I, bifurcation parameters the dynamics of the model are sensitive to variations to other parameters: ε being treated as a rate of activation for some current, and being viewed as a control parameter delaying and advancing the activation of the slow current in the modeled neuron. In response to variations of intrinsic parameters, or an external applied current, like I in (1), a neuron model should demonstrate, migrate, and switch flexibly between various types of activities such as quiescence, tonic spiking and bursting.
In this section we will brief the core of the numerical techniques employed in the analysis of the HR model. We will start with the specifics of the numerical integration of the differential equations of the model (1).
There are a plethora of high quality numerical integrators that have been created by numerical ODEs specialists. This study utilizes a recently developed freeware library TIDES (Taylor Integrator for Differential EquationS) available at http://gme.unizar.es/software/tides webcite[31]. TIDES is a highly adaptive software package for numerical simulations of ODE systems. While the Taylor method is one of the oldest numerical methods for solving ordinary differential equations, it is scarcely used nowadays but its use is growing in the computational dynamics community. The formulation of the method is quite simple. First consider the initial value problem . The value of the solution at can be evaluated as of the nth degree Taylor series of at (f is to be a smooth or analytical function). So, denoting
Therefore, the problem is reduced to the determination of the Taylor coefficients . This can be done efficiently by using automatic differentiation techniques (see details in [32]).
The Taylor method has several unique features [32,33]. One of its features is that it explicitly provides a dense output in the form of a power series, which becomes highly useful for detecting various instantaneous events, for example moments at which a solution hits a Poincaré crosssection, reaches a voltage maximum, for example, if one need to counts spikes in bursts, and so on. In addition, the method can be formulated using the interval arithmetics, which is often employed in the computerassisted proofs of chaos nowadays. The Taylor method also provides highprecision solutions of ODEs so demanded in studies of systems exhibiting multiple timescale dynamics in moderately stiff systems. In this work we use a Taylor series method of order 15 with an error tolerance set to for the most of the simulations. At this tolerance level the TIDES software is as fast and slightly more accurate [31] than the code DOPRI853 developed by Hairer and Wanner [34]. Note that TIDES is a general purpose software, and so it can be apply to general ODE systems, and not only to the HR model. We note also that we have employed the Taylor series method for solving variational equations and computations of the Lyapunov exponent spectrum, which are viewed as nonstandard options for the method [35]. As remarked, in the numerical simulations it is possible to use several good general ODE solvers, but the main advantage of the Taylor method in this kind of studies is that it provides most of the requirements we need for the problem, accuracy when needed, easy events detection, direct dense output as a power series and easy implementation of variational equations.
As mentioned above, the continuous output generated by the integrator based on the Taylor series method is able to detect accurately and effectively various instantaneous events such as whether the phase point hits a crosssection or reaches some critical value like voltage maximum/minimum, or the number of spikes per burst approaches the sought value, which is the underlying idea for the spikecounting (SC) technique [27]. In combination, the methods allows us to classify the solutions of the HR model in neuroscience terms: no spike  quiescence (convergence to a stable equilibrium point); single tonic spike  a round periodic orbit (like one around the manifold ); multiple spikes within a train  bursting orbit composed of alternating tonic spike and pseudoquiescent periods, as well as distinguishing chaotic behaviors characterized by wide variations of spike numbers in burst trains. Besides, the SCtechnique allows for indirect evaluations of the dutycycle [DC] of the orbit, which is a fraction of the burst period (that is, the ratio ) of regular, periodic dynamics. Long bursting implies a dutycycle close to one, meaning that the neuron is active most of the time; on the other hand if DC is close to zero the neuron generates rare spikes, mostly staying in the quasiquiescent state. Finally, the continuous output of the Taylor series is also used in constructing bifurcation diagrams based on the variability of intervals between the spikes.
The other technique that we use in this paper is the computation of the Lyapunov exponents for solutions of the HR model. In bounded systems the existence of a positive Lyapunov exponent is associated with chaotic motions, and any orbit on a compact set that does not converge to an equilibrium point has at least one zero Lyapunov exponent. Accordingly, computation of the largest Lyapunov exponents yields important information about the kind of orbits present in the system. The corresponding algorithm for the computation of the Lyapunov spectrum is a combination of the classical methods proposed by [36,37] and the ability of TIDES to compute solutions of the variational equations directly. The first variational equations of the system are integrated with the identity matrix as an initial condition, creating a canonical orthonormal basis which is mapped onto a new set of vectors. In a chaotic system each vector tends to expand along the local direction of most rapid growth, so to avoid this problem the GramSchmidt orthonormalization process is applied. Later the Lyapunov exponents are calculated from the growth of areas given by the different propagated vectors (for more details see [37,38]).
3 Screening the HR model in the plane
The HR model can exhibit a plethora of dynamical activities at different parameter values. Consequently, obtaining a comprehensive understanding of the multiparametric evolution of a system like the HR model is no easy task, and hence most of the parameters are fixed in the model. In this section, the free bifurcation parameters to be varied are b and I; both are in charge of transformations for the intrinsic structure of the fast subsystem in the HR model. Next we perform the twoparameter sweep or screen of the model to collect vital data representing the time and parameter transformations of the single ‘voltage’ xvariable of the neuronal model. Further data mining will be carried out to extract quantitative and qualitative information about the dynamical variability of the model, bifurcations of its solutions, and so on. In Figures 3 and 4 we use a homogeneous grid comprised of points within the given parameter range. In short, that means that this scan represents 10^{6} simulations.
Fig. 3. (A)parameter sweep of the HindmarshRose model based on the spikecounting approach. The colorcoded bar to the right gives the spikenumber range. The diagram clearly shows the boundaries of the spikeaddition sequence, and the border between squarewave and plateaulike bursting. It also reveals the cloveshape structure of the zones of chaotic bursting which adjoin to the regions of tonicspiking. (B) Samerange screening based on the evaluation of the the dutycycle of bursting. The duty cycle value comes close to one near the boundary between bursting and tonicspiking and drops close to zero near the border of the spiking region. Compare (A) and (B) with the screening diagram based on the Lyapunov exponents for solutions of the model in Figure 4 below.
Fig. 4. Parametric sweeping of a Lyapunov exponent spectrum: orangecolored zones indicate chaotic dynamics in the model, while a regular dynamics region is painted in grey colors of varying tint corresponding to a second Lyapunov exponent of zero value, and darker shades for negative values.
Figure 3(A, B) presents the parameter spikecounting (SC) and dutycycle (DC) diagrams of the HR model. Similar SC diagrams for the model were previously reported in [27]. The colorscale bar on the right in Inset (A) yields the number of spikes within a burst. The diagram for the dutycycle evolution is shown in Inset (B). By combining these diagrams, we can partition the parameter plane into regions of different kinds of behaviors and classify the regimes: tonicspiking (single spike), squarewave and plateaulike bursting, quiescence, and chaotic behaviors with the variability of spikes exceeding some preset limit. It is easy to see that both diagrams give consistent results. They reveal with clarity the region of tonicspiking where both DC and SN take minimum values, below which there is the bursting region at the right bottom corner of the diagram. Bursting emerges from tonic spiking through the spikeaddition cascade in two different ways: one is regular and reversible; the corresponding transitions are foliated by the bifurcation curves. The second kind of transitions is due to the cloveshaped regions (shown in red in Figure 4) corresponding to welldeveloped chaotic dynamics in the model.
Another interesting behavioral phenomenon reshaping the type of bursting occurs in the topleft corner of the diagrams. In this region of bursting with a high number of spikes transforms into bursting with a drastically lower number of spikes per burst. To elucidate what happens on the border between these regions, we visually examined the orbits of the model. We found the border corresponding to the transition between squarewave and plateaulike bursting (see the waveforms in Figure 1 corresponding to the selected points in Figure 3(A)). The corresponding planar bifurcations underlying the robust bursting configurations of both types are well described in the literature, see [6,21]. The bursting type depends on the way the slowmotion tonic spiking manifold terminates in the phase space of the fast subsystem of the HR model. In the squarewave bursting case, the burst termination is due to a homoclinic bifurcation, also known as  fold/homoclinic, whereas in the case of plateaulike, or fold/Hopf, bursting is due to the reverse supercritical AndronovHopf bifurcation (Figures 12). In essence, this means that the parameters b and I change the structure of the fast subsystem of the HR model so that the homoclinic bifurcation is no longer transverse in the zparametric cutaway in the singular limit [21]. Figure 2 elaborates on the metamorphoses of the structure transformations. One can see that plateaulike bursting takes the place of squarewave bursting after the spiking manifold, , becomes tangent to the saddle branch of in the middle of and further terminates on the upper depolarized branch of through the supercritical AndronovHopf bifurcation. In general, this is not a cod1 bifurcation, but a degeneracy due to loss of transversality.
The HR model can exhibit complex, chaotic bursting with a large number of spikes, especially near transitions to hyperpolarized quiescence (see Figure 3). In context of the dynamics of the model, bursting is technically treated as chaotic if trains have more than 25 unlike spikes with distinct interspike intervals. A large number of spikes can also be generated by periodic bursting. In order to differentiate between regular and chaotic bursting behaviors, we have also used another computational technique. The complete spectrum of the Lyapunov exponents was evaluated for the orbits of the model as the two parameters were varied within the same range. In the simulations we discard a transient time of 10^{3} and we integrate till 10^{5} with the algorithm to compute the exponents and using as initial conditions the last value of the previous simulation. The corresponding sweeping diagram is shown in Figure 4. In the diagram, shown in the yelloworange scale refers to the regions where the first Lyapunov exponent is positive. This means the occurrence of chaotic dynamics in the model. The graycolored regions are where the second Lyapunov exponent is negative while the first Lyapunov exponent remain zero on periodic orbits. The Lyapunov exponent based diagram also reveals the spike adding transitions, and the corresponding bifurcation lines can be drawn where the 2nd Lyapunov exponent reaches the maximal value of zero. This implies that the bifurcating bursting orbit is about to disappear and will be replaced by the successive bursting orbit with an extra spike in each train. Spike adding transitions were observed and studied in several neuron models, including the Chay and HindmarshRose mathematical models [15,39,40], and the leech heart interneuron [17]. Note that there are several universal scenarios for such cascades, including saddlenode bifurcations, homoclinic bifurcations of saddle equilibria [19,41] and periodic orbits [18], as well as through the blue sky catastrophe [21,22].
We have pointed out the changes in the color representation (value) of the dutycycle in Figure 3(B). This variability can be due to two reasons: one is the change in the number of spikes per burst; the other reason is due to noticeable changes in time intervals between the spikes (given that the spike duration itself does not vary much). The answer is given by Figure 5 demonstrating the evolution of the temporal characteristics of bursting as the bparameter is moved along the parametric cut . More specifically, using Figure 5 we examine the dependence of the first and second Lyapunov exponents in Inset (A1), the evaluations of the interspike interval bifurcationdiagram in (A2), and the dutycycle of bursting in (A3); finally, the spikecounting diagram in (A4) gives the number of spikes per burst as b is varied. The top of Figure 5 shows a magnification of the SC diagram in Figure 3(A), depicting the regular dynamics of the model, as the Lyapunov exponent stays nonpositive on this pathway, corresponding to some periodic oscillatory activity in the model. The maximal Lyapunov exponent remains at zero, while the second Lyapunov exponent shows intervals of growth to zero alternating with those of rapid decrease. Observe that the peaks in the value of the second Lyapunov exponent occur at spike adding transitions, when the bursting orbit comes close to the saddle which is the threshold between the spikes and the quiescent segment of the former. Spike adding can be detected by the increase of the last interspike interval in the bifurcating bursting orbit that results automatically with an increase of the value of the the duty cycle orbit, and falls right after every spike addition transition. This tendency is clearly interpreted when one examines the changes in the voltage xtraces shown in Figure 6. Inset on the left depicts the evolution of the bursting orbit gaining an extra spike after it comes close by the saddle, and leaves along by the other unstable separatrix of the saddle. One notices from Figure 6 that the first two bursting orbits have six spikes in each burst, while the trace of the second (B) shows a prolonged interspike interval at the end of the burst (also revealed in Inset (A2) in Figure 5). The time interval prior the last spike grows to a point where the isolated spike disappears and is substituted by a short action potential (see the third orbit, where the black point indicates the local xmaxima of the short oscillation). This bursting orbit has five large spikes followed by a single short spike. After this short spike disappears the bursting orbit steadily exhibits five spikes. The process occurs at every spikedeletion bifurcation, or spike adding bifurcation, if instead the parameter b is decreased.
Fig. 5. (A) Magnification of the SC diagram in Figure 3(A). (A1) First and second Lyapunov exponents plotted against the parameter b. Note that the second Lyapunov exponent of the bursting orbit raises to zero at the spikeadding transition and drops after. (A2) Interspike interval vsb at . Last interspike intervals in a burst increase at the spike addition. This is an indication of the homoclinic bifurcation. (A3) Evolution of the dutycycle and (A4) spike variability as b is varied.
Fig. 6. Left: threedimensional [3D] projections of three bursting orbits passing by the saddle. An extra spike is gained by the orbit after it has switched the outgoing directions determined by the unstable separatrices of the saddle. (A)(D) Waveforms of four bursting orbits, shown in the phase space of the model on the left, at spikeadding or deletion bifurcations for the indicated parameter values of I and b.
Indeed, the fact that the interspike interval grows by the end of the burst is a signature of squarewave bursters. Recall that these bursters are also codenamed fold/homoclinic meaning that the spiking, slowmotion manifold is terminated through a homoclinic bifurcation of the saddle that occurs in the fast subsystem of the model, see Figure 2. The dwelling time of the phase point along the bursting orbit grows logarithmically fast the closer the point comes to the saddle [22]. The increase of the dwelling time is a generic phenomenon for all systems in a neighborhood of a saddle. What makes this phenomenon special for slowfast systems is that the time scale of dynamics of the fast subsystem near the saddle turns out to be that of the slow subsystem, which gives rise to another peculiar phenomenon of ‘delayed loss of instability’ such that the phase point, previously turning around the spiking manifold, can be dragged along the middle, saddle branch of equilibria possibly all the way to the upper fold, provided that the timing is right, that is, the phase meets the saddle point right on the edge of the spiking manifold (lousily speaking, we encounter another kind of solutions broadly called ‘canards’, commonly characterized by the fact that a canard can follow an unstable branch of a slowmotion manifold). If the phase point reaches the edge before the saddle, it falls down to the hyperpolarized branch of to start a new cycle of bursting. If the phase point travels past the saddle, then it goes up along the other leading unstable separatrix of the saddle, makes another turn around , resulting in the addition of an extra spike in the bursting orbit. Note that when the phase point does not approach the saddle, the model generates bursts with same number of spikes. Again, let us emphasize that such a spike adding mechanism is typical for squarewave or fold/homoclinic bursters; however, underlying mechanisms for spike adding can be completely different even in square wave bursters, and other neuronal models [19,24,40,42], including the leech heart interneuron model [17,18].
4 Screening the HRmodel in the and planes
In this section we examine the dynamics of the model in response to variations of the slow parameter ε. In the neuroscience context, ε can be treated as the reciprocal of τ, which determines the (in)activation rate of the slow current in a neuronal model. For the sake of consistency we will first screen the model in the parameter plane, while fixing and . Recall that the parameter moves the slow nullcline of the model up and down in the xdirection (see Figure 7). As the slow equation in (1) contains no yvariable, the plane in the phase space of the HR model, where the time derivative vanishes, is the slow nullcline. One can see that and below and above this nullcline, respectively. Note that a simple round periodic orbit on the tonic spiking manifold, , corresponds to regular tonic spiking activity in the model. The position of the periodic orbit on depends on where the slow nullcline cuts through . By changing , we make the periodic orbit shift along the spiking manifold. More specifically it can be found around the intersection points of the slow nullcline with the average branch , see details in [2]. Tonic spiking remains regular until the periodic orbit stays away from the ‘homoclinic’ edge of .
Fig. 7. (3D version of Figure 2) Intersection point of the branch with the slow nullcline yields the equilibrium state of the HR model at given . The dark blue point is the center of the gravity of the stable periodic orbit of the HR model, which is depicted on the tonic spiking manifold at . It is located around the intersection point of the slow nullcline with the space curve of the average xvalues on each periodic orbit that foliate the spiking manifold . The phase point, while turning around , moves slowly toward the homoclinic edge () as long as it stays above the slow nullcline, and goes backward () after it lowers below the slow nullcline . When these opposite forces are canceled out over the revolution period, the phase point spins around the ‘center of the gravity’, that is, stays on the same periodic orbit. The variations of move the slow nullcline and thus make the periodic orbit slide along the manifold. When the slow nullcline cuts through the unstable section of between HB, standing for homoclinic bifurcation in the fast subsystem and the fold AH/SN, the model becomes a burster.
As it was said previously the HR model describes one of the most typical configurations of slow manifolds for squarewave bursting oscillations. First of all, the configuration needs the distinct Zshape for the quiescent manifold, , with the lower branch corresponding to a hyperpolarized quiescent state of the neuron, while the upper unstable branch is surrounded by the spiking manifold, , foliated by the stable limit cycles of the fast subsystem in the square bursting case. The branch regains stability in the case of plateaulike bursting. The manifold terminates through the homoclinic bifurcation that occurs in the fast subsystem in the square bursting case. Between the lower fold and this homoclinic point, the system has a hysteresis which gives rise to bursting. In the bursting regime, the phase point of the HR model switches repeatedly between the spiking, , and quiescent, , manifolds when it reaches their ends. In addition, both manifolds must be transient for the passing solutions of (1), that is, must be cut by the slow nullcline through the middle, saddle branch below and above the hyperpolarized fold point. This guarantees that is also transient for the trajectories of the model that coil around while translating slowly towards the edge, which corresponds to the aforementioned homoclinic bifurcation. Thus, the rapid jump from the lower point on towards indicates the beginning of the spiking period of a burst followed by the resting phase when the phase point drifts slowly along towards the fold, onto which it lands right after the homoclinic bifurcation. The number of complete revolutions of the phase point around the spiking manifold, , gives the number of spikes within a burst, see Figures 2 and 6. Bursting in the model takes place as long as the slow nullcline hits between the points labeled HB, corresponding to the homoclinic bifurcation, and AH/SN standing for the singular AndronovHopf bifurcation in the full system, see Figure 7.
Thus, by varying we can make the model generate trains of bursts with various number of spikes. It is easy to see that the value of the small parameter ε determines the slow passage along both manifolds. So, halving ε should make bursts twice as long at least, with a doubled number of spikes. Note, too, that the duration of the quiescent phase should increase proportionally. As for variations of I are concerned, I moves, geometrically, the manifold horizontally in the 3D phase space of the model, in particular due to linearity of the slow equation in both variables. We show that because of that both I and act similarly on dynamics; of special interest here are transitions between the activity type of the model, which this section is focused on.
Figure 8 demonstrates the 2D spike counting diagram with and the same picture using the first and second Lyapunov exponents. The diagram reveals a diagonal plot structure foliated by homogeneous bands. This suggests that variations of parameters I and cause similar responses in the model. On the band structure there is a thin band of chaotic motion located inside the band of high number of spikes, as remarked by the positive value of the maximum Lyapunov exponent. To gain insight into the band structure, we examine the evolutions of dynamics of the model as only is varied at fixed current input . At smaller values of the parameter, the model exhibits tonic spiking first, see Insets (C) and (E) of Figure 8 presenting the interspike bifurcation diagram (IBC), and the spike counting (SC) diagram. As the parameter is decreased further toward the bursting zone, the model enters a perioddoubling cascade leading to chaos [20,43]. In Figure 8(E) we see groups of spikes but without a clear bursting structure (bursting orbits with n spikes that doubles the period are denoted by ). Next, the model has a chaotic orbit (in region Ch), evolving into a compact chaotic region, which undergoes a boundary crisis, widening drastically the size of the chaotic attractor (also reported in [44,45]). As is decreased further, chaos terminates through another boundary crisis due to intermittency originating from a fold bifurcation. The model now generates regular trains of bursts. The corresponding bursting orbit further undergoes a series of perioddoubling and periodhalving bifurcations, before it steps into a cascade of spike deletion bifurcations, eventually leading to quiescence (Q) on the left hand side of of the parametric pathway, after the AndronovHopf bifurcation.
Fig. 8. (A) 2D Lyapunov exponents diagram at ; (B) 2D spikecounting diagram at ; (C) 1D interspikeinterval diagram vs for ; (D) Bursting dutycycle dependence on ; (E) Spike variability per burst is plotted against .
Figure 9 demonstrates the waveforms of the bursting solutions at selected points (shown in magenta) in Inset (C) of Figure 8. One complete period of each waveform is squared out within a red box. Thus, we have the following orbits: a spiking orbit for , two bistability points for and −0.9909, orbit (bursting with two spikes) for orbit for , chaotic orbit for and again a regular spiking orbit at of a smaller period compared to that at . A quite interesting observation can be deduced from Figure 8(C). Namely, there are several regions of bistability where there are coexisting stable periodic orbits. These regions, labeled as CR1, CR2 and CR3, are marked by using two colors in the IBDdiagram in Inset (C). The middle bistability region, CR2, is expanded in Figure 10(A), and it reveals that there coexist two distinct bursting orbits at . The 3D phase projections and the corresponding waveforms, with single spikes and spike duplets (B(2)), are shown in Figures 10(B) and 9, respectively. Figure 9 presents two plots of the waveforms for the coexisting orbits at ; the black point in Figure 6(C), indicates the local xmaxima of the short waveform. As was pointed out earlier such a coexistence of spike and spike bursts is a typical phenomenon for squarewave bursters at the spike adding or deletion transitions due to ‘the delay loss of instability’ that occurs along the saddle, threshold branch of that separates depolarized and hyperpolarized states of the neuron model [41]. Recall that theoretically, due to the equal time scale dynamics of the fast subsystem near the saddle and the slow equation because of small ε, the phase point can be dragged along the saddle branch up to the upper fold on . Interestingly enough, the range of bistability zones will shrink as the value of ε is increased. We would like to point out that multistability is a byproduct of the nonlinearity in the system. Multistability has been reported in several neuronal systems both experimentally and computationally, including individual neurons and their models, as well as in neural networks and multifunctional central pattern generators [1,46]. Multistability is of great interest for neuroscience as it can potentially enhance the flexibility of the nervous systems, decision making processes, and explain various nervous pathologies caused by sudden changes in system’s states.
Fig. 9. Tonic spiking and bursting xtraces at selected points on the bparameter passway (shown in magenta) across the diagram in Figure 8(C) showing several stages of spike addition transitions.
Fig. 10. (A) Magnification of the interspike bifurcationdiagram (IBD) for in the bistable region CR2 from Figure 8. It reveals two coexisting bursting orbits (shown in red and blue) at in the phase space of the model: the split between the solutions takes place near the saddle after the orbits departs in the direction determined by the unstable separatrices. (B) Two coexisting orbits corresponding to single and duplet spiking.
Another peculiar observation related to the band structure in the parametric plane is concerned with the regions of high sensitivity to small variations of the control parameter , whereas the overall band structure seems to be quite robust, or selfsimilar in I. This brings in one remaining question we would like to addresses in the paper: will this property of selfsimilarity of the band structure persist for smaller values of ε? Our findings are summarized in Figure 11 which demonstrates the bifurcation diagrams for two values: and . They confirm that the band structure does persist and show predictably regions with larger number of spikes per bursts, especially for (compare the SC plot in 1D and 2D diagrams). To study the genesis of the band structure we also built the corresponding diagrams shown in Figure 12. Although both diagrams represent the same data of the spikecounting, the diagram on the right in Figure 12 is given in a logarithmic scale to demonstrate that the selfsimilarity property of the band structure is exponential in ε. We can deduce from the diagrams that the most interesting, in terms of dynamics, region is contained between two curves (indicated by white dots). Within this region there are several diagonal bands corresponding to bursting orbits with different numbers of spikes. This plot explains the fact that a small εparametric cut will nevertheless reveal the band structure, observed in Figures 8 and 11, thus confirming that either bifurcation parameter I or can be equality singled out to perform the bifurcation analysis of the HindmarshRose model.
Fig. 11. (A) 2D spikecounting diagram projected to the plane at ; (A1) interspikeduration bifurcation diagram for . (B) 2D spikecounting plot in the plane for ; (B1) interspikeduration bifurcation diagram for and (B2) 1D spikecounting diagram for .
Fig. 12. spikecounting diagram in the linear and logarithmic scales in ε. Shown in blue and red are regions of passive silence and intense spiking activities.
5 Conclusions
As of today, the HindmarshRose remains justifiably one of the most popular mathematical models, that describes, qualitatively well, the dynamics of a certain class of neuronal models derived using the HodgkinHuxley formalism. The model has been carefully analyzed using various mathematical and computational tools, and has been viewed though the prisms of the advanced bifurcation theory, geometric methods of slowfast dynamical systems to reveal multiple peculiar qualitative features. Various brute force approaches [27,28] have been applied to the model to reveal its quantitative or metric properties though the examination of the Lyapunov exponent spectrum and the number of spikes per period.
In this paper we have tested other computationally effective tools tailored specifically for models originating in neuroscience, including 1D and 2D parametric screenings of the complex dynamics of the model aimed to bridge notions common for neuroscience with the accurate mathematicsbased findings. Our computational toolkit includes several methods for examining of temporal characteristics of a single variable, x, treated as a voltage across the cell membrane, or more specifically, the corresponding voltage waveforms. The list includes the spikecounting approach, the evaluations of interspike interval, and duty cycles of bursting, which is a ratio of the burst duration (active phase) over the total burst period. We confirmed our findings based on these methods with ‘calculus’based simulations for the whole spectrum of the Lyapunov exponents calculated for solutions of a system of ordinary differential equations like the HindmarshRose model. Our verdict is that both approaches consistently demonstrate very good agreements. The methods allow us to give detailed explanations for various global bifurcations in the model, including various spikeaddition/deletion cascades, phenomena of bistability, and various transitions between types of bursting: squarewave and plateaulike. This ensures that this toolkit of combined neurosciencenative methods and calculibased algorithms will yield effective and timely insights for pilot studies of new, previously unidentified, in sense of dynamics and bifurcations, models of individual neurons, as well as other cells, such as myocytes  cardiac tissue cell. We showed in [26] that the proposed methods could provide the bifurcation details and add even problemspecific nuances into regulation control of temporal characteristics of realistic interneuron models of the leech. The evident advantage of the approach is its nativeness for the neuroscience community. Last but not least, it should emphasized that a 2D parameter sweeping of the model shown in Figures 3 and 4 takes about 1020 times faster than the sweeping based on the Lyapunov exponents spectrum (Figure 5). A drawback of the approach is that it needs to be corrected in cases where the model is multistable; this remains common weakness of all straightforward methods unless one uses randomized initial conditions and longer transients that could overall substantially prolong the simulation time.
It is in our future plans to broaden the applicability of these proposed computational tools for studies of neuronal networks, especially for multifunctional central pattern generators comprised of several neurons [2,46]. Such a multifunctional CPG is capable of generating multiple bursting rhythms at quite different time scale. It was shown recently in [3] that bursting outcomes of a multistable 3cell network are determined by the duty cycle of bursting. Moreover, the longest bursting cell plays the role of the pacemaker of the network [47]. Note also that bursting network can be composed of individually tonic spiking cells that while inhibiting, even weekly, each other can create various bursting outcomes of the network as a whole. This indicates directly that such a cell, whether tonic spiking or bursting, is to be close to the boundary separating the activity types in the parameter space of distinct interneurons [25,26], in particular to control effectively the temporal characteristic of bursting, regular or chaotic. In the light of saying, it is evident that the tools native to neuroscience paradigms are suited more appropriately for studying bursting metamorphoses in networks and the sidebyside comparison of the results of mathematical and experimental studies using the common jargon. This computational toolkit shall bring us closer to the targeted goal  to build realistic and adequately responding models of concrete functional CPGs with specific time scales, phase locked states between synergetically coupled neurons with plausible bursting characteristics.
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
We would like to thank A. Neiman, J. Wojcik and B. Chung for very helpful comments. This work is supported by the Spanish Research project MTM200910767, and by NSF Grant DMS1009591, RFFI Grant No. 080100083, GSU Brain & Behaviors program, and MESRF ‘Attracting leading scientists to Russian universities’ project 14.740.11.0919.
References

Kristan, WB: Neuronal decisionmaking circuits. Curr. Biol.. 18(19), 928–932 (2008)

Shilnikov, AL, Gordon, R, Belykh, I: Polyrhythmic synchronization in bursting networking motifs. Chaos. 18(3), (2008)

Wojcik, J, Clewley, R, Shilnikov, AL: Order parameter for bursting polyrhythms in multifunctional central pattern generators. Phys. Rev. E. 83(5), (2011)

Rinzel, J: Bursting oscillations in an excitable membrane model. Lect. Notes Math.. 1151, 304–316 (1985)

Bertram, R, Butte, MJ, Kiemel, T, Sherman, A: Topological and phenomenological classification of bursting oscillations. Bull. Math. Biol.. 57(3), 413–439 (1995)

Izhikevich, EM: Dynamical Systems in Neuroscience. The Geometry of Excitability and Bursting, MIT Press, Cambridge, MA (2007)

Tikhonov, AN: On the dependence of solutions of differential equations from a small parameter. Mat. Sb.. 22(64), 193–204 (1948)

Pontryagin, LS, Rodygin, LV: Periodic solution of a system of ordinary differential equations with a small parameter in the terms containing derivatives. Sov. Math. Dokl.. 1, 611–619 (1960)

Fenichel, F: Geometric singular perturbation theory for ordinary differential equations. J. Differ. Equ.. 31, 53–98 (1979)

Mischenko, E.F., Rozov, N.K.: Differential Equations with Small Parameters and Relaxation Oscillations. Plenum Press (1980)

Arnold, V.I., Afraimovich, V.S., Ilyashenko, Y.S., Shilnikov, L.P.: Bifurcation Theory, Dynamical Systems. Volume V. Encyclopaedia of Mathematical Sciences. Springer (1994)

Terman, D: The transition from bursting to continuous spiking in excitable membrane models. J. Nonlinear Sci.. 2(2), 135–182 (1992)

Shilnikov, AL, Calabrese, R, Cymbalyuk, G: Mechanism of bistability: tonic spiking and bursting in a neuron model. Phys. Rev. E. 71(5), (2005)

Holden, AV, Fan, YS: From simple to simple bursting oscillatory behaviour via intermittent chaos in the RoseHindmarsh model for neuronal activity. Chaos Solitons Fractals. 2, 349–369 (1992)

Fan, YS, Holden, AV: Bifurcations, burstings, chaos and crises in the RoseHindmarsh model for neuronal activity. Chaos Solitons Fractals. 3, 439–449 (1995)

Shilnikov, AL, Cymbalyuk, G: Transition between tonic spiking and bursting in a neuron model via the blue sky catastrophe. Phys. Rev. Lett.. 94(4), (2005)

Channell, P, Cymbalyuk, G, Shilnikov, AL: Origin of bursting through homoclinic spike adding in a neuron model. Phys. Rev. Lett.. 98(13), (2007)

Channell, P, Fuwape, I, Neiman, AB, Shilnikov, AL: Variability of bursting patterns in a neuron model in the presence of noise. J. Comput. Neurosci.. 27(3), 527–542 (2009)

Innocenti, G, Genesio, R: On the dynamics of chaotic spikingbursting transition in the HindmarshRose neuron. Chaos. 19(2), (2009)

Wang, XJ: Genesis of bursting oscillations in the HindmarshRose model and homoclinicity to a chaotic saddle. Physica D. 62(14), 263–274 (1993)

Shilnikov, AL, Kolomiets, ML: Methods of the qualitative theory for the HindmarshRose model: a case study. A tutorial. Int. J. Bifurc. Chaos Appl. Sci. Eng.. 18(8), 2141–2168 (2008)

Shilnikov, L.P., Shilnikov, A.L., Turaev, D.V., Chua, L.O.: Methods of Qualitative Theory in Nonlinear Dynamics. Parts I and II. World Scientific Publishing Co. Inc. (1998) and (2001)

Feudel, U, Neiman, A, Pei, X, Wojtenek, W, Braun, H, Huber, M, Moss, F: Homoclinic bifurcation in a HodgkinHuxley model of thermally sensitive neurons. Chaos. 10(1), 231–239 (2000)

Kramer, MA, Traub, RD, Kopell, NJ: New dynamics in cerebellar purkinje cells: torus canards. Phys. Rev. Lett.. 101(6), (2008)

Wojcik, J, Shilnikov, AL: Voltage interval mappings for dynamics transitions in elliptic bursters. Physica D. 240, 1164–1180 (2011)

Shilnikov, AL: Complete dynamical analysis of an interneuron model. Nonlinear Dyn. (2011)

Storace, M, Linaro, D, de Lange, E: The HindmarshRose neuron model: bifurcation analysis and piecewiselinear approximations. Chaos. 18(3), (2008)

de Lange, E, Hasler, M: Predicting single spikes and spike patterns with the HindmarshRose model. Biol. Cybern.. 99(45), 349–360 (2008)

Hindmarsh, JL, Cornelius, P: The development of the HindmarshRose model for bursting. Bursting (2005)

Hindmarsh, JL, Rose, RM: A model of the nerve impulse using three coupled firstorder differential equations. Proc. R. Soc. Lond. B, Biol. Sci.. 221, 87–102 (1984)

Abad, R., Barrio, R., Blesa, F., Rodríguez, M.: TIDES: a Taylor series Integrator for Differential EquationS. ACM T. Math Software, (2011 in press). http://gme.unizar.es/software/tides

Barrio, R: Sensitivity analysis of ODEs/DAEs using the Taylor series method. SIAM J. Sci. Comput.. 27(6), 1929–1947 (2006)

Barrio, R, Blesa, F, Lara, M: VSVO formulation of the Taylor method for the numerical solution of ODEs. Comput. Math. Appl.. 50(12), 93–111 (2005)

Hairer, E, Nørsett, SP, Wanner, G: Solving Ordinary Differential Equations. I, SpringerVerlag, Berlin (1993)

Barrio, R, Rodríguez, M, Abad, A, Blesa, F: Breaking the limits: the Taylor series method. Appl. Math. Comput.. 167(20), 7940–7954 (2011)

Benettin, G, Galgani, L, Giorgilli, A, Strelcyn, JM: Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 2: Numerical application. Meccanica. 15, 21–30 (1980)

Wolf, A, Swift, JB, Swinney, HL, Vastano, JA: Determining Lyapunov exponents from a time series. Physica D. 16(3), 285–317 (1985)

Skokos, C: The Lyapunov characteristic exponents and their computation. Lect. Notes Phys.. 790, 63–135 (2010)

Chay, TR: Chaos in a threevariable model of an excitable cell. Physica D. 16(2), 233–242 (1985)

Yang, Z, Qishao, L, Li, L: The genesis of periodadding bursting without burstingchaos in the Chay model. Chaos Solitons Fractals. 27(3), 689–697 (2006)

Mozekilde, E, Lading, B, Yanchuk, S, Maistrenko, Y: Bifurcation structure of a model of bursting pancreatic cells. Biosystems. 63, 2–13 (2001)

Innocenti, G, Morelli, A, Genesio, R, Torcini, A: Dynamical phases of the HindmarshRose neuronal model: studies of the transition from bursting to spiking chaos. Chaos. 17(4), (2007)

Cymbalyuk, G, Shilnikov, AL: Coexistence of tonic spiking oscillations in a leech neuron model. J. Comput. Neurosci.. 18(3), 255–263 (2005)

GonzálezMiranda, JM: Observation of a continuous interior crisis in the HindmarshRose neuron model. Chaos. 13(3), 845–852 (2003)

GonzálezMiranda, JM: Complex bifurcation structures in the HindmarshRose neuron model. Int. J. Bifurc. Chaos Appl. Sci. Eng.. 17(9), 3071–3083 (2007)

Briggman, KL, Kristan, WB: Multifunctional patterngenerating circuits. Annu. Rev. Neurosci.. 31, 271–294 (2008)

Belykh, IV, Shilnikov, AL: When weak inhibition synchronizes strongly desynchronizing networks of bursting neurons. Phys. Rev. Lett.. 101, (2008)