Macro-and micro-chaotic structures in the Hindmarsh-Rose model of bursting neurons

We study a plethora of chaotic phenomena in the Hindmarsh-Rose neuron model with the use of several computational techniques including the bifurcation parameter continuation, spike-quantiﬁcation, and evaluation of Lyapunov exponents in bi-parameter diagrams. Such an aggregated approach allows for detecting regions of simple and chaotic dynamics, and demarcating borderlines—exact bifurcation curves. We demonstrate how the organizing centers—points corresponding to codimension-two homoclinic bifurcations—along with fold and period-doubling bifurcation curves structure the biparametric plane, thus forming macro-chaotic regions of onion bulb shapes and revealing spike-adding cascades that generate micro-chaotic structures due to the hysteresis. V C 2014 AIP Publishing LLC . [http://dx.doi.org/10.1063/1.4882171]

Understanding common dynamical principles underlying an abundance of widespread brain behaviors is a pivotal challenge in the new century. The bottom-up approach to the challenge should be based on solid foundations relying on detailed and systematic understanding of dynamical functions of its basic components-neurons-modeled as plausibly within the Hodgkin-Huxley framework as phenomenologically using mathematical abstractions. Such one is the Hindmarsh-Rose (HR) model, reproducing fairly the basic oscillatory activities routinely observed in isolated biological cells and in neural networks. This explains a wide popularity of the HR-model in modern research in computational neuroscience. A challenge for the mathematics community is to provide detailed explanations of fine aspects of the dynamics, which the model is capable of, including its responses to perturbations due to network interactions. This is the main focus of the bifurcation theory exploring quantitative variations and qualitative transformations of a system in its parameter space. We will show how generic homoclinic bifurcations of equilibria and periodic orbits can imply transformations and transitions between oscillatory activity types in this and other bursting models of neurons of the Hodgkin-Huxley type. The article is focused specifically on bifurcation scenarios in neuronal models giving rise to irregular or chaotic spiking and bursting. The article demonstrates how the combined use of several state-of-the-art numerical techniques helps us confine "onion"-like regions in the parameter space, with macro-chaotic complexes as well as micro-chaotic structures occurring near spike-adding bifurcations.

I. INTRODUCTION
The wide-range assessment of brain and behaviors is one of the pivotal challenges of this century. To move further in understanding how an incredibly sophisticated conglomerate such as the brain per se functions dynamically, at least it is imperative to know the dynamics of its elements-neurons. Having the exhaustive knowledge of behaviors of neurons should help us deepen our understanding the synergy of cooperative behaviors of populations of neurons-neuronal networks in this bottom-up approach. Dynamically, neurons are categorized and characterized by kinds and shapes of membrane potential oscillations. Waveforms of such oscillations have been proved to be useful for prediction and identification of various neurological deceases such as Parkinson or epilepsy. 1 A neuron is often viewed as a dynamical system, 2 which can be coupled into a multi-parametric family of nonlinear ordinary differential equations (ODEs) to describe neuronal activity in time. As such, the neuronal dynamics and its core elements such as robustness to perturbations and/or transformations depending on parameters become the subject of the bifurcation analysis. As control parameters of the system are varied, its dynamics undergo formal bifurcations in terms of the dynamical system theory that match qualitative metamorphoses in terms specific to neuroscience. The range of the non-stationary activity types is broad and includes regular and irregular tonic spiking, bursting and mixed-mode oscillations and combinations of them, as well as oscillatory transients toward quiescent states. In terms of the dynamical system theory, these would correspond to stable periodic and deterministically chaotic orbits in the phase space of the model. Individual neurons are modeled plausibly within the Hodgkin-Huxley framework, 2 as well as phenomenologically using mathematical abstractions. The 3D Hindmarsh-Rose model, 3 a) Author to whom correspondence should be addressed. Electronic mail: rbarrio@unizar.es reproducing fairly the basic oscillatory activities routinely observed in isolated biological cells and in neural networks, is such a popular abstraction among mathematicians.
In 1984, Hindmarsh and Rose proposed a 3D extension of the original tonic-spiking system to model and describe two-time scale oscillatory activity in bursting neurons. In this model, the x state variable should be treated as the voltage across the cell's membrane, while the y and z variables would describe kinetics of some ionic currents, respectively, fast and slow due to a small parameter 0 < e ( 1. Therefore, the HR model is a slow-fast system with the (x, y)-fast subsystem and single slow equation in z, in which e controls the time scale of z, while x 0 controls the rest potential of the system. Through the article, we will use the following values for these parameters: x 0 ¼ À1.6 and e ¼ 0.01 and e ¼ 0.001. The other parameters, b and I, will be the primary bifurcation parameters of the model in this study such that, by varying them, the model exhibits the whole spectrum of the activity types typical for endogenous bursters. A bifurcation study of the HR model sheds a light on qualitative properties of other neuronal systems of resemblant fast-slow dynamics. [4][5][6] Irregular or chaotic dynamics in neurons 7,8 and corresponding slow-fast models occur at transitions between neuronal activities such as spike adding, and between bursting and tonic spiking, which must be properly interpreted in terms of the theory of nonlocal bifurcations. 5,9,10 The goal of this article is to outline a general structure of chaotic regions in the HR and similar neuron models, which was first described by Linaro et al., 11 with a particular emphasis on a chain of narrow islands in the parameter plane, where chaotic attractors and stable periodic orbits can coexist. The occurrence of such islands related to spike-adding cascades is a typical feature of fast-slow systems with bursting dynamics. In what follows we will study and differentiate between two kinds of regions and islands, which are referred to as "macro" and "micro" chaotic structures, respectively, throughout the article. The macro structures due to spike adding cascades are sculpted by saddle-node (SN) (fold) and period doubling bifurcation curves originating in codimension-two homoclinic bifurcations. The chain of macro chaotic regions has a characteristic "onion"-like structure. An island of the micro chaos features the bi-stability and hence the hysteresis phenomenon.
This article is arranged as follows: Section II introduces briefly the numerical methods used in the article. Section III presents an organization of such large-scale regions and their generation. In Sec. IV, we investigate small-scale chaotic islands and coexistence regions. Summary is presented in Sec. V of the article.

II. NUMERICAL METHODS
An introductory classification scheme of bursting types in low-order neuronal models relies, in the very first approximation, on the slow-fast decomposition paradigm, where slow and fast dynamics in the subsystems do not overlay. Its core is a catalogue of codimension-one bifurcations that initiate or terminate slow motions, or critical manifolds composed of equilibria and foliated by limit cycles of the fast, (x, y)-subsystem in the "frozen"-singular limit e ¼ 0, at which the z-variable becomes the swiping parameter. All such bifurcations in a plane have been known for nearly a century. 12 The dynamics of the combined system singularly perturbed at small e 6 ¼ 0 concentrate around normally hyperbolic sections of the slow motion manifolds, [13][14][15] which are made of stable (hyperbolic in general) equilibria and periodic orbits staying away from bifurcations in the fast subsystem. Such sections constitute bare bones of most bursting patterns in slow-fast neuron models. A typical model of a (squarewave) burster possesses a pair of such manifolds, 16,17 respectively: quiescent, M eq , and tonic spiking, M lc (Fig. 1).
The slow-fast decomposition approach warrants a drastic reduction for dynamics of a singularly perturbed system. By combining the modulatory slow dynamics with the topology of the manifolds in the fast subspace, we are able to adequately interpret a plethora of dynamical phenomena observed in individual neurons and their plausible models. The approach, however, does not provide full explanations of (global) bifurcations underlying activity transitions in slow-fast neuronal systems, as the decomposition does not account for reciprocal interactions between the slow and fast dynamics that give rise often to complex dynamical phenomena only occurring in the whole system. 10,[18][19][20] This article is aimed to shed a light on an organization of such global bifurcations at transitions between regular and chaotic tonic spiking and bursting in the HR model, as well as to analyze the structure of its parametric space by combining several computational methods. Several recent works [4][5][6]11 were specifically focused on detailed studies of global bifurcations of tonic spiking and bursting orbits giving rise to chaotic dynamics in the HR model. Most computational studies of complex dynamics in neuron models typically consolidate the numerical continuation 21,22 for bifurcations of homoclinic and periodic orbits and neuroscience-native spike-quantification (SPQ) approaches to reveal universal structures in parameter spaces.
The SPQ is based on the location of the attracting periodic orbit after a transient time (in our case 10 3 time units) and counting the number of spikes of the resulting periodic orbit. Note that earlier two papers 4,6 introduced the concept of biparametric spike-quantification sweeps to identify and locate subregions of the occurrence of quiescent states, spiking, and bursting oscillations including chaos. To that goal, we use a robust and precise numerical method for ODEs: the free software TIDES, 23 that is based in the Taylor series method. 24 This method is quite adaptive and provides dense output, quite useful for event location used in the detection of periodic orbits, as it is the case because we obtain the maxima of the time-series and we search for very close approaches (distance lower than 10 À6 ). In the numerical integration, we use the variable-stepsize and variable-order integrator provided by TIDES, by means of an error tolerance of 10 À12 . The bi-parametric sweeping is done via a regular grid of 1000 Â 1000 points in the parameter space using as initial conditions for each point the final position of the previous point simulation.
The spike-quantification approach for a neuron model works well when a bursting solution follows closely the slow motion manifolds M lc and M eq of the fast subsystem and makes pronounced rapid jumps between them, thus defining the number of spikes per bursts in the voltage traces. Figure  1 illustrates square-wave bursting activity with robust five spikes for x 0 ¼ À1.3, b ¼ 3, I ¼ 5, and e ¼ 0.002 in the HR model. Indeed, the spike number within a burst is that of the complete revolutions of the bursting orbit around the spiking manifold M lc . One can observe from this figure that the interspike interval increases by the end of each burst: this is a signature of the square-wave bursting, which indicates the homoclinic bifurcation terminating the manifold M lc , when it touches a saddle point on the middle section of M eq . The location of this point can be used as a threshold separating the hyperpolarized and depolarized, stationary or oscillatory states of a neuron. In the spike-quantification technique, a fixed number of spikes per bursts is an indication of regular bursting, while unpredictably varying numbers are associated with chaotic bursting, which can be further supported by the computational approach utilizing the evaluations of the Lyapunov exponents.
The Lyapunov exponents are obtained using the classical algorithm proposed by Wolf et al. 25 The numerical integration is done again using the software TIDES with a tolerance 10 À12 and output each time unit for using the Wolf's algorithm. A transient time of 10 3 time units is used and the final simulation time is 10 4 .
Those two sweeping techniques are combined with the parameter continuation approach, which is bifurcation theory native, to uncover a fine structure of the bifurcation set of the HR model. The bifurcation parameter continuation done in this article is based on classical continuation theory using the well known free available software AUTO 21 for bifurcations of equilibrium points and periodic orbits. This software permits to locate the curves of codimension-one bifurcations and the main codimension-two bifurcation points in such a curve. Later, a detailed study of the influence of these curves and points gives key information of the system.

III. MACRO-CHAOTIC STRUCTURES
This article addresses the question concerning the global organization of the bifurcation structures of the parameter space of the HR model. Earlier two papers 4,6 introduced the concept of biparametric sweeps to identify and locate subregions of the occurrence of quiescent states, spiking, and bursting oscillations including chaos. Figure 2 represents two biparameter sweeps of the model in the (b, I)-plane at e ¼ 0.01 and 0.001 that are done with the spike-quantification approach. The parameter plane is clearly demarcated into regions corresponding to periodic tonic spiking, chaotic and regular bursting of the square-wave and plateau-like types. The obtained maps are color-coded so the spike numbers are associated with specific colors. The resulting diagram can be easily read and interpreted: the region shown in a dark blue is for stable spiking activity which can be treated as bursting with a single spike. The diagram reveals a global organization FIG. 2. Two (b, I)-parametric sweeps of the HR model, which is based on the spike-quantification approach, at (a) e ¼ 0.01 and (b) e ¼ 0.001. The parameter plane is clearly demarcated into regions corresponding to periodic tonic spiking, chaotic and regular bursting of the square-wave and plateaulike types. Stripes of gradually changing colors correspond to bursting with incrementally varying numbers of spikes due to a spike adding cascade. Bursting becomes chaotic near the transitions to tonic spiking in a chain of "onion"-like regions.
of spike adding bifurcations occurring on borderlines between the corresponding stripes in the blue hue, which all correspond to square-wave bursting on the model. Stripes of gradually changing colors correspond to bursting with incrementally varying numbers of spikes due to a spike adding cascade. Bursting becomes chaotic near the transitions to tonic spiking in a chain of "onion"-like regions. A sudden change in the number (colored in an orange hue) of spikes per bursts is associated with a transition from the square-wave to plateau-like bursting that occurs in the top-left corner of the diagram. While the terminal of the active phase of square wave bursting is due to a homoclinic bifurcation in the fast subsystem of the HR-model, that of plateau-like bursting is caused by saddlenode bifurcation of the depolarized equilibrium state in it after the spiking manifold M lc becomes tangent and stops touching the middle, saddle branch of the quiescent manifold M eq . 5 Decreasing the value of e ¼ 0.001, which determines the dynamics of the slow z-variable, does not change qualitatively the structure of the parameter plane but compresses it. Therefore, for simplicity without loss of generality we consider e ¼ 0.01.
The large scale diagram in Fig. 2 for e ¼ 0.01 helps us identify an organization of a multilayer area secluded around b 2 ½2:5; 3:3 and I 2 ½2:2; 4:5, which should carry over the typical structures for all such regions in the "onion"-like chain. Recently, the first small chaotic region around b ¼ 3.2 and I ¼ 2.5 has attracted some attention in some experimental studies. 26 In addition, we will get a closer look in regions bordering transitions from square-wave to plateau bursting, which is white boxed. The magnifications of the region are presented in Fig. 3.
Panel (a) in Fig. 3 represents the bi-parametric sweep using the Lyapunov exponents, k 2 k 1 . The colormap of the diagram is designed as follows: the colors, from blue to red are associated with the increasing value of the first, maximum Lyapunov exponent, k 1 > 0, thus indicating chaotic dynamics and quantifying the disorder degree. Whenever k 1 ¼ 0 on a periodic orbit within its existence region, we evaluate the second Lyapunov exponent, k 2 that determines its stability. Negative values of k 2 are associated with colors of grey shades, so that a black means that k 2 is close to zero from the left, which means that the corresponding multiplier of the periodic orbit is close to þ 1 or À1 and it is about to undergo a period doubling or a saddle-node bifurcation. 27 So, a single black line passing through symmetrically grey area of the diagram should be interpreted as the stable periodic orbit loosing stability through period doubling bifurcation. In the case of a saddle-node bifurcation, such a black line is surrounded by adjoin regions of some asymmetric hue on sides, which would be an indication of bistability in the system, through which the phase point switches to another attractor. It is seen clearly from panel (a) in Fig. 3 that the transition from tonic spiking to square-wave bursting must always passes through a strip of chaotic dynamics. The chaotic region produces a chain of shrinking extensions following the curves of spike adding cascade. It becomes evident that such extensions carry over the same properties, it suffices to figure out the bifurcation overlay in one such a "bulb." Fig. 3 represents bi-parametric sweep using the spike-quantification. The linear colormap of the diagram must be interpreted as follows: blue stands for tonic-spiking, or single-spike bursting. As the spike number increases, the color shifts towards the hot side of the spectrum. Circled numbers indicate the quantities of spikes per burst in the chosen locations in the parameter plane. We can clearly see the line demarcating transitions from square-wave to plateau bursting after the spike number jumps at once to 25 and greater, which are represented by orange-red in the diagram. The number of spikes fluctuating substantially over the time and reaching a threshold value for some single parameter values is associated with a dark red, thus indicating the onset of chaotic bursting in the model. There is a perfect agreement between both diagrams, Figs. 3(a) and 3(b), both revealing in detail the complex organization of the parameter space along with chaotic regions at the tonic-spiking and bursting transitions, as well as bursting metamorphoses including spike adding cascades. 6 Panel (c) in Fig. 3 represents a spike-quantification sweep of the HR model at e ¼ 0.001. Decreasing e leads automatically to proportional increases of spikes per regular bursts so that while the biparametric diagram becomes overpopulated and it persists qualitatively similar structures though considerably condensed.
In conclusion of this section, we recapitulate on the approaches: sweeps based on the spike-quantification organically complement the suit of tools for computational examinations of the dynamics and bifurcations in slow-fast models of neurons. This approach, natural for neuroscience, allows for both pilot and detailed studies of individual and networked models of oscillatory cells, detections of principal bifurcations underlying various intrinsic transformations (spike-adding) and transitions between activity types such as quiescence, regular and chaotic tonic-spiking and bursting of different kinds. While qualitatively the spikequantification approach works as effectively as the tools based on Lyapunov exponents, for swiping bifurcation diagrams, the former one is significantly faster and besides provides quantitative/statistical measurements available as well.

A. Bifurcation skeleton
In what follows we will augment the sweeping techniques with the parameter continuation approach for the HR model at e ¼ 0.01. 11 The results are shown in the two panels of Fig. 4: the full size diagram (a) and the magnification of its particular subregion (b). In these figures, the red and yellow lines correspond to period-doubling (PD) and SN bifurcations, respectively. The homoclinic bifurcation curve is drawn in the black; it bends into a characteristic U-shape. 11 As was said above, the consideration will be focused around the two largest (chaotic) regions in the chain.
The feature of square-wave bursters is the homoclinic bifurcation of a saddle. The homoclinic loop of the saddle ends up the stable slow-motion manifold, M lc , in the extended space of the fast subsystem of the HR model, while the saddle sets a voltage threshold between tonic spiking and hyperpolarized quiescent in the waveforms (Fig. 1). The dynamical consequences of this plain bifurcation in the full model are no longer that simple. The reason is that as soon as the static z becomes a dynamic variable, the saddle changes the leading unstable direction, which is due to the slow equation (1). This implies that a separatrix loop at e 6 ¼ 0 makes an orbit flip and approaches the saddle not in the (x, y) subspace but along the middle branch of M eq , see more details in Ref. 5. Such a transformation is called a homoclinic orbit flip bifurcation, 9 which is in general of codimension-two, and hence corresponds to an isolated point on the homoclinic bifurcation curve in the parameter plane. A few such identified points are represented by the green dots, OF, in Fig. 4. Depending on the magnitude of the characteristic exponents of the saddle, the unfolding of the orbit-flip bifurcation may vary. 9 In either case (Fig. 5), the loci of the unfolding includes curves corresponding to saddle-node and period-doubling bifurcations, as well a duple-pulse (secondary) homoclinic bifurcation (bifurcation type B 28 ), which is referred to as a homoclinic doubling bifurcation nowadays. In other cases, like one under consideration (type C), it can be more populated by a countable set of such bifurcation curves for longer, n-pulses homoclinics and matching periodic orbits. 9,28,29 Such an orbit-flip bifurcation gives rise to the onset of complex dynamics of a finite-shift type due to a formation of topological Smale horse-shoes in the return maps. Such dynamics occur near tonic-spiking and bursting transitions in square-wave bursters. 30 We will outline below how such codimension-2 points can organize globally the bifurcation set of the HR model.
In the model under consideration, countable sets of subsidiary bifurcation curves are originated in a vicinity of each orbit-flip bifurcation point of the type C. 11 Moreover, all the period-doubling and saddle-node bifurcation curves end up at the same (or several nearby) inclination flip (IF) bifurcation point at the top-left area. Due to abundance and density of such points, it is computationally hard to determine very accurately the terminals for some bifurcation curves. Together, they form a "palm" or "onion"-like structure (Fig. 6). It is seen from Fig. 4 that the formation of a chaotic region (in dark red) begins with a period doubling cascade (red curves), and it terminates with the saddle-node bifurcation (yellow curve) that generates a stable periodic orbit existing in the stability region (blue) between two successive chaotic layers, due to a spike increment.
The very tip of the chaotic region gives rise to the beginning of a line corresponding to the spike-adding in bursting. We call the chain of such chaotic organizations the macrochaotic structure of the HR-model.
Let us examine construction details of these macro-chaotic structures with the use of a single parameter pathway transversally passing throughout the chain of the chaotic regions. Such a pathway can be a line segment given by I ¼ (1 À 0.265 b)/0.0691, where b 2 ½2:6; 3:2, which is drawn in a brown in  Figures 8 and 9.
FIG. 5. Two unfoldings of homoclinic orbit-flip bifurcation of codimensiontwo. 9,28 Unfolding of type B contains three primary curves corresponding to SN and PD bifurcations of periodic orbit and to a homoclinic bifurcation of a double (secondary) separatrix loop. Unfolding of type C has a complex locus of curves for bifurcations of various n-pulse homoclinic orbits and PD and SN bifurcations of matching periodic orbits. The type of the unfolding depends on the magnitudes of the characteristic exponents of the saddle equilibrium state. exponent on attractors plotted against b. Both diagrams show alternation between periodic windows and chaotic regions. Chaotic dynamics occurring within each region is universally organized: it begins with a period-doubling cascade of a (n)-spike bursting orbit on one side and ends up by a saddle-node bifurcation underlying intermittency of type I transitioning to a stable (n þ 1)-spike bursting orbit on the other side of the region. As the value of b decreases, bursting orbits with larger spike numbers will be observed in the stability windows and so forth. The detailed examination confirms that the orbits, being created on one side, continue till they vanish in reverse period-halving cascades. This supports the results in Ref. 7 that the organization of macrochaotic structures in the Hindmarsh-Rose model is due to a sequence of forward period-doubling and reverse periodhalving cascades. 31 Considered next are chaotic dynamics and their metamorphoses on the single-parameter pathway I ¼ (1 À 0.265 b)/0.0691 through the chain of chaotic regions in the (b, I)plane. Fig. 8 represents a showcase of chaotic attractors for designated points (black dots in Fig. 7) with positive values of the maximum Lyapunov exponent. Plotting the attractors next to each other lets one assess the evolution of their shapes in qualities of square-wave bursters. The evolution starts at b ¼ 2.626 with chaotic tonic spiking. With increasing b to 2.635, bursting incorporates a large number of spikes thus increasing its period. Increasing b further makes the bursting attractor loose spikes. One way to quantify the shape of a chaotic attractor is to find its bare-bone made of minimal unstable periodic orbits (UPOs) (all periodic orbits of low multiplicity); this is a core of UPOs approach.
The numerical algorithm for the search of such periodic orbits is based on the stability transformation (ST) method combined with the damped Newton method. 32 The ST method 33 was specifically developed to locate unstable periodic orbits in chaotic attractors of discrete dynamical systems. The approach was generalized for time-continuous systems through reductions of the original continuous flow to a Poincar e return map, and then to detection of the desired fixed points. 34 An advantage of the method is that it let one investigate less constrained and hence wider basins for locations of fixed points in question compared to the Newton method. 33 The number of all UPOs for values of b marked in Fig. 7 is summarized in Table I. We present, for each selected value b (column 1), the number of periodic orbits with m spikes, up to m ¼ 5, inside the corresponding chaotic attractor (columns 2 to 6). The number of UPOs embedded on the chaotic attractor and their distribution determines its structure. 35,36 Taking into account the "onion"-like structure of the set of chaotic layers, chaotic attractors with the same topological structure exist with big and small values of b. For example, we can see in this table that attractors in the extremes (b ¼ 2.626 and b ¼ 3.05) present the same number of UPOs for each quantity of spikes (up m ¼ 5). This fact is not enough to say that their topological structure is the same, but it guarantees some similarity in the attractors. 37 The chaotic attractors sampled at b ¼ 3.05 and 2.69 are superimposed by their skeletons (formed by the minimal unstable periodic orbits) in Fig. 9. We can see how the distribution of the orbits along the attractor is very different in both situations. While for b ¼ 2.69, the periodic orbits cover regions of chaotic attractor with smaller values of z when the number of spikes, m, increases. As for b ¼ 3.05, from m ¼ 2, unstable periodic orbits go over all range of values of z where the chaotic attractor exists. It should be noted that in the first case (b ¼ 2.69), all orbits pass close to each other. This makes it more difficult to find them. Actually blue and black cases in m ¼ 4 have been found using only ST method since the Newton method does not converge even for points very near to the solution.  Fig. 7).

IV. MICRO-CHAOTIC STRUCTURES
In this section, we examine less pronounced regions of the co-existence of chaotic bursting with stable periodic orbits. Such regions featuring bistability are proposed to be termed micro-chaotic structures.
By revisiting Fig. 4, we note that each chaotic layer starts from the vicinity of an orbit flip bifurcation point. If we pay attention to panel (b) of the same figure, we can observe that bifurcation curves related to each layer come from different OF points in the homoclinic curve. The first one (PD 1 ) comes from an OF point with smaller b values than the OF point generating the countable pencil of PDs.
This situation generates that the first period-doubling bifurcation curve crosses the rest of bifurcation curves of the chaotic layer. Panel (a) of Fig. 10 displays a magnification of the region where these crossings occur for the first chaotic layer. In panel (b), we can see the bifurcation diagram along the green segment of (a). This diagram shows that two attractors coexist for values of parameter b higher than the corresponding value for the first period-doubling bifurcation (an example is plotted in panel (d)). In fact, slightly lower values also show coexistence of attractors. This is because the fold bifurcation (yellow dotted) curve is almost touching the first period-doubling bifurcation and makes one of the attractors becomes unstable. These two bifurcation curves correspond to attractor with fewer spikes (in red). During coexistence, this attractor has a basin of attraction larger than the other, being the dominant. After the fold, it becomes unstable and the other attractor dominates the evolution of the system. The attractor with the largest number of spikes (in black) experiences the rest of bifurcations and becomes chaotic (see panel (c)). Therefore, the coexistence region contains an island where a chaotic attractor with a small basin of attraction exists. This is which we have called a micro-chaotic structure. Figure 11 documents the case of bistability that is due to a hysteresis loop. Such hysteresis of coexisting attractors, occurring in a variety of nonlinear systems is organized by two saddle-node bifurcations interconnected in a closed loop. It is seen from the schematic explanation in Fig. 12 that in between two turning points corresponding to saddlenodes, the system exhibits bistability of two attractors whose basins are separated by an unstable threshold. Such attractors can be as simple equilibria and periodic orbits as a result of a sequence of bifurcations, like period doubling or spike adding, producing, respectively, chaotic tonic spiking or bursting.
The HR model can exhibit several types of bistability that occur in a narrow parameter regions close to the OF points. A typical bistability is the coexistence of two stable periodic orbits, one of which with an extra spike, as shown in Figure 11. 6,11 It is predictably observed within a relatively broad range near every spike-adding bifurcations.
Spike-adding cascades have happened to be endogenous bifurcation phenomena for various fast-slow systems. For the HR at I ¼ 2.15, panels (b) and (e) of Fig. 11 present the b-parameter continuation diagram of the PSS of two bursting orbits with two (in blue) and three (in red) spikes. One can see from the diagram that the bistability of these orbits is due to the hysteresis.
Another less common situation is the hysteresis underlying bistability of a periodic orbit and weakly chaotic attractor. It occurs from a quick period-doubling cascade and hence looks like a narrow band in panels (a) and (c) of Fig. 11. The corresponding pathway is taken on the line segment, given by I ¼ À25/9(b À 2.92) þ 2.7, that cuts through the spike-adding region and parallel to the primary homoclinic bifurcation curve (in black). Closeness of the chaotic attractor to the threshold saddle orbit (dashed line) explains that out of two, the stable orbit will dominant in this bistability competitions. As before, we refer to the existence region of the other attractor as a micro-chaotic structure. Such micro structures are located nearby the OF homoclinic bifurcations on the parameter plane. We remark that in this region of spike-adding cascades another micro-chaotic behavior (MMBOs) have been discovered recently 38 related with canard phenomena.
The summarizing Fig. 13 discloses schematically the likely global wiring by the curves corresponding to saddlenode, period doubling, and homoclinic bifurcations that partition the (b, I)-parameter plane of the HR-model. 11 It singles out the relative locations of spike addition transitions as well as the islands of microscopic and macroscopic chaos, which are caused by and constrained within the unfolding of the orbit-flip codimension-two bifurcations, which are essential for understanding the complex dynamics at a transition between tonic spiking and bursting in square-wave bursters. 5,11 FIG. 12. Schematic explanation of hysteresis loops with two lines involving a pair of periodic orbits co-existing (a) or one stable periodic orbit and a micro-scaled chaotic attractor caused by period doubling cascades (b).
FIG. 11. Sketch of infolding structures near two orbit-flip (OF) homoclinics including SN, PD, and spike adding bifurcations. Two consecutive SN bifurcations set frames for the hysteresis loops of coexisting bursting orbits with 2 and 3 spikes (panels (b) and (e)) on the yellow dashed-segment on I ¼ 2.15, or 3-spike bursting and a chaotic attractor (panels (a) and (c)) due to PD-bifurcations on the pathway I ¼ À25/9(b À 2.92) þ 2.7, close to the homoclinic bifurcation curve (in black).

V. SUMMARY
Two-time scale bursting activity has been long modeled by slow-fast odes. To understand bursting, one has traditionally used the decomposition approach dismantling the system to basic components-critical or slow motion manifolds. The approach has proved to work well for spiking and bursting models unless one is interested in particular bifurcations and transition between activity types. For square-wave bursters, the key transformations that shape the dynamics are due to homoclinic bifurcations. The first is a simple homoclinic bifurcation that ends up the cylinder-shape manifold spanned by periodic orbits corresponding to tonic spiking in neuron models. This homoclinic bifurcation becomes much less trivial with an addition of slow-dynamics necessary to describe square-wave bursting. The slow-fast scales make the bifurcation degenerate and drastically complicate its unfoldings. A feature of such bifurcations is that they can cause chaotic dynamics with overlapping complex bifurcation structures including saddle-node and period doubling bifurcations of periodic orbits. As such, a multilayer passage, through such structures in the parameter space of the neuron model, reveals a plethora of dynamical phenomena including the aforementioned bifurcations of tonic spiking and bursting orbits, a hierarchy, based on self-similarity of small and large scale chaotic dynamics, various kinds of bistability of the periodic and chaotic attractors, spike adding and more. A comprehensive understanding of such metamorphoses of homoclinic bifurcations and impending consequences on overall dynamics in singular perturbed system is yet to come.
This article is aimed to disclose the main components assembling the bifurcation structure of the Hindmarsh-Rose model exhibiting all above dynamic features. Given the breadth and ranges of non-local bifurcations underlying transformations of oscillatory activity including chaos onsets, we rely heavily on computational explorations of the dynamics of the model. We have proposed and combined several numerical techniques to uncover and explain the organization of its bifurcation structure. They include but not limited to biparametric spike-quantification diagrams, Lyapunov exponent sweeps, parameter continuation to locate individual bifurcations of periodic and homoclinic orbits. With this consolidated approach, we were able to identify the modulatory bifurcation structure in the chain form of chaotic islands centered around the codimension-two orbitflip homoclinic bifurcations.
At the aftermath of our detailed numerical studies, we come up with the symbolic representation of the bifurcation organization that we believe would be typical for models of square-wave bursters other than the Hindmarsh-Rose model under considerations. As in most cases of systems with complex chaotic dynamics, it is impossible to embrace all, but principle elements of constructions of infinite complexity.