Open Access Highly Accessed Research

Parameter-sweeping techniques for temporal dynamics of neuronal systems: case study of Hindmarsh-Rose model

Roberto Barrio1 and Andrey Shilnikov2*

Author Affiliations

1 Departamento de Matemática Aplicada and IUMA, Universidad de Zaragoza, E-50009, Zaragoza, Spain

2 Neuroscience Institute and Department of Mathematics and Statistics, Georgia State University, Atlanta, Georgia, 30303, USA

For all author emails, please log on.

The Journal of Mathematical Neuroscience 2011, 1:6  doi:10.1186/2190-8567-1-6


The electronic version of this article is the complete one and can be found online at: http://www.mathematical-neuroscience.com/content/1/1/6


Received:1 February 2011
Accepted:11 July 2011
Published:11 July 2011

© 2011 Barrio, Shilnikov; licensee Springer

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 two-parameter screening of dynamics in neuronal models. We test a ‘brute-force’ 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 Hindmarsh-Rose model of a bursting neuron and compare the results obtained by calculus-based 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 two-time scale bursting, is done by revealing generic properties of mathematical and realistic models of neurons; the latter are derived through the Hodgkin-Huxley formalism for gating variables. Either bursting model falls into a class of dynamical systems with at least two time scales, known as slow-fast 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 so-called slow motion manifolds composed of the limiting solutions, such as equilibria and limit cycles, of the fast subsystem of a model [7-11]. These manifolds constitute the backbones of bursting patterns in a neuronal model. A typical Hodgkin-Huxley model possesses a pair of such manifolds [4]: quiescent and tonic spiking. The existing classifications of bursting are based on codimension-one bifurcations that initiate or terminate the fast trajectory transitions between such one-dimensional [1D] and two-dimensional [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 Hopf-fold; square-wave burster, or fold-homoclinic; parabolic, or circle-circle class describing top-hat 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 Hindmarsh-Rose 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 Andronov-Hopf bifurcation, along with a saddle-node 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.

thumbnailFig. 1. (A) Square-wave and (B) plateau-like bursting traces in the Hindmarsh-Rose model at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M1">View MathML</a>, respectively. Transformations of bursting can be detected quantitatively by a sudden change of the number of spikes per burst.

thumbnailFig. 2. Top: a plateau-like or fold/Hopf bursting starts after the spiking manifold <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> becomes tangent to the middle, saddle branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> and terminates further through the reverse supercritical Andronov-Hopf bifurcation on the upper depolarized branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a>. Bottom: the primary feature of the square-wave bursting activity in the HR model also referred to as of fold/homoclinic type is the termination of the spiking manifold <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> by the homoclinic bifurcation in the phase space of the fast subsystem. In both cases: fold stands for a saddle-node bifurcation at the turning point (SN) on the lower, hyperpolarized branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a>.

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 slow-fast 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 slow-fast dissection has been proven to work very well for a low-order 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 saddle-node 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 square-wave burster like the Hindmarsh-Rose model [13-15]. In addition [12] gives an explanation of common spike adding cascade in classical square-wave 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 codimension-two homoclinic bifurcations, including inclination-switch and orbit-flips, 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 slow-fast dynamics include various homoclinic of saddles and saddle-foci, the blue sky catastrophe, bistability due non-transverse homoclinics to saddle-node periodic orbits, canard-tori [13,16-18,23-26]. The range of bifurcations and dynamical phenomena giving rise to bursting transcends the existing static classification schemes based solely on slow-fast dissection.

In-depth 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 multi-stability of several co-existing 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 bi-parametrically 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 ‘brute-force’ 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 brute-force 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 AUTO-based 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 Hindmarsh-Rose 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 calculus-based 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 non-transient 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 bi-parametric 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 plateau-like transitioning, as well as forward and backward sequences of spike-addition 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:

<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M7">View MathML</a>

(1)

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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M8">View MathML</a>. The parameters in (1) are typically set as follows <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M9">View MathML</a><a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M10">View MathML</a><a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M11">View MathML</a><a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M12">View MathML</a><a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M13">View MathML</a> and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M14">View MathML</a>, so that regular bursting oscillations in the model at an ‘applied current’ <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M15">View MathML</a>, which belongs to the square-wave type at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M16">View MathML</a>, and transforms to a plateau-like bursting at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M17">View MathML</a>, 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M19">View MathML</a>. The value <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M20">View MathML</a> of the solution at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M21">View MathML</a> can be evaluated as <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M22">View MathML</a> of the nth degree Taylor series of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M23">View MathML</a> at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M24">View MathML</a> (f is to be a smooth or analytical function). So, denoting <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M25">View MathML</a>

<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M26">View MathML</a>

Therefore, the problem is reduced to the determination of the Taylor coefficients <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M27">View MathML</a>. 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é cross-section, 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 computer-assisted proofs of chaos nowadays. The Taylor method also provides high-precision solutions of ODEs so demanded in studies of systems exhibiting multiple time-scale dynamics in moderately stiff systems. In this work we use a Taylor series method of order 15 with an error tolerance set to <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M28">View MathML</a> 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 cross-section 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 spike-counting (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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>); multiple spikes within a train - bursting orbit composed of alternating tonic spike and pseudo-quiescent periods, as well as distinguishing chaotic behaviors characterized by wide variations of spike numbers in burst trains. Besides, the SC-technique allows for indirect evaluations of the duty-cycle [DC] of the orbit, which is a fraction of the burst period (that is, the ratio <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M30">View MathML</a>) of regular, periodic dynamics. Long bursting implies a duty-cycle 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 quasi-quiescent 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 Gram-Schmidt 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31">View MathML</a>-plane

The HR model can exhibit a plethora of dynamical activities at different parameter values. Consequently, obtaining a comprehensive understanding of the multi-parametric 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 two-parameter sweep or screen of the model to collect vital data representing the time and parameter transformations of the single ‘voltage’ x-variable 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M32">View MathML</a> points within the given parameter range. In short, that means that this scan represents 106 simulations.

thumbnailFig. 3. (A)<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M33">View MathML</a>-parameter sweep of the Hindmarsh-Rose model based on the spike-counting approach. The color-coded bar to the right gives the spike-number range. The diagram clearly shows the boundaries of the spike-addition sequence, and the border between square-wave and plateau-like bursting. It also reveals the clove-shape structure of the zones of chaotic bursting which adjoin to the regions of tonic-spiking. (B) Same-range screening based on the evaluation of the the duty-cycle of bursting. The duty cycle value comes close to one near the boundary between bursting and tonic-spiking 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.

thumbnailFig. 4. Parametric sweeping of a Lyapunov exponent spectrum: orange-colored 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31">View MathML</a>-parameter spike-counting (SC) and duty-cycle (DC) diagrams of the HR model. Similar SC diagrams for the model were previously reported in [27]. The color-scale bar on the right in Inset (A) yields the number of spikes within a burst. The diagram for the duty-cycle 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: tonic-spiking (single spike), square-wave and plateau-like 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 tonic-spiking 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 spike-addition 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 clove-shaped regions (shown in red in Figure 4) corresponding to well-developed chaotic dynamics in the model.

Another interesting behavioral phenomenon reshaping the type of bursting occurs in the top-left 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 square-wave and plateau-like 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 slow-motion tonic spiking manifold terminates in the phase space of the fast subsystem of the HR model. In the square-wave bursting case, the burst termination is due to a homoclinic bifurcation, also known as - fold/homoclinic, whereas in the case of plateau-like, or fold/Hopf, bursting is due to the reverse supercritical Andronov-Hopf 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 z-parametric cutaway in the singular limit [21]. Figure 2 elaborates on the metamorphoses of the structure transformations. One can see that plateau-like bursting takes the place of square-wave bursting after the spiking manifold, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, becomes tangent to the saddle branch of in the middle of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> and further terminates on the upper depolarized branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> through the supercritical Andronov-Hopf bifurcation. In general, this is not a cod-1 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 103 and we integrate till 105 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 yellow-orange scale refers to the regions where the first Lyapunov exponent is positive. This means the occurrence of chaotic dynamics in the model. The gray-colored 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 Hindmarsh-Rose mathematical models [15,39,40], and the leech heart interneuron [17]. Note that there are several universal scenarios for such cascades, including saddle-node 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 duty-cycle 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 b-parameter is moved along the parametric cut <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M38">View MathML</a>. 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 bifurcation-diagram in (A2), and the duty-cycle of bursting in (A3); finally, the spike-counting 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 non-positive 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 x-traces 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 x-maxima 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 spike-deletion bifurcation, or spike adding bifurcation, if instead the parameter b is decreased.

thumbnailFig. 5. (A) Magnification of the <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M31">View MathML</a>-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 spike-adding transition and drops after. (A2) Interspike interval vsb at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M38">View MathML</a>. Last interspike intervals in a burst increase at the spike addition. This is an indication of the homoclinic bifurcation. (A3) Evolution of the duty-cycle and (A4) spike variability as b is varied.

thumbnailFig. 6. Left: three-dimensional [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 spike-adding 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 square-wave bursters. Recall that these bursters are also code-named fold/homoclinic meaning that the spiking, slow-motion 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 slow-fast 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> 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 slow-motion manifold). If the phase point reaches the edge before the saddle, it falls down to the hyperpolarized branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, 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 square-wave 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 HR-model in the <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a> and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M45">View MathML</a>-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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a>-parameter plane, while fixing <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M47">View MathML</a> and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M48">View MathML</a>. Recall that the parameter <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> moves the slow nullcline of the model up and down in the x-direction (see Figure 7). As the slow equation in (1) contains no y-variable, the plane in the <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M50">View MathML</a>-phase space of the HR model, where the time derivative <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M51">View MathML</a> vanishes, is the slow nullcline. One can see that <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M52">View MathML</a> and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M53">View MathML</a> below and above this nullcline, respectively. Note that a simple round periodic orbit on the tonic spiking manifold, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, corresponds to regular tonic spiking activity in the model. The position of the periodic orbit on <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> depends on where the slow nullcline <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56">View MathML</a> cuts through <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>. By changing <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>, 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M59">View MathML</a>, see details in [2]. Tonic spiking remains regular until the periodic orbit stays away from the ‘homoclinic’ edge of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>.

thumbnailFig. 7. (3D version of Figure 2) Intersection point of the branch <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> with the slow nullcline <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M62">View MathML</a> yields the equilibrium state of the HR model at given <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>. 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M65">View MathML</a>. It is located around the intersection point of the slow nullcline <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56">View MathML</a> with the space curve <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M59">View MathML</a> of the average x-values on each periodic orbit that foliate the spiking manifold <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>. The phase point, while turning around <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, moves slowly toward the homoclinic edge (<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M53">View MathML</a>) as long as it stays above the slow nullcline, and goes backward (<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M52">View MathML</a>) after it lowers below the slow nullcline <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56">View MathML</a>. 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> move the slow nullcline and thus make the periodic orbit slide along the manifold. When the slow nullcline <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M56">View MathML</a> cuts through the unstable section of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> 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 square-wave bursting oscillations. First of all, the configuration needs the distinct Z-shape for the quiescent manifold, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a>, with the lower branch corresponding to a hyperpolarized quiescent state of the neuron, while the upper unstable branch is surrounded by the spiking manifold, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, foliated by the stable limit cycles of the fast subsystem in the square bursting case. The branch regains stability in the case of plateau-like bursting. The manifold <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> 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, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, and quiescent, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a>, manifolds when it reaches their ends. In addition, both manifolds must be transient for the passing solutions of (1), that is, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> must be cut by the slow nullcline through the middle, saddle branch below <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> and above the hyperpolarized fold point. This guarantees that <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> is also transient for the trajectories of the model that coil around <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> while translating slowly towards the edge, which corresponds to the aforementioned homoclinic bifurcation. Thus, the rapid jump from the lower point on <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> towards <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a> indicates the beginning of the spiking period of a burst followed by the resting phase when the phase point drifts slowly along <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> 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, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M2">View MathML</a>, 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> between the points labeled HB, corresponding to the homoclinic bifurcation, and AH/SN standing for the singular Andronov-Hopf bifurcation in the full system, see Figure 7.

Thus, by varying <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a>-spike counting diagram with <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> is varied at fixed current input <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a>. 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 period-doubling 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M97">View MathML</a>). 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> 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 period-doubling and period-halving bifurcations, before it steps into a cascade of spike deletion bifurcations, eventually leading to quiescence (Q) on the left hand side of of the <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>-parametric pathway, after the Andronov-Hopf bifurcation.

thumbnailFig. 8. (A) 2D <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a> Lyapunov exponents diagram at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93">View MathML</a>; (B) 2D <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a> spike-counting diagram at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M93">View MathML</a>; (C) 1D interspike-interval diagram vs<a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a>; (D) Bursting duty-cycle dependence on <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>; (E) Spike variability per burst is plotted against <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>.

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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M108">View MathML</a>, two bistability points for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M109">View MathML</a> and −0.9909, <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M110">View MathML</a>-orbit (bursting with two spikes) for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M111">View MathML</a><a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M112">View MathML</a>-orbit for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M113">View MathML</a>, chaotic orbit for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M114">View MathML</a> and again a regular spiking orbit at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M115">View MathML</a> of a smaller period compared to that at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M108">View MathML</a>. 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 IBD-diagram in Inset (C). The middle bistability region, CR2, is expanded in Figure 10(A), and it reveals that there co-exist two distinct bursting orbits at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M117">View MathML</a>. 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 co-existing orbits at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M118">View MathML</a>; the black point in Figure 6(C), indicates the local x-maxima of the short waveform. As was pointed out earlier such a coexistence of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M119">View MathML</a>-spike and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M120">View MathML</a>-spike bursts is a typical phenomenon for square-wave bursters at the spike adding or deletion transitions due to ‘the delay loss of instability’ that occurs along the saddle, threshold branch of <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a> that separates depolarized and hyper-polarized 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M3">View MathML</a>. 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 by-product 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.

thumbnailFig. 9. Tonic spiking and bursting x-traces at selected points on the b-parameter passway (shown in magenta) across the diagram in Figure 8(C) showing several stages of spike addition transitions.

thumbnailFig. 10. (A) Magnification of the interspike bifurcation-diagram (IBD) for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a> in the bistable region CR2 from Figure 8. It reveals two coexisting bursting orbits (shown in red and blue) at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M124">View MathML</a> 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a>-parametric plane is concerned with the regions of high sensitivity to small variations of the control parameter <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a>, whereas the overall band structure seems to be quite robust, or self-similar in I. This brings in one remaining question we would like to addresses in the paper: will this property of self-similarity of the band structure persist for smaller values of ε? Our findings are summarized in Figure 11 which demonstrates the bifurcation diagrams for two values: <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M127">View MathML</a> and <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128">View MathML</a>. They confirm that the band structure does persist and show predictably regions with larger number of spikes per bursts, especially for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128">View MathML</a> (compare the SC plot in 1D and 2D diagrams). To study the genesis of the band structure we also built the corresponding <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M130">View MathML</a>-diagrams shown in Figure 12. Although both diagrams represent the same data of the spike-counting, the diagram on the right in Figure 12 is given in a logarithmic scale to demonstrate that the self-similarity 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 <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M18">View MathML</a> can be equality singled out to perform the bifurcation analysis of the Hindmarsh-Rose model.

thumbnailFig. 11. (A) 2D spike-counting diagram projected to the plane <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a> at <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M127">View MathML</a>; (A1) interspike-duration bifurcation diagram for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a>. (B) 2D spike-counting plot in the plane <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M44">View MathML</a> for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M128">View MathML</a>; (B1) interspike-duration bifurcation diagram for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a> and (B2) 1D spike-counting diagram for <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M96">View MathML</a>.

thumbnailFig. 12. <a onClick="popup('http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.mathematical-neuroscience.com/content/1/1/6/mathml/M130">View MathML</a> spike-counting 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 Hindmarsh-Rose 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 Hodgkin-Huxley 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 slow-fast 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 mathematics-based 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 spike-counting approach, the evaluations of inter-spike 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 Hindmarsh-Rose 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 spike-addition/deletion cascades, phenomena of bistability, and various transitions between types of bursting: square-wave and plateau-like. This ensures that this toolkit of combined neuroscience-native methods and calculi-based 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 problem-specific 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 10-20 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 3-cell 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 side-by-side 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 MTM2009-10767, and by NSF Grant DMS-1009591, RFFI Grant No. 08-01-00083, GSU Brain & Behaviors program, and MESRF ‘Attracting leading scientists to Russian universities’ project 14.740.11.0919.

References

  1. Kristan, WB: Neuronal decision-making circuits. Curr. Biol.. 18(19), 928–932 (2008)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  19. Innocenti, G, Genesio, R: On the dynamics of chaotic spiking-bursting transition in the Hindmarsh-Rose neuron. Chaos. 19(2), (2009)

  20. Wang, XJ: Genesis of bursting oscillations in the Hindmarsh-Rose model and homoclinicity to a chaotic saddle. Physica D. 62(1-4), 263–274 (1993)

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

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

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

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

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

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

  27. Storace, M, Linaro, D, de Lange, E: The Hindmarsh-Rose neuron model: bifurcation analysis and piecewise-linear approximations. Chaos. 18(3), (2008)

  28. de Lange, E, Hasler, M: Predicting single spikes and spike patterns with the Hindmarsh-Rose model. Biol. Cybern.. 99(4-5), 349–360 (2008)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  44. González-Miranda, JM: Observation of a continuous interior crisis in the Hindmarsh-Rose neuron model. Chaos. 13(3), 845–852 (2003)

  45. González-Miranda, JM: Complex bifurcation structures in the Hindmarsh-Rose neuron model. Int. J. Bifurc. Chaos Appl. Sci. Eng.. 17(9), 3071–3083 (2007)

  46. Briggman, KL, Kristan, WB: Multifunctional pattern-generating circuits. Annu. Rev. Neurosci.. 31, 271–294 (2008)

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