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
Fig. 2. Top: a plateaulike or fold/Hopf bursting starts after the spiking manifold
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
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
Therefore, the problem is reduced to the determination of the Taylor coefficients
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
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
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
(
b
,
I
)
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
Fig. 3. (A)
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
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,
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
Fig. 5. (A) Magnification of the
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
4 Screening the HRmodel in the
(
x
0
,
I
)
and
(
x
0
,
ε
)
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
Fig. 7. (3D version of Figure 2) Intersection point of the branch
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,
Thus, by varying
Figure 8 demonstrates the 2D
Fig. 8. (A) 2D
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
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
Another peculiar observation related to the band structure in the
Fig. 11. (A) 2D spikecounting diagram projected to the plane
Fig. 12.
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)