Homoclinic organization in the Hindmarsh–Rose model: A three parameter study

Bursting phenomena are found in a wide variety of fast–slow systems. In this article, we consider the Hindmarsh–Rose neuron model, where, as it is known in the literature, there are homoclinic bifurcations involved in the bursting dynamics. However, the global homoclinic structure is far from being fully understood. Working in a three-parameter space, the results of our numerical analysis show a complex atlas of bifurcations, which extends from the singular limit to regions where a fast–slow perspective no longer applies. Based on this information, we propose a global theoretical description. Surfaces of codimension-one homoclinic bifurcations are exponentially close to each other in the fast–slow regime. Remarkably, explained by the specific properties of these surfaces, we show how the Hindmarsh–Rose model exhibits isolas of homoclinic bifurcations when appropriate two-dimensional slices are considered in the three-parameter space. On the other hand, these homoclinic bifurcation surfaces contain curves corresponding to parameter values where additional degeneracies are exhibited. These codimension-two bifurcation curves organize the bifurcations associated with the spike-adding process and they behave like the “spines-of-a-book,” gathering “pages” of bifurcations of periodic orbits. Depending on how the parameter space is explored, homoclinic phenomena may be absent or far away, but their organizing role in the bursting dynamics is beyond doubt, since the involved bifurcations are generated in them. This is shown in the global analysis and in the proposed theoretical scheme.


ARTICLE
scitation.org/journal/cha is one of the fields where these models are more abundant. In neuroscience, to understand how an incredibly sophisticated system such as the brain per se functions dynamically, it is imperative to study the dynamics of its constitutive elements-neurons. Since Hodgkin and Huxley developed the first model of action potentials in the membrane, 1 the design of mathematical models for neurons has arisen as a trending topic in science for a few decades, and a lot of models and variations describing different kinds of neuron cells in numerous animals have been proposed in the literature. What all these systems have in common is the existence of fast-slow dynamics, 2 that is also quite usual in a lot of other practical applications, like in chemical reactions 3 and laser dynamics. 4 In all these models, one of the key magnitudes is the time that a neuron, or other dynamical system, is active, and this is related to the number of oscillations (spikes) in the fast subregime.
In order to help in the analysis of neuron models simulated realistically within the Hodgkin-Huxley framework, 1 a common approach is to use some simplified models. In particular, the 3D Hindmarsh-Rose (HR) model 5 reproduces fairly well the basic oscillatory activities routinely observed in isolated biological cells and in neural networks. It fulfills the two basic conditions of being computationally simple but, at the same time, able to reproduce the main behavior (the rich firing patterns) exhibited by the real biological neuron. The HR model is described by three nonlinear ordinary differential equations, where x is the membrane potential, y the fast, and z the slow gating variables for ionic current. In our study, we will consider a typical choice of parameters: a = 1, c = 1, d = 5, s = 4, and x 0 = −1.6. Parameters b and I determine the bursting or spiking behavior and their values are considered in specific ranges where such phenomena are present. Parameter ε governs the fast-slow behavior and we will study dynamics for ε small, but including scenarios far from the singular limit ε = 0. In the sequel, we consider (1) as a family of vector fields depending on parameters (b, I, ε), and, say, fast subsystem to refer to the z-family obtained after taking ε = 0.
Roughly speaking, we can say that a fast-slow system exhibits bursting when orbits exhibit periods of fast spiking followed by periods of quiescence. When the jump between these two different regimes can be explained by a fold bifurcation of equilibria and a homoclinic bifurcation of periodic orbits (both bifurcations occurring in the fast subsystem), we say that the bursting is of fold/hom type. 6 In Sec. II (see Fig. 3), we will describe how fold/hom bursters arise in the HR model.
One of the big challenges regarding bursting phenomena is to understand the mechanisms explaining the variation in the number of spikes ( Fig. 4 in Sec. II B provides an illustrative example in the HR model). These spike-adding processes have been studied for several mathematical neuron models (see, for example, Refs. [7][8][9], and also in other contexts such as laser dynamics, chemical reactions, or discrete maps, with the alternative name of periodadding. [10][11][12][13][14] This process is quite important in that it progressively modifies the spectrum of periodic orbits of the system and the structure of chaotic attractors. [15][16][17][18] As argued by Terman, 18 these transitions may be either continuous, with the period of the bursting solution increasing along the process or they may involve chaotic behaviors (see also Ref. 19). Recently, these transitions have been studied in detail 20 providing a theoretical scheme for ε fixed. The relevance of fold bifurcations of periodic orbits in this process was pointed out earlier in Ref. 21. Dealing with fold/hom bursting, the spike-adding process has also been related to the existence of canard orbits [22][23][24][25] and with the existence of certain codimensiontwo homoclinic bifurcations. 15,26,27 Working with a fixed value ε = 0.01, the role of homoclinic bifurcations of codimension-one and -two in the spike-adding mechanisms was discussed in Ref. 26 and some preliminary results were advanced. In particular, bifurcations of periodic orbits around flip and Belyakov bifurcations (see Sec. II A for background) were identified as crucial ingredients to understand some spike-adding transitions that are present in the HR model. Again, working with that fixed value of ε, codimension-two homoclinic bifurcations were again considered in Ref. 27, but providing a much more thorough study. Different homoclinic curves were discussed and their sharp fold points were already detected in that reference and linked to the spike-adding processes. Codimension-two homoclinic bifurcations in Refs. 26 and 27 are also organizing centers of chaotic regions in the bifurcation diagram. All these chaotic phenomena were discussed in Ref. 15.
What is missing in the literature is a global study of how homoclinic bifurcations are organized, and to that goal we need, at least, to describe them in a three-parameter space. Note that it is intrinsic to the notion of bifurcation the possibility of observing its effects without the bifurcation point being present. In the HR model, one can explore the parameter space without detecting homoclinic bifurcations (see Fig. 4), although their consequences (fold and period-doubling bifurcations) are exhibited. The organizing points (the codimension-two homoclinic bifurcations) may be placed far away in the space or parameters, and even, they may be outside a particular set of parameters that we are visualizing, but they continue being the organizing centers. Taking all of this into account, the goal of this article is to provide a model of the homoclinic organization that explains all these facts.
As already mentioned, previous work in the literature was focused in studying, for some ε fixed, the curves of homoclinic bifurcation at equilibria displayed by the system. 15,26,27 A bifurcation diagram in a three-parameter space, including variation of ε, was first considered in Ref. 46. Changes in the spike-adding structures and the underlying bifurcations were observed. Moreover, foldings in the curves of inclination flip (IF) bifurcation were already detected. In Ref. 20, a theoretical scheme giving a complete scenario of bifurcations involved in the spike-adding processes in fold/hom bursters was introduced. This theoretical scheme provides a complete description of the connections of the different codimension-two points and the organization of the homoclinic curves for ε fixed. Also, in that paper, the validity of the scheme is checked for a pancreatic β-cell neuron model.
In this article, we are interested in understanding the global structure of the homoclinic surfaces in the three-parameter space. To that goal, a detailed numerical study with continuation techniques is required (we use the well-known software AUTO 28,29 )  Supported by numerical evidences, we conjecture that the intersection of each homoclinic surface with horizontal planes (with ε fixed) produces isolas in the plane of parameters (compare with results in Ref. 31 for the FitzHugh-Nagumo system), that is, simple closed curves in the corresponding slice. We show how, for each ε fixed, the model exhibits a finite number (number that grows when the small parameter decreases) of isolas corresponding to primary homoclinic bifurcations. Isolas are not only exponentially close to each other but they also exhibit a pair of extremely sharp folds so that the width of each isola is also exponentially small. These folds allow two sides of the isola to be distinguished (and also two faces of the surface of homoclinic bifurcations). On one of the faces, the corresponding homoclinic orbits on the fold/hom regime exhibit n spikes and, on the other face, n + 1. It is because of this fact that, from now on, we use the notation hom (n,n+1) to refer to the different homoclinic bifurcation surfaces (or isolas if working with two-parameter plots  27, the authors use the notation hom (n) , where we use hom (n,n+1) . Nevertheless, one should note that when required (see Figs. 4,5,and 7 in Ref. 27), they also use two different notations for a unique curve of homoclinic bifurcation, changing the label from hom (n) to hom (n+1 a) after a sharp fold of the curve is crossed, pointing out that the number of spikes changes from n to n + 1.
Homoclinic surfaces are the main focus of this article. We show how they are disposed in the parameter space, taking into account that, as numerics show, they are exponentially close to each other when ε → 0. Because of their tubular shape and the proximity of the surfaces, we can compare the whole structure with a "mille-feuille" pastry. There, we observe pencils of curves of fold and period-doubling (PD) bifurcations of periodic orbits generated on codimension-two bifurcation points. Moving ε, each of the curves in the pencil gives rise to a surface. Hence, we can compare the codimension-two bifurcation curves with the "spines-of-abook" with pages correspondent to surfaces of bifurcations of periodic orbits. Besides, the ε-level reached by each surface hom (n,n+1) decreases as n increases. This allows us to explain the simplification mechanisms (bursting with a lower number of spikes) that can be observed as ε increases.
The article is organized as follows. In Sec. II, we provide general information about the HR model: fast-slow decomposition, spikeadding process linked to fold/hom bursters exhibited by the model, and a discussion about existence of equilibria in the full system. A short survey about homoclinic bifurcations is also provided in Sec. II. In Sec. III, we pay attention to some particular slices (with ε fixed) inside the three-parameter space. Here, we show how the base shape of the homoclinic curves evolves as ε varies, but much more significantly, how the codimension-two homoclinic bifurcation points move on the homoclinic curves and, in fact, how they disappear when ε grows. Section IV presents a three-parameter study explaining some of the phenomena that are observed when ε is fixed and shows isolas of codimension-one homoclinic bifurcations. Section V introduces the global theoretical scheme creating the structure that we call "homoclinic mille-feuille," bearing in mind the codimension-one bifurcation surfaces. In them, we find "spines-ofa-book," bearing in mind the codimension-two bifurcation curves, holding pencils of bifurcations of periodic orbits. Both structures give rise to the theoretical model proposed in this article. Finally, we present some conclusions in Sec. VI.

II. BACKGROUND
In this section, we recall some basic aspects about homoclinic bifurcations and fast-slow dynamics, including a description of foldhom bursters, one of the mechanisms exhibited by the HR model for the creation of bursting orbits. In addition, to facilitate further discussions, the equilibrium points displayed by the full system (1) are explained. In our revision on bifurcations, only those that play a relevant role in the global organization of dynamics in the HR model are included. Consider a smooth family of vector fields X µ on R 3 with µ ∈ R k and suppose that there exist µ 0 ∈ R k and p 0 ∈ R 3 such that p 0 is a saddle-type hyperbolic equilibrium of X µ 0 . Without loss of generality, we assume that µ 0 = 0 and p 0 = 0. Let W s (0) [respectively, W u (0)] be the stable (respectively, unstable) invariant manifolds of X 0 at 0. Up to time reversal, we can assume that dim(W s (0)) = 1.

A. Homoclinic bifurcations
To state certain conditions, we will need to use the notions of strong unstable manifold and center stable manifold. Assume that DX 0 (0) has real eigenvalues λ s , λ u and λ uu with λ s < 0 < λ u < λ uu . The strong unstable manifold W uu (0) is a one-dimensional invariant manifold whose tangent space at 0 is given by the eigenspace corresponding to the eigenvalue λ uu (the so called strong unstable direction). On the other hand, the center stable manifold W cs (0) is a two-dimensional invariant manifold whose tangent space at 0 is given by the eigenspace corresponding to the eigenvalues λ u and λ s . Let 0 ⊂ W s (0) ∩ W u (0) be a homoclinic orbit. In the sequel, we assume that the family X µ unfolds 0 generically. To understand this condition, consider a cross section at a point in 0 and define the distance (µ) between the point W s (p µ ) ∩ and the curve W u (p µ ) ∩ , where p µ denotes the saddle type hyperbolic equilibrium of X µ , which exists close to 0 for µ to be small enough. We say that 0 is generically unfolded with respect to µ if D µ (0) ̸ = 0. Under this generic assumption, there always exists a hypersurface H in the parameter space such that 0 ∈ H and X µ has a homoclinic orbit asymptotic to p µ for all µ ∈ H.
There exist four classes of codimension-one homoclinic orbits.
Conditions λ s + λ u ̸ = 0 and λ s + ρ u ̸ = 0 are non-resonance hypothesis. Condition (H1) implies that 0 is tangent to the weak unstable direction, that is, the direction given by the eigenspace associated with the weak unstable eigenvalue λ u . Condition (H2) is a "non-inclination" property.
In Case 1 and Case 3, a single unstable (repelling) periodic orbit is born from the homoclinic connection for parameter values on one side of the hypersurface H. In Case 2, a saddle periodic orbit emerges from the homoclinic orbit. Its stable manifold is orientable or not, depending on the orientability of W u (0). In Case 4, there exist infinitely many saddle-type periodic orbits in any neighborhood of the homoclinic orbit. In fact, as argued in Ref. 44, there exist infinitely many horseshoes in any neighborhood of the homoclinic orbit 0 . When the connection is destroyed, finitely many of the horseshoes persist and hence it follows the existence of an infinite number of periodic solutions. The appearance or disappearance of horseshoes is accompanied by unfoldings of homoclinic tangencies of saddle-type periodic orbits and hence, strange repellers should emerge. 35,36,45 Readers can find more extended explanations about all these bifurcation results in Refs. 37 and 38.
Regarding codimension-two homoclinic bifurcations, we only pay attention to the inclination flip, orbit flip (OF), and Belyakov bifurcations because they are the only cases that we will discuss in the context of the Hindmarsh-Rose model. So, we distinguish the cases below.
Inclination Flip (IF): Assume all conditions in Case 2 except (H2), that is, we assume that the intersection between W cs (0) and W u (0) is non-transversal along 0 . Orbit Flip (OF): Assume all conditions in Case 2 except (H1), that is, we assume that 0 ⊂ W uu (0). Belyakov Point: Assume that the equilibrium point is a saddlenode (SN) with eigenvalues λ s and λ u with λ s < 0 < λ u . The eigenvalue λ u has geometric multiplicity one and algebraic multiplicity two.
To characterize different types of inclination and orbit flip bifurcations, we need to introduce the following ratios between eigenvalues: Note that α > β.
Bifurcation diagrams corresponding to IF and OF bifurcation points are quite similar and they can be described simultaneously. First, we observe that the hypersurface H of homoclinic bifurcation splits into two regions separated by a manifold of codimension-two homoclinic bifurcations. The orientation of the unstable invariant manifold at the equilibrium point reverses when such manifold is crossed. For either IF or OF bifurcations, there are three cases (see Fig. 1):

Inclination flip
Orbit flip We are only interested in Case C because the other two cases are not detected in our exploration of the HR model. Homoclinic flip bifurcations in Case C require additional generic assumptions. In particular, for inclination flips, we assume: Fig. 1), the homoclinic orbit does not lie in the unique smooth leading unstable manifold. (I3) If β < 1 2 α (region C 2 in the left panel of Fig. 1), there is a quadratic tangency between W cs (0) and W u (0) along the homoclinic orbit.
On the other hand, for orbit flips in Case C, we assume: Hypothesis (I2) [respectively, (I3)] makes sense in the region C 1 (respectively, C 2 ) depicted in Fig. 1. We do not extend in details about these two cases because they make no difference in the unfoldings. The essential distinction has to do with the way in which the unstable manifold approach the origin when it is followed along the homoclinic orbit by the forward flow (see Fig. 2 in Ref. 40).
There are two possible bifurcation diagrams in case C. In both cases, horseshoes exist in a region of the parameter space. We remark that chaotic regions have been observed in the HR model 46 connected with the infinite fans of period doubling and fold bifurcations of periodic orbits generated at these codimension-two points. Depending on how they are formed, cases C (in) and C (out) are distinguished (see Fig. 2). In both, infinitely many one-sided curves of secondary N-homoclinic orbits emerge for each N ≥ 2 from the flip point on the branch of primary homoclinic orbits (labeled hom in Fig. 2). These are homoclinic orbits that follow N times the primary one before closing up. Also, in both cases, the bifurcation diagram exhibits an infinite fan of bifurcation curves corresponding to period doublings and folds of periodic orbits. The horseshoe dynamics appear in between that cascade and the infinite fans of N-homoclinic orbits. In case C (in), shift dynamics and the homoclinic cascade are separated by the curve hom, whereas, in case C (out), the homoclinic cascade, the shift dynamics, and the fan of bifurcations of periodic orbits are located on the same side of the curve hom (see Fig. 2). A complete description of the bifurcation diagrams can be found in Refs. 37, 40, and 41.
Regarding Belyakov bifurcations, we remark that the hypersurface H of homoclinic bifurcation splits into two regions separated by a manifold of codimension-two homoclinic bifurcations. Saddles change from saddle-node type to saddle-focus (SF) type when such manifold is crossed. Additional generic conditions include global assumptions on the behavior of the invariant manifolds (see Refs. 37 and 43 for a complete description and particularly Fig. 14 in the second reference).
If λ s + ρ u < 0, a unique unstable limit cycle bifurcates from the homoclinic orbit (see Ref. 43). Otherwise, a two-parameter bifurcation diagram is quite similar to those in Fig. 2. Infinitely many one-sided curves of secondary N-homoclinic orbits emerge for each N ≥ 2 from the Belyakov point and they are tangent at the flip point to the branch of primary homoclinic orbits corresponding to saddle-focus. The bifurcation diagram also exhibits infinite fans of bifurcation curves corresponding to period doublings and folds of periodic orbits, but, in contrary to what is shown in Fig. 2, they accumulate on the branch of saddle-focus homoclinic orbits from both sides (see Fig. 14 Codimension-three homoclinic bifurcations have been studied in Ref. 40. In particular, transitions from Case A to Case B and also from Case B to Case C were discussed and conjectural bifurcation diagrams were provided. See also Ref. 42 regarding the case of the coalescence of resonances between eigenvalues with an orbit flip degeneracy. In both references, particular attention is devoted to the existence of homoclinic doubling cascades. Our study of the homoclinic phenomena in the HR model focuses on codimension-one and codimension-two bifurcations, but, as expected in a three-parameter study, higher codimension configurations do exist. For instance, coalescence between IF and Belyakov bifurcations and transitions from C 1 to C 2 in Fig. 1 (left) are expected in the HR model. Nevertheless, although these codimension-three phenomena have not been previously considered in the literature, they are out of the scope of this paper. Despite this, none of the scenarios considered in Refs. 40 and 42 have been detected in the HR model, but the bifurcation diagrams there proposed should inspire our future analysis of such configurations. These diagrams show pencils of codimension-one bifurcations connecting codimension-two bifurcation points. This is similar to what is shown in Fig. 6 in Ref. 20.

B. Slow-fast dynamics and fold-hom bursters
Equilibrium points in the full system (1) are given, after substituting the parameter values, by the intersection of the plane and the curve They do not depend on ε, but there can be one, two, or three equilibrium points depending on the values of parameters b and I. Projections of the plane (3) and the curve (4) on the plane (z, x) are illustrated in Fig. 3 for b = 2.7 and I = 2.2; see the brown colored straight line and the green-red colored Z-shaped curve, respectively. For these parameter values, there is a unique equilibrium point in the full system (1).
A detailed discussion about local bifurcations was given in Ref. 26. In particular, the description provided in Ref. 26 (Fig. 3) is similar to the information given at the bottom panels in our Figs. 5 and 6. As a reference, we use the bottom panel in Fig. 6 for the value ε = 0.08. For parameters in the purple region, there are three equilibrium points. Outside this region (at least in the range of parameters under consideration) there is only one equilibrium point that is attracting for parameter values on the green region until it undergoes a Hopf bifurcation (yellow line). The pale blue region correspond to saddle-focus (SF) equilibria with stability index 1, that is, equilibria where the linear part has eigenvalues λ s and ρ u ± ω, with λ s < 0 < ρ u and ω ̸ = 0. The transition from the pale blue to the red region means the change from SF to saddle-node (SN) equilibria (with stability index 1), that is, equilibria where the linear part has eigenvalues λ s , λ u , and λ uu such that λ s < 0 < λ u < λ uu , with λ s < 0 < λ u < λ uu . Note that λ u = λ uu for parameters on the borderline between the pale blue and the red regions. This transition is related to the existence of Belyakov bifurcations, which were described in Sec. II A. The difference between red and white regions-labeled SN1 and SN2, respectively-has to do with conditions on the eigenvalues, which are used to characterize specific cases of flip bifurcations (see Remark 2). In any case, both regions correspond to SN equilibria with stability index 1.
The Hindmarsh-Rose model is a prototypical example of a fast-slow system. The bifurcation diagram of the fast subsystem, obtained when ε = 0 is crucial to explain the dynamics when ε is small. 2 It should be remarked that each time that we refer to the fast subsystem (5), z is considered as an additional parameter. Fixing b and I, the model analysis provides two invariant objects: a curve of equilibrium points, with equations given in (4), and a manifold of limit cycles. As illustration, in Fig. 3, we show a partial bifurcation diagram of (5) with b = 2.7 and I = 2.2. The Z-shaped curve corresponds to equilibrium points: solid green lines correspond to stable equilibria, whereas dashed red lines correspond to unstable points. Note that the displayed curve corresponds to the projection of the curve with Eq. (4) on the plane (z, x). Stability along the lower branch is lost at a fold bifurcation point. There is also a second fold where the equilibria recover their stability to become again unstable when they undergo a Hopf bifurcation. The emerging limit cycles disappear in a homoclinic bifurcation to emerge again for lower values of z through an additional homoclinic bifurcation. This second family of limit cycles disappears at a Hopf bifurcation point, which is not displayed in the figure. We also show the maximum and minimum values of the x variable along the periodic orbits with solid blue lines. So, in general, we identify two invariant manifolds. On the one hand, the fast manifold M fast , also named spiking manifold, is given by the second family of attracting limit cycles of the fast subsystem (5) and, on the other hand, the slow manifold M slow , formed by the equilibrium points of the fast subsystem (5).
It follows from the Fenichel theory that for values of z where these manifolds are normally hyperbolic, they perturb to invariant manifolds M ε fast and M ε slow , which exist for ε small enough in the full system.
Bursting in the full system emerges because orbits repeatedly switch between M ε slow and M ε fast . An example of a bursting orbit with 5 spikes for ε = 0.01 is shown in Fig. 3. Top panel shows the bursting orbit projected on the plane (z, x) and superimposed on the picture of the fast-slow decomposition. The time series of the x component of the solution is displayed in the bottom panel. Note that the active regime begins close to a fold bifurcation of equilibria and finishes at a homoclinic bifurcation of limit cycles in the limit case. Due to this reason, following the Izhikevich 6 classification of bursting types, we say fold/hom bursting (also named square-wave bursting) to refer to the case illustrated in Fig. 3 In the literature, there are a lot of papers devoted to the study of the variation in the number of spikes that can be observed when one parameter is changed. Thus, plots similar to Fig. 4 are obtained (see also, for instance, Fig. 4 of Ref. 34), where the number of spikes in the neuronal response increases from two to six as a parameter is varied, and where each spike-adding transition is characterized by a strong increase in the L 2 integral norm of the orbit. By spikeadding process, we mean any mechanism leading to the formation of extra excursions around the tubular invariant manifold M fast (and therefore the addition of one spike to the bursting orbit).
In Fig. 4, we use the HR model to exemplify a process of spike-adding. We fix ε = 0.01 and I = 2.2 and let b vary as the continuation parameter of a periodic orbit. It is clear from the picture that a sequence of fold bifurcations (blue dots in the figure) is involved, giving rise to hysteresis phenomena and the appearance of bistability regions (in Fig. 4, we show two examples of coexisting stable periodic orbits). Although they are not shown, period-doubling bifurcations may also be present. As shown in Refs. 15, 20, 26, and 27, at least in the case of the HR model, all these bifurcations of periodic orbits are related to homoclinic phenomena.

III. ANALYSIS WITH ε FIXED
In this section, we begin our analysis by describing all the information provided by a selection of horizontal slices with the small parameter ε fixed. These selected slices will show us the different scenarios that we can find by changing ε, and it will help us later to develop a complete three-dimensional bifurcation diagram in the parameter space (b, I, ε) shown in Secs. IV-VI. Also, these twoparameter plots will show the connection of the spike-adding process with the "far-away" codimension-two homoclinic bifurcation points. Recall that the notation hom (n,n+1) was already introduced in Sec. I to refer to codimension-one homoclinic bifurcation curves.
As a first analysis, Figs All the ingredients that we need in our description of dynamical and topological changes are shown in Figs. 5 and 6. The displayed bifurcations are the following: black lines correspond to hom (1,2) bifurcation curves; red lines represent period-doubling bifurcation curves; yellow lines stand for Hopf bifurcation curves; red points are Belyakov bifurcation points and green and gray points represent, respectively, IF and OF bifurcation points. When displayed all together, the homoclinic bifurcation curves hom (n,n+1) are not distinguishable because for low values of ε, they are exponentially close and the largest is hom (1,2) , the one shown. Therefore, the IF and Belyakov bifurcation points corresponding to different homoclinic curves are superimposed (they are in different homoclinic curves but at a very small distance). The OF bifurcation points also correspond to several homoclinic curves (to be studied later), but they are clearly distinguished. In Fig. 7, we provide an alternative schematic view. Taking four representative values of ε, we show separately the homoclinic curves hom (1,2) , hom (2,3) , and hom (11,12) and some connected bifurcations. These figures illustrate the changes that can be expected in our global study and that we should explain.
In each lower panel of Figs. 5 and 6, the parameter plane is partitioned in different regions corresponding to different types of equilibrium points. As already explained in Subsection II B, this classification does not depend on ε. There is either a unique equilibrium point (purple region labeled 1EP and only displayed for ε = 0.07 and ε = 0.08) or three equilibrium points (3EP). In fact, we only need to pay attention to regions where the unique equilibrium point is a saddle-focus (region SF in the plots) or a saddle-node (regions SN1 and SN2 in the plots). Distinction between regions SN1 and SN2 has to do with two different cases for IF bifurcations characterized in Subsection II A. In particular, if a Case C of IF bifurcation is detected for parameter values on SN1 (respectively, SN2), hence eigenvalues correspond to the region C 1 (respectively, C 2 ) shown in Fig. 1 (left). Moreover, eigenvalues at the saddle-node point for parameter values in regions SN1 and SN2 correspond to region C in Fig. 1 (right), where the cases for the OF bifurcations are shown. In short, all IF and OF bifurcations are in Case C. Lower panels also display the curve hom (1,2) to understand all different types of homoclinic bifurcations: saddle-focus homoclinic orbits along sections contained in region SF and saddle-node homoclinic orbits along sections contained in regions SN1 and SN2. Several changes can be observed as ε increases. First of all, as we have already noted in Ref. 46, there is an evolution in the shape of the homoclinic bifurcation curves. For lower values of ε, the homoclinic bifurcation curves have a C-shape, with just one visible fold (as we can see in the case ε = 0.01). For intermediate values of ε, the C-shape transforms into a Z-shape, with two visible folds (see ε = 0.03). Last, for higher values of ε, the homoclinic bifurcation curves have no visible folds (ε = 0.07). As shown in Refs. 22, 23, 46, and 47, the C-shape is typical of the homoclinic bifurcation curves in the fast-slow regime. Another apparent change is the disappearance of some codimension-two bifurcations. Regarding IF points, when ε is small enough (for instance, ε = 0.009 18), there is only one IF point. When ε increases a little (ε ≈ 0.01), there are two IF points. When ε = 0.01 (see Fig. 5), the uppermost IF point is superimposed to the Belyakov point. For smaller values of ε, the role of the IF point is taken by the Belyakov bifurcation point. Besides, for large values of the small parameter (ε ≥ 0.02), there are no IF points. Obviously, these facts need a more detailed analysis provided by the three-parameter study done in Sec. IV as one may ask him/herself about codimension-three bifurcation points. Regarding OF bifurcation points, for ε = 0.015 (see Fig. 5), we show four OFs, one for each homoclinic bifurcation curve hom (1,2) , hom (2,3) , hom (6,7) , and hom (11,12) (there are more OF points on each curve but we just present one to show a scheme). For ε = 0.03 (see Fig. 6), only three OFs remain, due to the disappearance of the one on hom (11,12) . In fact, the complete homoclinic curve hom (11,12) disappears, together with the strip corresponding to 11 spikes per burst. For ε = 0.05, there are two OF points placed on hom (1,2) and hom (  (although there are some bands with bursting dynamics). Again, all these changes ask for a detailed three-parameter study. Recall that attending to the lower panels of Figs. 5 and 6, we can conclude that all OF and IF bifurcations are in Case C. This fact implies the birth of an infinite number of fold and period-doubling bifurcation curves emerging from these points, as well as infinitely many secondary homoclinic bifurcation curves with extra passages close to the equilibrium point (see Fig. 2).
The bifurcation diagrams in Figs. 5 and 6 also show the disappearance of the Belyakov bifurcation points. As ε increases, the distance between the two Belyakov points shrinks until they collapse; for ε = 0.08, there are no Belyakov bifurcation points. Lower panels help to understand how the Belyakov bifurcation points disappear. As ε increases, the homoclinic bifurcation curve has a smaller portion in regions SN1 and SN2. Note that the Belyakov bifurcation points appear when the homoclinic bifurcation curve intersects the borderline between regions SN1 and SF.
As it can be observed in the upper panels of ε = 0.018, 0.02, 0.03, qualitative changes in the period-doubling (PD) bifurcation curves occur for values of ε near to the value for which IF bifurcation points disappear (ε ≈ 0.0197). For ε = 0.015, we have plotted just one of the PD bifurcation curves emerging from each IF bifurcation point and for each one of the homoclinic bifurcation curves (in fact, the theory 37 regarding IF bifurcation points shows that infinitely many one-sided PD bifurcation curves emerge; see Fig. 2). A continuation of these curves in the plane (b, ε) shows that pairs of PD bifurcation curves are transformed into a single curve that persists for higher values of ε. This fact is a direct consequence of the disappearance of IF points where the pencils of PD and fold bifurcation curves are born. Therefore, the curves do not have a mechanism to finish and so they have to continue connecting both branches. Effects of this type have been already reported in the literature in other contexts (e.g., Refs. 48 and 49).
In order to summarize all the previous results, we show in Fig. 7 the complete global schemes with the different possibilities on the parameter plane (b, I) when the parameter ε changes. The schemes correspond to the results obtained for particular values of ε, but each bifurcation diagram is persistent, that is, it is qualitatively equivalent on any close enough horizontal slice. In the figure, we show a table in which each row corresponds to a certain transition from n to n + 1 spikes, while each column corresponds to a given value of ε. For each n and for each value of ε, we show the corresponding homoclinic bifurcation curve(s), the codimension-two homoclinic bifurcation points and some PD bifurcation curves. Color codes are those used in Figs. 5 and 6. When two adjacent boxes share the same diagram, we mean that the corresponding two cases are qualitatively the same. When a certain box appears crossed out, it means that there is no homoclinic structure for the corresponding transition in the number of spikes and for the given value of ε. This organization allows the reader to have a clear sight of all the different situations and to understand how the homoclinic structures vary as ε moves and different number of spikes are considered.
The first row of the table, i.e., the cases associated with 1 spike, has been already discussed. As it can be easily observed, the main difference between the case n = 1 (change from 1 to 2 spikes) and the other cases is that in the latter cases there is no longer a unique homoclinic curve for all values of ε, but two homoclinic curves exist for low values (this is the first time this fact is observed in the HR Chaos ARTICLE scitation.org/journal/cha model). Second, it is also important to note that the number and the type of codimension-two bifurcation points vary significantly with n. In the case n = 2, for all the values of ε, the codimensiontwo points present a similar situation to their analogs of 1-2 spikes. However, in the case n = 11, some of the codimension-two points that appear in the former cases do not exist (see, for example, the Belyakov points for ε = 0.009 18 and 0.015). Last, the case n = 11 reveals that the persistence of the homoclinic structure as ε increases depends on the number of spikes to which it is associated (see the fourth column, corresponding to ε = 0.08). This fact suggests the existence of a mechanism of disappearance of the global structures for large number of spikes when ε grows. All these numerical findings and hypothesis underlying these differences will be discussed in Secs. IV-VI. Note that all the previous discussions make clear that when dealing with fast-slow systems, the understanding of the mechanisms of creation and destruction of spikes requires studies in spaces of parameters which include the "small parameters." It is essential to have a global view of the bifurcations and Secs. IV-VI will stress the relevance of this goal.

IV. GLOBAL ANALYSIS CHANGING ε
As shown in Sec. III, a higher dimensional analysis is needed in the parameter space in order to explain the changes in the bifurcation diagrams observed in planes (b, I). In this section, we will discuss the three-dimensional structures associated to the different homoclinic bifurcation curves we have observed.
In Figs. 8-10, we provide bifurcation diagrams in the threeparameter space (b, I, ε). Codimension-one homoclinic bifurcations are shown in black, Belyakov bifurcations in magenta, IF bifurcations in green, and OF bifurcations in gray, as in previous pictures of this article. We have calculated curves of codimension-one homoclinic bifurcations with a step of 0.001 in the parameter ε using AUTO software, in order to visualize surfaces. For each case, the three-dimensional diagram is shown, as well as projections in the  planes (b, I) and (I, ε). These representations allow us to understand the mechanisms of appearance or disappearance of the different codimension-two bifurcation curves. It must be remarked that we have found difficulties for the continuation of OF bifurcation curves with AUTO in the HR model. For that reason, the continuation of OF curves is only partial in Figs. 8 and 9. In the parametric zones, where we have been able to obtain the OF points, we provide an interpolated curve in gray color. We conjecture, taking into account the points already calculated and the rest of bifurcation curves, that the full OF bifurcation curve in these two cases will be similar in shape to the IF curve. They will show a fold for large ε values, and for ε ↘ 0, they can continue or they can end in either a codimensionthree point (such as the IF curve in Fig. 8) or at one turning point of the homoclinic codimension-one curves when they have two components (such as the IF curve in Figs. 9 and 10 and the OF curve in Fig. 10). In any case, the numerical results show us a complete picture of the global dynamics of the system. Looking at the first two cases in Fig. 5, we observe how a IF bifurcation point appears close to the upper Belyakov point. If we observe now Fig. 8, we clearly see that it seems that the IF and Belyakov bifurcation curves collide at the numerically obtained parameter values, ε ≈ 0.009 189, b ≈ 3.102, I ≈ 4.713.
This "collision" would give rise to a codimension-three point that it is not studied in the literature, but it is out of the scope of this article. Besides, it is clear that, in the case hom (1,2) (Fig. 8), the Belyakov bifurcation points and also the IF bifurcation points disappear due to a folding of the bifurcation curve with respect to ε (the maxima we can observe in the 3D plots) of their corresponding bifurcation curves in the three-dimensional parameter space. Specifically, the Belyakov bifurcation curve has its folding point at ε ≈ 0.0748 and the IF bifurcation curve at ε ≈ 0.0197. In the case of hom (2,3) (Fig. 9), the Belyakov bifurcation curve presents a similar behavior to the hom (1,2) case. However, there is a very important difference in the way the IF bifurcation curve disappears. Note that curves forming the surface hom (2,3) have two disconnected components for (fixed) low values of ε. In addition, the system ceases to exhibit homoclinic connections in one of the regions in the parameter space where the geometry of the flow is the appropriate for the formation of IF bifurcations. This situation appears again in all the codimension-two curves in the case of hom (11,12) (Fig. 10). Therefore, we can observe a clear difference between hom (1,2) and all the other cases. This change in the topology of the homoclinic surfaces will be explained in more detail in Sec. V.
There is also another remarkable difference regarding the values of the small parameter for which each homoclinic surface disappears. In the cases hom (1,2) and hom (2,3) , it can be seen that the homoclinic curves clearly persist for all the values of ε we have studied, namely, up to ε = 0.08. Note that for larger values, we cannot consider the system as a fast-slow one. However, in the case hom (11,12) , the homoclinic surface has disappeared at ε ≈ 0.038. Using the SC technique, we discover band structures in the parameter planes with ε fixed, as shown in Figs. 5 and 6. Each band is associated to a given number of spikes per burst. The spike-adding process in fold/hom bursters was connected recently 25,50 with saddletype canards. 51,52 Besides, the necessary fold bifurcations of periodic orbits of the spike-adding process for hold/hom bursters were also recently connected with codimension-two homoclinic bifurcation points and also the homoclinic orbits experiment canard phenomena on one turning point of the homoclinic bifurcation curves. 15,27 Our numerical findings also support this idea, as they show clearly that the disappearance of a band corresponding to n spikes is linked to the disappearance of the corresponding homoclinic curves (surfaces) hom (n,n+1) . This is a quite important consequence of the three-parameter plots, as they explain the simplifications that are observed in the band structure of the fold/hom regime as ε increases, giving rise to burst phenomena with a small number of spikes (see in Figs. 5 and 6 how the number of color stripes decreases when ε grows). All the above mentioned features, together with the SC sweeps, suggest that the bigger the number n of spikes is, the smaller is the value of ε for which the corresponding homoclinic curve vanishes. Moreover, the numerical results show that the different homoclinic curves are stacked in a certain direction, being hom (1,2) the first one, providing an upper bound for "length and shape." The other homoclinic surfaces are disposed, exponentially close each other, as slabs in increasing order with respect to number of spikes per burst, but decreasing their size.
We have checked that Belyakov and IF bifurcation curves of different numbers of spikes overlap with each other in all the points in the (b, I, ε) where they coexist (they are exponentially close each other, like the homoclinic bifurcation surfaces). One can understand that the magenta (Belyakov) and green (IF) curves are placed in a fixed location in all the diagrams due to the requirements for their existence, and the existence or not of bifurcation points for some of the ε values depends if the corresponding homoclinic bifurcation curves (black curves) cut them. However, OF bifurcation curves corresponding to different numbers of spikes do not coincide with each other, and in fact they are quite far. This behavior is consistent with the role of OF bifurcation points in the spike-adding process as stated in Refs. 15, 20, and 27. What remains in the numerical tests is to reveal what is the aspect of the homoclinic surface in all cases, that is, if it is just a one leave surface or it has folds and it is a two (or more) leaves surface. This is in fact a relevant question as it will give the global structure of the homoclinic leaves. We are going to show the structure of isolas displayed by the different homoclinic bifurcation curves, once the parameter ε is fixed. We do not pay much attention to explain the transitions from n to n + 1 spikes on a given curve or surface (for details of this process, see Refs. 27 and 31) on both sharp folds of the isolas. Isolas are isolated closed curves of solution branches; hence, the curve is homotopic to a circle. In the literature, there are several examples of isolas of equilibria 53,54 or limit cycles. [55][56][57] Computing many isolas is tedious and requires an adequate strategy. For instance, in Ref. 53, the authors develop a strategy for locating families of isolas of equilibria. In this article, we focus on the detection of isolas of homoclinic orbits (see also Ref. 58) in the parameter space.
By performing sections on the surface hom (2,3) and using AUTO, with a large number of points and steps to guarantee some numerical precision in the computations, we have obtained the results given in Fig. 11. The pictures show codimension-one homoclinic isolas in the parameter plane (b, I) for ε = 0.03 (panel A) and ε = 0.07 (panel B). In the case ε = 0.03, the AUTO software is not able to connect one side of the isola, and adjusting different parameters of the software, just slight increments in the length of For ε = 0.03, the passage through the milder visible folds (compared with the sharp U-turns of both extremes of the isolas) of the homoclinic curve exhibit no bifurcations as the plots xz along the isola show (-3-to -4-and -5-to -6-). It is important to remark that taking the homoclinic orbits close to the values of the parameter where the continuation software stops for ε = 0.03, subplots -2and -4-, the different orbits show exactly the same behavior, with just small modifications (as it also shows the intermediate subplot -3-for one side). Therefore, it is perfectly logical to conjecture in this case that both sides of the curve are connected giving an isola, even more taking into account the results for ε = 0.07 where the isola is fully obtained. Note that in Ref. 27, the homoclinic isolas and the homoclinic organization were not detected as their main interest was the spike-adding and canard process of the homoclinic orbits on the lower-right sharp fold of the homoclinic bifurcation curve for ε fixed. In panel C (Fig. 11), we show two magnifications of the lower sharp fold of the isola for ε = 0.07. In these zooms, instead of plotting on the parametric plane (b, I), we use the plane with b and the AUTO norm L 2 to get a clearer image of the fold, showing two curves, and thus it illustrates one extreme of the isola.
In any case, the numerics only can give strong evidences of the existence of the isola. This fact is shown in the theoretical scheme shown in Fig. 12. The black curve is our conjectured isola (based in our numerical results), but, as the observed phenomena is on a small distance in the parameter space (the isola is very "thin," with a width about 10 −8 ), other options can be possible, like the existence of foldings in both sides and also some extra homoclinic codimensiontwo points, that is, two connected isolas, that are able to give rise to the folds (one option can be the dotted curve in Fig. 12). In any case, all of our numerical results show that it seems that we really have isolas, that is, the topological structure of the black curve in Fig. 12.
If one looks at the theoretical unfolding of the OF and IF codimension-two points shown in Fig. 2, there is an infinite fan of secondary codimension-one homoclinic bifurcation curves. None of the numerical simulations on the system (our studies in this article and in Refs. 15,20,and 46,and in Refs. 26 and 27 of other authors) show any of these bifurcations and any dynamical effect that can be related to them. This fact allows us (as also done in Ref. 27) to conjecture that the secondary homoclinics are inside the very thin homoclinic isola, and, therefore, it is not computationally possible to observe any of them. With these elements, we propose in Fig. 13 a theoretical scheme of the secondary homoclinic bifurcation curves and their connections (in a similar way as in Ref. 40) in the cases of having an even or odd number of pairs of codimension-two points.
As already remarked, it is apparent that there is an overlap between the different hom (n,n+1) bifurcation curves (in fact, they are exponentially close to each other as commented above), except for the higher values of ε where a slight separation can be observed. This separation of the curves occurs progressively as ε increases, and it can be appreciated for ε > 0.07. In Fig. 14, we show superimposed the three homoclinic isolas hom (1,2) , hom (2,3) and hom (11,12) for ε = 0.036 and ε = 0.07 to show that the isolas are outside of each other but exponentially close.

V. THEORETICAL SCHEME: THE HOMOCLINIC "MILLE-FEUILLE"
In Sec. IV, we have explored the three-dimensional parameter space of the HR model considering in detail the homoclinic structure. What it remains is to provide a complete theoretical scheme that connects all the basic ingredients of the spike-adding process in fold/hom bursters. That is, on the one hand, we have that in the parameter space the system experiments the spike-adding process far from the homoclinic bifurcations. On the other hand, the spike-adding process requires of two fold bifurcations to give rise to hysteresis phenomena and canards on one side to generate the extra spike (see Refs. 20,25,and 50). However, where are these fold bifurcation points generated? These points form bifurcation curves that Homoclinic isolas hom (1,2) , hom (2,3) , and hom (11,12) for ε = 0.036 and ε = 0.07 showing their relative position.
are born at codimension-two bifurcation points located on the "faraway" homoclinic bifurcation lines. All the bifurcation lines, in fact pencils of fold and PD bifurcation lines, are born, like the "pages-ofa-book" at the OF and IF points of the hom (n,n+1) curves as shown in Figs [15][16][17]. First, in Fig. 15, we show the different homoclinic surfaces. All of them are composed of one or two tubular structures.

ARTICLE scitation.org/journal/cha
As the number of spikes of the homoclinic orbit grows, we distinguish three types, either a tubular surface (hom (1,2) ) or two tubular surfaces connected (hom (2,3) , . . . , hom (k,k+1) ) or, finally, surfaces that disappear when ε grows (hom (k+1,k+2) , . . .). Note that Figs. 8-10 also illustrate numerically each one of these three types of surfaces. In the scheme, the different homoclinic surfaces are clearly separated from one another, but in the real parameter space, they are extremely close when ε is small, being organized in shape and size by the hom (1,2) surface. When ε is large, the separation becomes evident, showing that, indeed, these homoclinic surfaces have no contact point when ε > 0 (see Fig. 14).
If we take a section fixing the value of ε, we find three different situations, already partially described in Ref. 46, depending on the value of ε. When ε is large [O(1)], the slices just show a few homoclinic isolas corresponding to a small number of spikes and without visible folds. For intermediate values of ε, the isola corresponding to hom (1,2) has Z-shape with two visible folds. The other isolas complete a Z-shape, or not, depending on their length. Finally, for small ε, that is, in the generic situation when we are concerned with fast-slow systems, the principal isola for hom (1,2) has a C-shape with one visible fold. The curves corresponding to hom (n,n+1) , with n ≥ 2, split into two isolas also disposed in such a way that they are adapted to the C-shape of the principal isola. In this case, all the homoclinic curves have two components (isolas) but the first one hom (1,2) and all of them have folds with branches exponentially close to each other.
Due to the fact that, from a certain point of view, homoclinic surfaces are piled up one upon another, we refer to this conjectured global theoretical structure as the fold/hom homoclinic "mille-feuille" organization. Note that for fixed ε, we have a finite number of homoclinic curves, but the number of them grows as ε decreases. 25,33 Codimension-one homoclinic bifurcations that form each surface hom (n,n+1) must be understood as primary bifurcations. These surfaces contain curves of codimension-two homoclinic bifurcation: IF, OF, and Belyakov points. Emerging from these curves, there exist surfaces of bifurcation of periodic orbits: PD or folds, some of them involved in the spike-adding process. Also attached to these curves there are surfaces of secondary homoclinic bifurcations arising in the inner side of the surface, that is, separated from the surfaces of bifurcation of periodic orbits by the surface of primary homoclinic bifurcations [see case C(in) in Fig. 2]. Note that this scenario is covered by the classical unfolding theory of codimension-two homoclinic bifurcations. 37,38 We remark that these unfoldings have to be "glued" to the homoclinic surfaces given by the "homoclinic mille-feuille." Figure 16 illustrates the described scenario. Each of these curves of codimension-two bifurcations behaves as the "spineof-a-book" located on the homoclinic surface (like the "bookshelves" of a "bookcase"), whose "pages" consist of surfaces of bifurcations of periodic orbits and secondary homoclinic bifurcations. Plot 16.A provides the theoretical scheme of a homoclinic surface with the curve of codimension-two bifurcation points that form the "spineof-a-book" structure creating the pencils of surfaces of fold and PD bifurcations. In plots 16.B and 16.C, we show some numerical results illustrating such a theoretical scheme. The plot 16.B presents a projection of the homoclinic structure for three values of ε. Also, in plot 16.C, we see the global three-parametric view illustrating the theoretical scheme proposed in 16.A.
Finally, Fig. 17 illustrates the complete "mille-feuille" organization together with the "books" of bifurcation of periodic orbits. Now,

FIG. 16.
Theoretical and numerical illustration of the "spines-of-a-book" structure on the hom (1,2) homoclinic surface. Each of the curves of codimension-two homoclinic bifurcations is identified with the "spine-of-a-book" gathering "pages" of fold bifurcations, period-doubling (PD) bifurcations, and also (not showed) secondary homoclinic bifurcations. Panel A shows this theoretical model in the case of a "spine" of orbit flip (OF) points. Panels B and C show numerical results illustrating typical "pages" of one of these "books." In particular, panel B shows numerical slices of a "book" projected on the (b, I) plane. A three-dimensional view is given in panel C. Attached to each "spine," we see two "pages" of fold bifurcation and one "page" of period-doubling.

ARTICLE
scitation.org/journal/cha FIG. 17. Complete "mille-feuille" and "spines-of-a-book" theoretical structure. In panel A, we recall the unfolding of the bifurcation diagram associated to a OF bifurcation: there are pencils of PD and fold bifurcation of periodic orbits and also a pencil of secondary homoclinic bifurcations. In panel B, we see how these pencils are attached along a primary homoclinic curve. The isola has an exponentially small width d. Panel C illustrates a collection of isolas for a small value of ε. Finally, a three-dimensional scheme is provided in panel D. We see three "bookshelves" (homoclinic surfaces) and with some "books" (codimension-two points and the bifurcations generated) on them.
we can identify each layer of the "mille-feuille" with a "bookshelf " keeping as many "books" as "spines" of codimension-two homoclinic bifurcations it contains. So, we have a complete "bookcase" of bifurcations of periodic orbits. Moreover, we must notice that each surface in the "mille-feuille" has their own collection of "spines," that is, their own collection of "books." This figure gives an idea of how much entangled the bifurcations involved in the spike-adding process is. As illustrated in Fig. 17 (panel B), there are "pages" of the "books" involved in the spike-adding process. We remark that Fig. 17 provides a complete theoretical explanation of all the numerical findings obtained in this article (and in the literature).
Our conjectured theoretical structure permits to link the global three-parametric structure (the homoclinic surfaces) with the spikeadding phenomena that can be observed on parameter regions that are quite far from the homoclinic curves. In addition, if we use another set of parameters, we can also observe the fold/hom spikeadding processes, even without homoclinic bifurcations in the entire parametric plane. This is easily explained from Fig. 15, as if our parameters do not cut the homoclinic surface we cannot observe the homoclinic orbits themselves. But what remains are the fold and PD surfaces generated on the codimension-two points attached to the homoclinic surfaces, as shown in Figs. 16 and 17. Following with the "bookcase" analogy, this will be the case if we have "books" wider than the "bookshelves," and we observe it without seeing the bookcase. Obviously, our theoretical scheme is necessarily a partial one, as other bifurcations and phenomena may be present on the complete global picture, but it englobes all the current numerical and theoretical analysis in literature. This article provides new insights into the spike-adding process and the global parametric study of the Hindmarsh-Rose model. We hope that it may be applied to other fold/hom bursters, and this is part of our future work.

VI. CONCLUSIONS
In this article, we have presented a three-parameter study of homoclinic bifurcations in the canonical Hindmarsh-Rose neuron model when it evolves in the fold/hom bursting regime. We have introduced a new structure, the homoclinic "mille-feuille" connected with the fold/hom spike-adding process. Fold/hom bursting is found in numerous fast-slow models, and we expect that most of the findings of this article will be present in many similar problems.

ARTICLE scitation.org/journal/cha
Exploration of other fold/hom bursters is a goal for our future work, but a preliminary study as well as the theoretical scheme of the spike-adding process was introduced in Ref. 20. Our numerical analysis using different techniques allows us to conjecture the global theoretical homoclinic organization. There exists a "mille-feuille" structure of tubular-like homoclinic surfaces. Each of them corresponds to a transition where the homoclinic orbit increases the number of spikes by one, that is, taking the appropriate paths of parameters, one could observe in the phase-space how the orbits pass from n to n + 1 spikes for certain n. Moreover, as ϵ increases, the disappearance of a homoclinic surface associated to the transitions from n to n + 1 spikes means the "de facto" disappearance in the surroundings of the band of periodic orbits with n + 1 spikes. This structure provides a theoretical explanation of why there is not a regular fold/hom bursting regime with a large number of spikes when the small parameter grows. Moreover, due to the tubular structures, an analysis for fixed values of the small parameter gives rise to the appearance of isolas of homoclinic bifurcation points.
Note that previous relevant studies in the literature 15,26,27 focus their attention on the spike-adding and canard process of the homoclinic orbits on the lower-right sharp fold of the homoclinic bifurcation curve for ε fixed. The other sharp fold, the isolas, and also the complete bifurcation scheme where not identified and studied.
Located on each homoclinic surface, we find curves of codimension-two homoclinic bifurcation. These curves act as the organizing centers for the framework of fold and period-doubling bifurcations of periodic orbits, which is behind one of the main spike-adding mechanisms. The discovery of the global structure of orbit flip, inclination flip, and Belyakov bifurcations is one of our main motivations. Homoclinic surfaces can be compared with "bookshelves," where the "books" of bifurcation of periodic orbits are kept. Hence, curves of codimension-two homoclinic bifurcations can be compared with the "spines-of-a-book." The global structure (homoclinic "mille-feuille" + "spines-of-abook"), which is revealed in the three-parameter space, is a motivation for further study of higher codimension bifurcation points, which appear on the homoclinic bifurcation surfaces. In fact, the global structure we have uncovered gives clues about part of the bifurcations, which should be expected when dealing with such bifurcation points (and their connections, in a similar way as some codimension-three phenomena provides a global theoretical picture in Ref. 40). These relevant open problems are out of the scope of this article but they are part of our current research.