Wild oscillations in a nonlinear neuron model with resets: (I) Bursting, spike adding and chaos

In a series of two papers, we investigate the mechanisms by which complex oscillations are generated in a class of nonlinear dynamical systems with resets modeling the voltage and adaptation of neurons. This first paper presents mathematical analysis showing that the system can support bursts of any period as a function of model parameters. In continuous dynamical systems with resets, period-incrementing structures are complex to analyze. In the present context, we use the fact that bursting patterns correspond to periodic orbits of the adaptation map that governs the sequence of values of the adaptation variable at the resets. Using a slow-fast approach, we show that this map converges towards a piecewise linear discontinuous map whose orbits are exactly characterized. That map shows a period-incrementing structure with instantaneous transitions. We show that the period-incrementing structure persists for the full system with non-constant adaptation, but the transitions are more complex. We investigate the presence of chaos at the transitions.

Introduction. Neurons are excitable cells that communicate with each other through stereotyped electrical impulses, called action potentials or spikes. Because of the almost invariable shape of spikes, it is widely believed that the neural code is contained in the times at which spikes are fired and in the types of spike patterns fired. In addition to phasic responses (returning to resting potential after a few spikes in response to an input) and tonic firing (spiking repeatedly), other ubiquitous patterns of neuronal activity include bursting, in which several spikes are produced in rapid succession, followed by a quiescent phase lacking spikes, and mixed-mode oscillations, in which subthreshold oscillations alternate with active periods of one or more spikes. In this paper and the sequel, we provide a fine description of bursting patterns, spikeadding transitions, and mixed-mode oscillations in hybrid integrate-and-fire systems, a form of neuron model exhibiting a range of desirable features in a framework that is simple enough to allow for mathematical analysis. In particular, in this paper, we provide a novel rigorous demonstration of the existence of a period-incrementing cascade that extends away from the limit of timescale separation.
Hybrid integrate-and-fire systems combine nonlinear subthreshold dynamics accounting for the excitable properties of nerve cells together with discrete resets of the voltage, abstracting the complex yet stereotyped mechanisms of spike emission and reset that generally occur at a much faster timescale. These models follow a tradition introduced more than a century ago [37,9] and have attracted important interest of the computational neuroscience community in the past decade or so. Indeed, by decoupling the input integration and cellular excitability from the processes related to the emission of the spikes, they provide a mathematically convenient description of the neuronal dynamics that yields a precise representation of spike timing, conserves the main excitability properties of neurons, exhibits a wide repertoire of behaviors [25,27,7,67,32], can be precisely fitted to data [33], and allows for efficient computational implementation, which, for example, led to the first simulation of a network on the order of the mammalian brain size [29].
In this paper, we focus on bursting patterns and transitions between them involving period adding and chaos, while in the companion we move on to mixed-mode oscillations. We have two main motivations for studying bursting and associated transitions in hybrid integrate-and-fire systems. The first is that, in light of their utility and popularity, it is essential to harness the mathematical simplicity of hybrid integrate-and-fire systems in order to elucidate their mathematical properties in general. The second is that bursting is an important biological phenomenon in a variety of neural contexts, for which certain theoretical aspects have proved difficult to study analytically in smooth dynamical system models. As reviewed in [28], bursting is found both in brain areas ranging from the neocortex (e.g., in pyramidal neurons of layer 5) to the hippocampus, the thalamus and other subcortical areas. Bursting is hypothesized to contribute to many brain functions and states, as reviewed and detailed in [24,30,50]; for instance, it may enhance the reliability and flexibility of information transmission [36,35], support synchronization [4] in contexts including slow synchronized sleep rhythms [16], drive automated behaviors such as respiration and locomotion [42,43,40], promote hormone or neuromodulator release [63,65], and play a role in pathologies such as epilepsy [52] and parkinsonism [58]. Moreover, the number of spikes per burst may itself be functionally significant and has been linked to the phase or slope of input signals received by bursting neurons [36,35,61].
Because of both its biological significance and its mathematical complexity, bursting has been the subject of significant attention from the theoretical community (e.g., [11]). A key approach for the study of these behaviors takes advantage of the vastly different timescales between spike emission and input integration to develop a slowfast decomposition of the dynamics. This method has led to the identification of minimal models of smooth differential equations supporting various possible forms of repetitive bursting and the classification of the geometry of distinct types of canonical bursting scenarios in smooth dynamical systems [55,24]. While the geometric structure of bursts is now relatively well understood in this context, it remains quite difficult to characterize transitions between different bursting patterns and to analytically determine the existence of bursts with a prescribed period or number of spikes per burst in such systems. Indeed, smooth dynamical systems displaying bursting require at least 3 dimensions and elucidating the features of bursting orbits relies on detailed slow-fast analysis. In particular, changes in the number of spikes arising in each cycle of a periodic burst pattern occur through complex transitions such as the spike-adding mechanism, in which a modulation of the value of a parameter leads to the addition of one spike per burst cycle. Evidence for the existence of transient chaotic behaviors has been reported when the fast dynamics has specific features (e.g., in the case of a slow-fast transition induced by a homoclinic bifurcation [38] and in the case of a double homoclinic bifurcation [15]). Yet the difficulty of fully characterizing the global return mechanisms induced by a nonlinear flow of 3 or more dimensions generally has precluded detailed quantitative elucidation of the underlying sequence of orbit bifurcations hidden in the spike-adding transition.
A fundamental work of Rinzel and Troy [56] initiated the rigorous mathematical study of bursting and spike incrementing in the Belousov-Zhabotinsky reaction by reducing the analysis to the study of a one-dimensional map approximated by a discontinuous, piecewise linear map. A similar idea was used by Levi [39] to provide a geometric explanation of the period-adding phenomenon observed experimentally in periodically forced neon tubes. Levi initiated the use of the properties of circle map for this study, which proved to be a particularly fruitful technique. These two approaches have been adapted in applications related to bursting in neurons [51,34,31,2] and in this vein, a variety of discrete map-based models of bursting neurons have been developed that provide a versatile and easily implementable description of neuronal dynamics [60,59,44,41] but generally remain abstract from the biological viewpoint (see the review [23]).
The present manuscript applies the map-based approach to a hybrid dynamical system and provides a detailed analysis of spike adding and transitions between burst patterns. Previous work reported that these neuron models can produce bursts of various periods [49,67,72,7,71,32], and numerical evidence suggested the presence of an underlying period-incrementing structure [72], in which bursts of distinct periods arise, separated by chaotic intervals, as reset voltage is progressively varied. This manuscript undertakes the rigorous study of the spike adding in the hybrid dynamical system framework. In particular, we establish here a deep relationship between the dynamics of nonlinear integrate-and-fire models and unimodal maps of the interval, whose dynamics has been extensively studied. Our approach combines previous studies on piecewise continuous maps with notoriously rich dynamics together with a slow-fast analysis.
The paper is organized as follows. We start in Section 1 by describing the model and summarizing relevant properties of nonlinear integrate-and-fire models, in particular those related to the adaptation map governing the sequence of values of the adaptation variable across successive resets. In Section 2, we consider the limit where adaptation is very slow compared to the voltage, in which case we show that the map governing consecutive resets of the hybrid model converges to a piecewise continuous map and we explicitly demonstrate the presence of spike adding (as in [56] and subsequent works). In the limit of timescale separation, the transition from bursts with k spikes to bursts with k + 1 spikes is instantaneous and no chaos appears. In section 3 we prove the persistence of the period-incrementing structure related to spike adding when the adaptation variable is slow but not constant along trajectories (as predicted by numerical evidence in a general context in [53]). Finally, we investigate in Section 4 the behavior of the system at the transitions between bursts of different periods, analyzing the emergence of different forms of chaos and establishing the presence of period doubling events.
1. Background and main results. In this section, we review the properties of the class of neuron models studied in this manuscript and in the sequel and summarize the main results that we attain in this paper.
1.1. Nonlinear bidimensional integrate-and-fire neuron models. The class of nonlinear bidimensional models used in the present manuscript describes the excitable properties of the membrane voltage of the nerve cell, v, coupled to an adaptation variable, w, according to the ordinary differential equation where ε > 0 accounts for the ratio of timescales between the adaptation variable and the voltage dynamics; b > 0 represents the steady state ratio of adaptation to voltage, and can be seen as the coupling strength between these two variables; I is a real parameter modeling the input current received by the neuron; and F is a real function accounting for the intrinsic dynamical properties of the cell membrane, typically provided by leak currents together with spike initiation currents that provide excitability. Specific models may differ in the form of F ; in particular, F was assumed to be quadratic in [25] and exponential in [7], while in [67] a quartic model was proposed that had the capacity to support stable subthreshold oscillations (i.e., periodic orbits with no spike).
Mathematically, it was shown that qualitative properties of these systems do not strongly depend on the precise choice of the nonlinearity F as along as a few minimal properties are satisfied [68,72], which we assume to be the case.
Assumption (A1). The map F : R → R has the following properties: • it is regular (at least three times continuously differentiable); • it is strictly convex; • its derivative diverges at +∞, i.e. lim v→∞ F (v) = ∞, and has a negative limit at −∞ (possibly also negative infinite) satisfying: • there exist η, α,v > 0 such that F (v)/v 2+η ≥ α for all v ≥v. The first three hypotheses control the shape of the v-nullcline and hence constrain the possible fixed points of the system. The last assumption ensures that, at least for some initial conditions, the membrane potential variable v blows up in finite time, while the adaptation variable remains finite. Hence, the need to introduce an artificial arbitrary spike threshold is eliminated (contrarily to the case of the linear or quadratic adaptive models, for instance 1 , see [69]. The blow up of v is interpreted as the time of a spike. Following a voltage blow up at time t * , the voltage is instantaneously reset to a constant value v R and the adaptation variable is updated as follows: with γ ≤ 1 and d ≥ 0, corresponding to the effect on the adaptation variable of the emission of a spike. With this reset mechanism, it is not hard to show that the system is globally well-posed, i.e. that one can define a unique forward solution for all times and initial conditions. This property requires that spikes do not accumulate in time, which can be demonstrated as done in [72] for γ = 1. Note that in all models in the literature it is assumed that γ = 1. However, going back to the biological problem, spikes are not Dirac masses but stereotypical electrical impulses s(t) = 1 δt S( t δt ) where S(t) is the typical spike shape rescaled on the dimensionless interval [0, 1], and δt is the spike duration, generally small compared to the input integration timescale: 0 < δt 1/ε. The adaptation variable integrates this sharp impulse: The classical nonlinear integrateand-fire neuron of [7,26,67] corresponds to the limit δt → 0 and hence γ → 1. In this paper, we take γ = 1, while in the companion we do not.
The excitability properties of the system governed by the subthreshold system (1.1) were investigated exhaustively in [67]. It was found that all models undergo a saddlenode bifurcation and a Hopf bifurcation, organized around a Bogdanov-Takens bifurcation, along curves that can be expressed in closed form. The different curves are depicted in the parameter space (I, b) in Fig. 1.1. They split the parameter space into a region in which the system has no singular point (yellow region), a region in which the neuron has two singular points, one a stable steady state (orange and blue regions), and a region in which the system has two unstable singular points, one a saddle and one repulsive (magenta and pink regions). In the latter case, the stable manifold of the saddle is made in part of a heteroclinic orbit connecting the saddle to the repulsive point. The heteroclinic orbit can either (i) wind around the unstable point (pink region B), or (ii) monotonically connect to the repulsive point (magenta region A), depending on whether the eigenvalues of the repulsive point are real (A) or not (B). The companion paper will specifically focus on the dynamics of the system in case B.
This paper treats the case where the input is large enough so that the system has no equilibrium, represented by the yellow region in Fig. 1  While numerics suggest that our results do not sensitively rely on this assumption, working under this condition simplifies the analysis, since the absence of equilibrium of the subthreshold dynamics implies that all initial conditions lead to spiking. In particular, this assumption allow us to define on the whole real line the adaptation map [72] that describes how the value of the adaptation variable evolves over consecutive resets on the reset line {v = v R }. Definition 1.1. The adaptation map Φ associates to any value w ∈ R the value of the adaptation variable after the spike and reset of the trajectory starting at initial condition (v R , w). Rigorously, if (V (·, v R , w), W (·, v R , w)) is the solution of equation (1.1) with initial condition (v R , w) and if V blows up at t * , then Φ(w) := W (t * , v R , w) = γW (t − * , v R , w) + d. In this paper, we take γ = 1. This characterization obviously defines a unique value Φ(w) for any w ∈ R, since under assumption (A2) every trajectory blows up in finite time and since the value of the adaptation variable remains finite by assumption (A1), see [70]. A fine characterization of the shape and regularity of the adaptation map Φ has been provided and illustrates the central role of the value w * = F (v R ) + I of the intersection of the reset line with the v-nullcline [72]. Similarly, by w * * = bv we will denote the intersection of the reset line with the w-nullcline. We recall a few properties useful to the present analysis [72]: P1. Φ is increasing and concave on (−∞, w * ) (with Φ (w) < 0 for w < w * ) P2. Φ is decreasing and bounded below on [w * , ∞) and thus has an horizontal asymptote (plateau) at infinity, provided that lim v→−∞ F (v) < −ε(b + √ 2). P3. Φ is at least C 3 (more generally, C k if F is as well) P4. Φ has a unique fixed point in R P5. For all w < w * * , we have Φ(w) ≥ w + d ≥ w.
The map Φ plays a central role in the analysis of the dynamics of the system since its iterations Φ n (w) define the sequence of values of the adaptation variable after spikes, from which one can infer the spike pattern fired by the neuron and distinguish regular spiking, spike frequency adaptation, bursting or chaotic spiking [72]. In the present case, the introduction of the adaptation map will reduce the analysis of the model to the study of a one-dimensional continuous unimodal map and thus allow the use of a number of well-developed tools from the theory of discrete dynamical systems, which we shall exploit in the present manuscript.
Before we continue, it is important to note that the study of this map in the context of bidimensional integrate-and-fire models has seen recent developments. In particular, Foxall and collaborators established a relation between the adaptation map and transverse Lyapunov exponents to characterize the stability of spiking periodic orbits in nonlinear integrate-and-fire systems [19]. In [32], a generalized linear integrate-and-fire system was investigated via a similar map that is locally contractive, either globally or in a piecewise manner; conditions for spiking and bursting dynamics were established and bifurcations underlying transitions between solution patterns were studied. In the accompanying paper [57], we address the question of the dynamics of the system in the presence of two unstable equilibria of the subthreshold dynamics; in that case, the map Φ is no longer continuous, and the rotation theory for discontinuous maps is applied to characterize the dynamics. The present study focuses on characterizing a sequence of period-incrementing bifurcations associated with spike-adding transitions in bursting solutions. The main results of the present study are summarized below.

Summary of the main results.
Numerical evidence suggests that, as v R increases, the adaptation map Φ undergoes a sequence of period-incrementing bifurcations characterized by the presence of chaotic transitions (see [72]), as depicted in Fig. 1.2 in the case of the quartic model F (v) = v 4 + 2av with the standard parameter set (used throughout the manuscript except otherwise specified): We explore this structure in detail in the present manuscript.
To establish the presence of a period-incrementing structure, we start by considering in Section 2 the dynamics of the system in the limit of perfect timescale separation (i.e., ε → 0). In that case, we will show that a number of properties of the map Φ can be obtained from the analysis of a simple piecewise linear map where (v F , w F ) are the coordinates of the unique minimum of the graph of F + I. We assume that w * + d > p 0 . This map is particularly easy to analyze mathematically: it has a simple dynamics with a unique globally attractive periodic orbit with period p = w * −p0 d + 2 (where denotes the integer part). Since w * is an increasing function of v R , the map will thus shows a period-incrementing bifurcation structure as a function of both parameters.
We show that the adaptation map Φ converges towards Φ 0 as ε → 0 in the Hausdorff distance (Proposition 2.1 in Section 2) using standard slow-fast properties of the flow, and in the C 1 distance on any bounded interval not containing w * using more refined slow divergence estimates (Lemma 3.7 in Section 3). This allows us to show in Section 3 that the adaptation map Φ inherits the period-incrementing structure of Φ 0 for small enough ε. However, in contrast with the instantaneous bifurcations of Φ 0 , the period transitions are complex for ε > 0. The study of the dynamics in the vicinity of the transitions is performed in Section 4. We show that Φ is topologically chaotic for a wide range of values of the reset voltage, for most reasonable definitions of chaos (e.g. chaos in the sense of Devaney or Block and Coppel, all of them being equivalent for the adaptation map). Of course, topological chaos is not necessarily reflected in the iterates of the map for non-specific initial conditions. This leads us in Section 4.3 to concentrate on the occurrence of the clearly visible chaotic bouts occurring at the transitions in the numerically obtained bifurcation diagrams. To this end, we characterize the presence of metric chaos, corresponding to situations where the map features an invariant measure that is absolutely continuous with respect to the Lebesgue measure.
2. Period-incrementing in the singular limit. We start by characterizing the bifurcations of the adaptation map in the singular limit of extreme separation of timescales ε → 0. In this section, we denote the map by Φ ε to emphasize its εdependence. Φ ε also depends on v R , but we omit this dependence from our notation. The limit of Φ ε as ε → 0 relies on a number of geometric invariants provided by Fenichel theory. Of particular importance, we will consider the critical manifold C formed by the singular points of the fast v-dynamics and given as a graph over v by w = F (v) + I. Thanks to assumption (A1), we know that C has the unique fold point (v F , w F ) ∈ C with F (v F ) = 0. This fold splits the critical manifold C into two branches, C − defined for v < v F and C + for v > v F . For any point (v,w) on C − (resp. C + ),v is an attractive (resp. repulsive) singular point for the fast dynamicṡ v = F (v) −w + I (considering w as a parameter,w).
By Fenichel theory, we know that any compact submanifold for small values of ε, into a (non-unique) normally attractive invariant manifold for the flow of (1.1) lying in a O(ε)-neighborhood of C − (for the Hausdorff distance). As usual, the attractive slow manifold C − ε is defined as the unique invariant manifold that is asymptotic to C − when taking the limit v → −∞ and extended by the flow for positive time. Hence, any compact submanifold is normally repulsive for system (1.1) and O(ε)-close to C + . The slow manifolds, together with the points and notation that will arise in our analysis, are depicted in Figure 2.1.
Since we are interested in the bursting regime, we shall assume that: . For ε small enough, for any v R , trajectories emanating from sufficiently large w values along the reset line accumulate on the stable slow manifold C − ε and hence p ε is simply the value of the adaptation variable after reset for an initial condition on C − ε . Thus, beyond its linear dependence on d, the value of p ε depends on ε but is independent of 1. An example of the adaptation map for small ε, annotated with notation used throughout the paper. Left : attractive (green, C − ε ) and repulsive (red, C + ε ) slow manifolds of the subthreshold system as perturbations of the critical manifold C (black curve). As stated in the main text (in order from largest to smallest w): w + ε denotes the w-coordinate of the upper intersection We start by showing that Φ ε converges, in a suitable sense, towards the piecewise linear map Φ 0 given in (1.6), with p 0 = w F + d. Indeed, for ε small enough, one notices heuristically that: • for any w ≤ w * , the value of w remains almost constant during the whole trajectory while v blows up, implying that the value of w at reset approaches w + d as ε → 0; • for any w > w * , the trajectory approaches the stable slow manifold (which is very close to the v-nullcline) and follows it until it fires, and thus the adaptation variable approaches p ε , which in turn converges toward p 0 when ε → 0. To rigorously make sense of the convergence of the family of continuous maps Φ ε towards the discontinuous piecewise linear map Φ 0 , we will rely on the notion of Hausdorff distance between the graphs of functions. In detail, let us denote by G(Φ ε ) the graph of the map Φ ε , i.e. {(w, y) ∈ R 2 : y = Φ ε (w)}, and G(Φ 0 ) the graph of the function Φ 0 augmented with the vertical segment at the discontinuity point w 0 : The Hausdorff distance between the sets X := G(Φ 0 ) and Z := G(Φ ε ) is defined by: where d(x, z) denotes the Euclidean distance between the points x, z ∈ R 2 . We show the following: We split the proof into three steps: (i) in Lemma 2.2, we show the C 0 -convergence of Φ ε towards Φ 0 for w < w * , (ii) in Lemma 2.3, we establish C 0 -convergence on w > w * + δ for some δ > 0, and (iii) we complete the proof by showing that as ε → 0, the value of δ can be chosen arbitrarily small. For the proof, we will use the notation w * v R rather than just w * to emphasize the fact that this point depends on v R . Note that we will strengthen this result in Lemma 3.7 by showing that this convergence is also true in a C 1 sense, i.e. that away from w * , the differential of the adaptation map converges to the differential of Φ 0 as ε → 0.
Proof. The proof proceeds by following the trajectories of (1.1) in the phase plane. For w < w * the solution starting from (v R , w) remains below the v-nullcline, thus v is stricly increasing along the whole trajectory. It is thus easy to show that we can express the trajectory of (1.1) as W ε = W ε (v; v R , w) for v ≥ v R , and where W ε is the solution of the following differential equation: (2. 2) The distance we aim at evaluating is thus given by: which is finite because of the integrability of u/F (u) at infinity (Assumption (A1)). The righthand side thus provides the uniform bound announced for , withv the value at which the trajectory W ε crosses the w-nullcline, and the part v ∈ [v, ∞): (2. 3) The second term is handled similarly as the case w < w * * . As for the first term, using the fact that v R > v F and v ∈ [v R ,v], we readily obtain: and thus conclude: In both cases, we find an upper bound proportional to ε, and these bounds are continuous with respect to v R . Hence, given δ > 0, we can find ε sufficiently small to obtain that for every Moreover, p ε → p 0 as ε → 0. Proof. We start by proving the convergence of p ε using the expression it follows that for any δ > 0, there existsε small enough so that for ε ≤ε It thus suffices to establish that the integral term vanishes as ε → 0, i.e. that for ε small enough we have which indeed follows as a simple consequence of the integrability at infinity of u/F (u) as in the proof of Lemma 2.2.
With this result in hand we now prove the convergence of Φ ε . Note that the backwards trajectory from For v R > v F and δ > 0, there existsε > 0 sufficiently small so that for any ε <ε, w * + δ/2 >w F . It follows that for any ε ≤ε and w ≥ For t > t ε , V ε is strictly increasing and the trajectory remains bounded between the v-nullcline and the forward trajectory from (v F , w F ). In particular, The same arguments as used above now apply to show that Φ ε (w) is within a δ/2-neigborhood of Φ 0 (w) = p 0 . Again, as the manifold C − ε and the v-nullcline do not depend on v R , this estimate holds uniformly in v R for v R varying in a compact interval.
[Proof of Proposition 2.1] Lemmas 2.2 and 2.3 taken together prove uniform conver- ; therefore, these parts of the graphs G(Φ 0 ) and G(Φ ε ) converge also in the Hausdorff metric. Moreover, the Hausdorff distance between the two maps for and hence its distance to the graph G(Φ 0 ) does not exceed δ. Altogether, we thus have: We emphasize that the limit (in Hausdorff distance) of Φ ε at ε = 0 does not correspond to a reset map for ε = 0 since the adaptation variable for the bidimensional system (1.1) with ε = 0 is not defined for w > w * . Indeed, the trajectories with initial condition (v R , w) with w > w * will simply converge towards some point (v, w) on the critical attractive manifold C − and therefore will not blow up and fire a spike. Thus in the statement of Proposition 2.1, Φ 0 is not the adaptation map for ε = 0, but just the limit of maps Φ ε for ε → 0. Now that the limit Φ 0 of Φ ε has been derived, we investigate the dynamics of Φ 0 and establish that this map exhibits a period-incrementing phenomenon.
Proposition 2.4. For any v R , the map Φ 0 has a unique periodic orbit, which is globally attractive and has period p given by: With the increase of v R and hence w * , the period of this orbit is incremented by 1 at each point w * = p 0 + (k − 1)d, k ∈ N. The map thus displays a period-incrementing structure with instantaneous transitions.
Proof. The dynamics of the piecewise linear map Φ 0 is particularly simple. The only possible fixed point of the map is p 0 . As soon as w * ≥ p 0 , this fixed point no longer exists, and the system displays a unique periodic orbit that contains p 0 : where k ≥ 2 is the smallest positive integer such that w * < p 0 + (k − 1)d and thus the orbit is periodic with period k.
It is easy to see that this orbit is globally attractive. Indeed, no orbit can be fully contained in the interval (−∞, w * ) as Φ is above the identity line there. Thus any trajectory will eventually be mapped to p 0 and absorbed by the periodic orbit O p0 given in (2.5). In other words, every orbit of this simple system is eventually periodic with period p. Figure 2.2 illustrates the shape of the map Φ 0 and displays its periodic orbit as well as the associated bifurcations occurring under variation of w * or equivalently v R .
3. Persistence of period-incrementing in the non-singular case. We work under assumptions (A1), (A2) and (A3) and build upon results obtained in the limit ε → 0 to show that bidimensional integrate-and-fire neurons governed by (1.1) with ε > 0 undergo a period-incrementing cascade as v R is increased. The main distinction between the singular limit and the system with ε > 0 is the fact that the adaptation map Φ ε is continuous and therefore the transitions from p-to (p + 1)-periodic behavior are not instantaneous, leaving room for much richer dynamics. Although Φ ε is continous and unimodal, and the values of w * and Φ(w * ) increase with v R , the distance Φ(w * ) − w * does not necessarily vary monotonically with v R . Therefore, the adaptation map does not follow the same bifurcation pattern as the logistic map.
Note that henceforth in this section, to lighten the notational burden, we will omit the index ε indicating ε-dependence except when particular emphasis on ε is important.
3.1. Attracting periodic orbits. It is known that necessarily, if Φ(w * ) ≤ w * , then there is a globally attracting fixed point w f ∈ (−∞, w * ] (see [72]). In this section, we shall thus concentrate on the case where Φ(w * ) > w * . In this case the unique fixed point w f of Φ is located in (w * , Φ(w * )) and one can justify that the interval [Φ 2 (w * ), Φ(w * )] is invariant and every trajectory enters this interval after at most a few iterates. Therefore, if additionally Φ 2 (w * ) ≥ w * , then the dynamics is trivial: (3.1) We will often make the assumption that the fixed point w f is unstable, i.e. when Φ (w f ) < −1. In this case, is well-defined.
Remark. Let us already mention that in Lemma 4.1 below we show that the assumption (3.1) holds for ε sufficiently small.
Before we study the periodic solutions of Φ, we need one more simple but useful result.
We now provide sufficient conditions for the existence of attracting periodic orbits. We associate to any k-periodic orbit not containing the point w * a signature (or itinerary) as follows: if the orbit contains m points smaller than w * followed by n = k − m points larger than w * , its itinerary is denoted L m R n .  Proposition 3.2. As discussed above, we assume that w * < w f < ξ < Φ(w * ) and Φ 2 (w * ) < w * . If, moreover, Φ 3 (w * ) < w * , then let k ∈ N be defined as If there existsw ≥ ξ such that (see Fig. 3.1):

3)
then Φ admits an asymptotically stable k-periodic orbit, with itinerary L k−1 R 1 , attracting the orbit of w * . Moreover, there is no other periodic orbit fully contained in the set (−∞, Thus interval [w, Φ(w * )] is Φ k invariant and contains a fixed point for map Φ k . Because of assumption (3.3), the Φ-orbit of this k-periodic point admits k − 1 points to the left of w * and one point to the right of ξ (where the derivative of Φ is smaller than 1 in absolute value). Hence, the derivative Φ along the orbit is negative and strictly greater than −1 on the account of Lemma 3.1, which proves that the associated fixed point of Φ k is stable. Now suppose that there is another periodic orbit fully contained in the set (−∞, w * ]∪ [w, ∞). Then this orbit necessarily contains points belonging to (−∞, w * ] as well as some points of interval [w, ∞). But the points of this orbit in [w, ∞) are necessarily within the interval [w, Φ(w * )], which is invariant for Φ k . For any w ∈ [w, Φ(w * )] and 0 < l < k, Φ l (w) ∈ [w, Φ(w * )]. Thus, any periodic point in (w, Φ(w * )) would necessarily be of period nk for some n ∈ N. Since Φ nk is a contraction on (w, ∞), Φ nk admits a unique fixed point in (w, Φ(w * )). This fixed point is necessarily the unique fixed point of Φ k in this interval that we identified before.
This result is quite general and will prove particularly useful for exhibiting the bifurcation structure of Φ for ε small enough: it will guarantee the existence of an attractive k-periodic orbit. It also allows proving that actually most points, in a certain sense, are attracted by this k-periodic orbit: Proof. The proof is elementary once we note that every point w ∈ [Φ 2 (w * ), Φ(w * )]\ H will after (at most) a few iterates enter the interval [w, Φ(w * )], and therefore any w ∈ [Φ 2 (w * ), Φ(w * )] \ H will be eventually attracted by the k-periodic orbit characterized in Proposition 3.2.
Remark. The assumption (3.3) for k ≥ 3 of Proposition 3.2 implies in particular that As proved below in Theorem 4.7, this inequality implies the existence of periodic orbits of all periods as well as chaos, at least on a compact subset of the dynamical core [Φ 2 (w * ), Φ(w * )].
Given these results, a natural question is whether the attracting k periodic orbit in the Proposition 3.2 is the only attracting periodic orbit of the map Φ. We remark than in general for unimodal maps, this need not be the case. However, we know that any other such an attracting periodic orbit will not attract the critical point w * and will not have a point in the interval [w, Φ(w * )].
To go beyond this result and address the question of the uniqueness of the stable periodic orbit, we would like to make use of the theory of negative Schwarzian derivative maps. We recall that: Definition 3.4. The Schwarzian derivative of a C 3 interval map f : I → I at x ∈ I such that f (x) = 0 is given by: Moreover, if a unimodal function f has a unique critical point c in I (i.e. f (c) = 0 and f (x) = 0 for x = c) and (Sf )(x) < 0 for all x = c, then we say that f is an S-unimodal map or, equivalently, that f has negative Schwarzian derivative, which we denote Sf < 0. In our case SΦ is well-defined everywhere except for w * , which is the only critical point, as we prove later in Theorem 4.3.
Corollary 3.5. Under the assumptions of Proposition 3.2 suppose that SΦ < 0. Then the k-periodic orbit established by Proposition 3.2 is the unique attracting periodic orbit of Φ.
The above follows immediately from the Singer Theorem (see e.g. [14]) as the S-unimodal map f : I → I with no attracting periodic points in ∂I can have at most one attracting periodic orbit, i.e. the one which attracts the critical point.
Since we aim to keep our analysis general, we do not restrict our choice of F to achieve SΦ < 0. Nonetheless, this property certainly holds for certain choices of F satisfying the general regularity properties that we require and for certain corresponding parameter sets. Thus, we will assume SΦ < 0 at some points to show the interested reader that when this condition holds, some of our results can be immediately strengthened.
3.2. Period-incrementing cascade. Proposition 3.2 remains quite abstract at this level of generality, and a natural question that arises is to characterize the parameter sets for which the proposition applies, and to what extent the conditions of this proposition are satisfied. We will prove that as ε → 0, the proposition applies for almost all reset values v R and as v R increases provides a sequence of asymptoticallystable k-periodic orbits with incrementing k that account for the numerically observed period-incrementing behavior (e.g., Figure 1.2).
Theorem 3.6. [Period-incrementing] For any integer N > 3, there exist ε > 0 and a sequence {J k } N k=3 of ordered intervals J k of reset values v R (i.e. v R k < v R k+1 for any v R k ∈ J k and v R k+1 ∈ J k+1 ) such that for any ε ≤ε and v R ∈ J k , k = 3, ..., N , the adaptation map Φ ε has an asymptotically stable k-periodic orbit with itinerary L k−1 R.
Furthermore, for any ζ > 0, we can pickε small enough so that for every ε ≤ε and any v R ∈ J k with k ∈ {3, · · · , N }, the set H ε of initial conditions w that might not be attracted by the k-periodic orbit of Φ ε has Lebesgue measure smaller than ζ.
Remark. This phenomenon is clearly visible in our numerical simulations. Indeed, while with the increase of vR the period-incrementing sequence persists, the intervals of values of vR associated to complex and seemingly chaotic transitions between intervals of periodic behaviors diminish significantly; see Fig. 1

.2.
Before we prove Theorem 3.6, we need a preliminary result relating Φ ε to Φ 0 . The results of Section 2 ensure that Φ ε approaches Φ 0 in the Hausdorff distance. The object of the following Lemma is to establish the C 1 -convergence of Φ ε to Φ 0 away from w * . Henceforth, we introduce the notation w * v R to emphasize that the value of w * is determined by the choice of v R .

Lemma 3.7. Given any bounded interval
Proof. The proof will use similar methods as that of Proposition 2.1. First we prove the convergence in the interval (−∞, w * v R − ν]. In any interval of this form we have (Φ 0 ) (w) = 1, while (Φ ε ) (w) satisfies by an application of Peano's Theorem (see e.g. Theorem 3.1 and Corollary 3.1 in chapter 5 of [21]). We want the above to be close to 1 and thus we need to show that the absolute value of the integral above is close to 0. But this in turn is equal to the sum of the integrals wherev denotes as before the value of v at which the solution W ε (v; v R , w) intersects the w-nullcline (for initial conditions w < w * * v R , while for w ≥ w * * v R we can just takȇ v = v R ). Althoughv depends on ε and v R it can be always overestimated asv < w * v R 2 /b. Thus for the first integral above we have where L 1 is the value of the integral of 0 < F (u)+I −w over the interval where L 2 denotes the value of the (convergent) integral of 1/(F (u) + I − bu) on the interval (v R1 , ∞). In this way we have obtained the desired estimate independently of To show the convergence for w ∈ [w * v R + ν, ∞), we notice that on this domain, the adaptation map can be expressed as the composition of two maps, where Θ ε assigns to w > w * v R the pointŵ on the reset line v = v R , below the vnullcline, such that the trajectory (V ε (t; v R , w), W ε (t; v R , w)) crosses the reset line at this point before spiking. Thus (3.6) The repulsive slow manifold C + ε (prolonged by the flow) intersects {v = v R } both above and below the v-nullcline. Having computedε for the first part of the proof, we can further assume thatε is so small that ε ≤ε implies that for every As a result, we can conclude that Θ ε (w) intersects {v = v R } below the lower intersection of C + ε with this line. Since the first factor (Φ ε ) (Θ ε (w)) in the multiplication in the formula (3.6) is always non-negative and bounded by 1, in order to show that (Φ ε ) (w) is small, it suffices to show that |(Θ ε ) (w)| ≈ 0 for small ε. For this purpose, consider initial conditions (v R , w 1 (0)), (v R , w 2 (0) := w 1 (0) + ∆w) with w 1 (0) > w * v R + ν, such that for both orbits, v (0), w (0) < 0. Denote trajectories from these initial conditions by Now, after these intersections, γ 1 is bounded between C − ε , namely the invariant slow manifold that perturbs from the attracting branch of the critical manifold for 0 < ε 1, and the critical manifold itself, while γ 2 is bounded between C − ε and γ 1 . Each trajectory crosses the v-nullcline and returns to intersect Σ a second time. Since the trajectories only traverse an O(ε) distance between intersections with Σ, the flow box theorem ensures that their w-coordinates differ by c 2 (ε)∆w for an O(1) constant c 2 (ε) at the second intersection.
Finally, we invoke the slow divergence integral [18,12] to quantify the change in the distance between the trajectories' w-coordinates as they evolve from Σ back to {v = v R } below the v-nullcline (and C + ε ), close to C − ε . Let (v Σ , w − ε ) denote the intersection of C − ε with Σ and W (v; v Σ , w − ε ) the w-coordinate of C − ε expressed as a graph of v for v ≥ v Σ . In our case, the slow divergence integral is given by dv.
Assume that v R is such that w * * < w F . Then we can take ε sufficiently small such that C − ε is bounded away from the w-nullcline between Σ and {v = v R }. Then this integral is negative and O(1), such that we have an O(1) exponential contraction in w from Σ to {v = v R }. Hence, (Θ ε ) (w) can be made arbitrarily small by shrinking ε, as desired. This completes the proof for w * * < w F . If w * * ≥ w F , then yet another section is needed, say at {v =ṽ R } with w * * v R < w F . The previous arguments bound the expansion from {v =ṽ R } to {v = v R } and hence the desired result still holds.
Proof. [Proof of Theorem 3.6] The proof relies on the fact that for the map Φ 0 it is relatively easy to satisfy the assumptions of Proposition 3.2; given this observation, we can then exploit the closeness of Φ ε > 0 to Φ 0 .
Indeed, given N > 3, there exists a sequence J 0 k , k = 3, ..., N , such that for any v R ∈ J 0 k the orbit of p 0 under Φ 0 is k-periodic, with itinerary L k−1 R. Since the set N k=3 J 0 k is bounded, we conclude that for any δ > 0 there existsε > 0 such that for any ε ≤ε, for any k ∈ {3, ..., N } and any v R ∈ J 0 k our convergence results imply: We henceforth assume δ < min{1, d/(N + 1)}, thus the above implies that |(Φ ε ) (w)| < 1/2 < 1 for every w ∈ [w * v R + δ/2, ∞) (and every choice of v R ∈ J 0 k ). Moreover, for any δ < 1 chosen, one can always, if necessary, instead of taking "large" intervals J 0 k (of the length equal to d), take their subintervals (which will be for simplicity denoted again J 0 k ) so that with any v R ∈ J 0 k we have which implies that Now having ε ≤ε, for any k ∈ {3, ..., N }, let J k := J 0 k . Note that by assumption δ < d/(N + 1) < 2d/5, and therefore the interval (3.2) for particular values of ε and v R ), and alsow With the above choice ofw v R , We conclude that the v R -dependent maps Φ 0 and Φ ε satisfy the assumptions of Proposition 3.2 for any v R ∈ J k , where k ∈ {3, .., N }, and the maps Φ ε have k-periodic orbits of signature L k−1 R, asymptotically stable. Notice that all the above reasoning stands for any δ ≤ min{1, d/(N + 1)}. We now prove the second statement. Choose ζ > 0 arbitrarily and select any k ∈ {3, ..., N } and v R ∈ J k . We define and introduce the sets which are well-defined for ε sufficiently small. As from the first part of the proof can be arbitrary, we can assume thatw v R < w * v R + δ, so that diamA 1 < δ. Simultaneously, we see that the (which corresponds to choosing δ smaller), the further (to the left) from w * v R is the right-endpoint of interval A 2 . Let B 2 < w * v R be the right-endpoint of interval A 2 . Then: Since we consider finite set of k values (k ∈ {3, .., N }) and sets J k are bounded, there exists a constantM ≥ M that is independent of the choice of k, v R ∈ J k and any ε ≤ε withε fixed in the first part of this proof for δ ≤ min{1, d/(N +1)}. Note in particular that taking ε smaller (e.g., by further lowering δ) makes the derivative (Φ ε ) (w) closer to 1, thus further from 0, for any fixed w < w * v R and therefore reinforces the above estimates.
Thus we have the explicit expression which yields that for every v R ∈ J k and k ∈ {3, 4, .., N }, the Lebesgue measure of H is bounded:

Chaos between period-incrementing transitions.
We have showed that the bifurcation diagram of the non-singular system shows a period-incrementing structure, in the sense that there exists a sequence of disjoint ordered intervals of values of v R for which the adaptation map features attractive periodic orbits of incrementing periods. In this section, we focus on the phenomena arising between two intervals J k and J k+1 , i.e. at the transition between periodic orbits of periods k and k +1. Chaotic period-incrementing transitions are expected in continuous maps well approximated by discontinous piecewise linear maps featuring pure period-incrementing (see e.g. [54]). Numerically, we observe chaos in transitions from period 2 to period 3, period 3 to 4, and possibly for additional period-incrementing transitions, preceded by one or a  Fig. 2.2 (we superimposed the plot of the map F (v) + I in orange to emphasize this similarity). As ε is increased, the map Φ slowly deviates from Φ 0 ; the associated bifurcation diagrams conserve the overall period-incrementing structure, but with larger transition regimes characterized by the presence of chaos. We also note that as v R increases, the bifurcation structure becomes more similar to the singular limit.
few period-doubling bifurcations, provided that one takes ε not too small (see Figure  3.2). As a more detailed example, in Fig. 4.1, we numerically illustrate the transitions from bursts of period 2 to 3. We observe that the period-2 orbit loses stability through a period-doubling bifurcation, yielding a period-4 orbit with points that progressively approach the region of instability where Φ is strictly smaller than −1; this period-4 orbit again loses stability, seemingly through a period-doubling bifurcation, and very rapidly progresses into a chaotic trajectory before suddenly stabilizing on a period-3 orbit. In fact, as we will see, chaos is present between any two intervals J k and J k+1 , k ≥ 2, but for greater values of v R and hence larger k, the chaotic transitions are more abrupt and therefore less visible in numerical simulations. Our main results of this section can be outlined as follows: Topological chaos (including the existence of periodic orbits of all periods and sensitive dependence on initial conditions) occurs for all parameter values v R big enough (Theorem 4.7). Moreover, at the transition between periodic orbits of types L k−1 R and L k R, one can expect a positive measure set of parameter values v R for which the system is strongly chaotic, with positive Lyapunov exponent and, presumably, with absolutely continuous invariant probability measure.
Before we proceed, let us make also one simple remark concerning the observed period-doubling bifurcations. Suppose that for some value v R , Φ has an attracting periodic orbit O of type L k−1 R (e.g., when the assumptions of Proposition 3.2 hold). Then for arbitrary w 0 ∈ O v R , we have −1 < (Φ k ) (w 0 ) < 0. Consequently, the orbit can lose its stability only when the derivative of Φ k at points on the orbit reaches −1. This observation allows us to predict the presence of period-doubling bifurcations in the transition to L k R orbits. , for values of v R spanning the period-incrementing transition between bursts with 2 and 3 spikes. A period-2 orbit undergoes a period doubling giving rise to a period-4 orbit, which itself loses stability, yielding chaotic spiking, before the system stabilizes on a period-3 orbit.
The purpose of this section is to rigorously explain why we necessarily observe chaotic behavior in such families of continuous maps and to describe precisely the chaotic properties of the adaptation map in corresponding regions of parameter space, using different notions of chaos. To this end, we will need to invoke some powerful results on the dynamics of unimodal maps, which, although today well-known to the specialist in the field, are certainly non-trivial. In the first part 4.1 we establish some general properties of the adaptation map such as uniqueness of the critical point (Theorem 4.3) and its non-degeneracy (Theorem 4.4). Next, in Section 4.2 we refer to various notions of topological chaos and show that for almost all parameter values v R the map Φ exhibits topological chaos, with any of the reasonable definitions of chaos, such as e.g. chaos in the sense of Devaney or Block and Coppel, all of them being equivalent for the adaptation map. However, since topological chaos is not always reflected in the observed behavior of the system, in Section 4.3 we aim to explain the occurrence of the chaos that is clearly visible in our bifurcation diagrams. We use one more notion of chaos, namely metric chaos, saying roughly speaking, that the map is chaotic when it admits an invariant measure, absolutely continuous with respect to the Lebesgue measure. This is one of the strongest notions of chaos.

4.1.
A few additional useful properties of the adaptation map. We start with a result that, for a range of v R values, establishes the non-monotonicity in the dynamical core of the adaptation map formed by the initial iterates of w * v R .
Remark. Note that from the above proof it follows that any given condition of the form is guaranteed to hold for vR sufficiently large to satisfy (4.1), for a large enough choice of l and ε sufficiently small.

Remark.
We introduced here the extra technical assumption that Φ ε (w) < −1 on (α, ξ); our numerical simulations suggest that this is always satisfied and that Φ has only one inflection point, located in (α, ξ).

Definition 4.2.
We say that the critical point We recall that under the current assumptions, Φ is at least C 3 since it is given by the flow of (1.1) with F being at least C 3 . In fact, in the most common cases, such as the adaptive exponential model F (v) = e v − v or the quartic model F (v) = v 4 + 2av, we can expect Φ ∈ C ∞ . Since the results below (Theorems 4.3 and 4.4) do not depend on ε > 0 but the dependence of Φ on v R is important in the proofs, we use the notation Φ v R for the adaptation map in the remainder of this subsection.
Proof. We want to show that Φ v R (w) = 0 for any w = w * . For w < w * the statement follows immediately from Lemma 3.1 (also from equation (3.5)). Now take w > w * . By P v R (w) < w * < w denote the w-coordinate of the crossing of the trajectory (V (t; w, v R ), W (t; w, v R )) with the reset line v = v R (below the v-nullcline). We have two possibilities: (a) P v R (w) > bv R (crossing above the w-nullcline), or (b) P v R (w) ≤ bv R (crossing below the w-nullcline).

In both cases we compute
where Φ v R (P v R (w)) = 0 by the previous (w < w * ) result. It remains to show that P v R (w) = 0.
First, consider case (a). We notice that the trajectory (V (t; w, v R ), W (t; w, v R )) between the point (v R , w) and (v R , P v R (w)) does not cross w-nullcline and thus can be seen as with the initial condition V (w) = v R . Hence we get the following implicit equation for P v R (w): We define a function where V (W ; w, v R ) is the solution of (4.4) with the initial condition V (w) = v R , such that (4.5) is equivalent to G(P v R , w) = 0. Since the point (v R , P v R (w)) lies apart from both the nullclines we compute and by Implicit Function Theorem we obtain that the mappingw → P v R (w) is a C 1 -function in the neighbourhood ofw = w. Consequently, P v R (w) exists. Moreover, by differentiating the equation (4.5) with respect to w we obtain that P v R (w) satisfies where we abuse notation by letting ∂V ∂w denote the partial derivative of V specifically with respect to the initial condition w, namely the second argument of V .
As the point (v R , P v R (w)) lies apart both the nullclines, it follows that Hence, based on equation (4.6), it suffices to show that ∂V ∂w (P v R (w); w, v R ) = 0. By an application of Peano's Theorem (see equation (3.4) of [21]), with (see Corollary 3.1 of [21]) where V = V (W ; w, v R ). Therefore, ∂V ∂w (P v R (w); w, v R ) = 0, as desired. The proof for case (a) is completed.
In case (b) one needs to notice that the trajectory (V (t; w, v R ), W (t; w, v R )) between the point (v R , w) and (v R , P v R (w)) crosses first the v-nullcline and then the w-nullcline. Therefore we cannot argue exactly as in case (a). Nevertheless, (b) can be reduced to (a) in the following way: There exists a neighbourhood U of (v R , w) and the reset valuev R < v R such that for every (v R ,w) ∈ U the trajectory (V (t;w, v R ), W (t;w, v R )) crosses the line v =v R exactly two times, say at points v R is such as in case (a): the crossing occurs between the two nullclines, not below both).
)) = 0. By expressing the trajectories (V, W ) as the function W = W (V ), we similarly obtain Theorem 4.4. Given ε > 0, if v R is sufficiently large such that F (v R ) > ε, then the point w * is non-degenerate; that is, Φ v R (w * ) < 0.
Proof. As in equation (3.5), for any v > v R , each trajectory W (v; v R , w) with initial condition (v R , w), w < w * , satisfies We already know that Φ (w) < 0 for w < w * . Using the flow-box theorem, it is also obvious that Φ (w) < 0 for w > w * close enough to w * . We aim at proving that Φ (w * ) < 0 by considering the limit of the difference quotient of Φ on the left of w * . Hence, since we already know that ∂W ∂w (v; v R , w * ) = 0 we want to prove that By assumption, F (v R ) > ε; that is, the reset line is bounded away from the knee of the v-nullcline. We introduce a parameter δ > 0 small but fixed and, for studying the w-limit, we only consider the values w ∈ [w min , w * ] for some w min < w * and Hence, any trajectory starting from (v R , w) with w ∈ [w min , w * ] remains above the w-nullcline at least for v ∈ [v R , v R + δ] and W (v; v R , w) decreases over this interval of v values. Then, for v > v R + δ, we split the integral and exponentials as follows: First note that the second factor is well-defined even around w = w * and v tending to +∞ since Hence, this factor converges uniformly in w (in the vicinity of w * ) and v (on On the other hand, the first term does not depend on v and its limit as w → (w * ) − is well-defined. Now, we focus on proving that this latter limit is strictly negative. We have, for w ∈ [w min , w * ), ds. (4.8) The first factor in the integrand is strictly positive for (s, w) ∈ [v R , v R +δ]×[w min , w * ]. Moreover, from (1.1) we have Using this expression, one obtains Finally, it follows that for any w ∈ [w min , w * ), And since, w − w * < 0, such that the limit for w → (w * ) − is strictly negative, which completes the proof.
Remark. Note that Theorems 4.3 and 4.4 do not require ε small, therefore they apply to the adaptation map in general, not only near the singularly perturbed limit.

Topological chaos.
Probably the most common definition of chaos is the one due to Devaney [17], which states that a continuous map f : X → X on a compact metric space X is chaotic if there exists a compact invariant subset Y ⊂ X (called a D-chaotic set ) such that f | Y is transitive, the set of periodic points of f | Y is dense in Y and f | Y has a sensitive dependence on initial conditions 2 (where f | Y denotes the restriction of f to the set Y ). A map satisfying these properties is called Devaney-or D-chaotic. We refer the reader e.g. to [1,3,20,64] for definitions of transitivity and sensitive dependence on initial conditions as well as results on redundancy of the last one or even sometimes the last two conditions in the definition of D-chaos.
The other notions of topological chaos include, among others, positive topological entropy, chaos in the sense of Block and Coppel (see e.g. [1,5]), Li-Yorke chaos (weaker than D-chaos and Block-Coppel chaos, see [1]), or even weaker notions such as the sensitive dependence on initial conditions itself or the existence of a period-3 orbit. It is not our aim to discuss here all these notions but just to make the observation that for the adaptation map all of them are very likely to occur.
To express the complexity of the dynamics, it is also useful to look for "horseshoes", defined for a one-dimensional map as follows.
A continuous map f : I → I (where I ⊂ R is a compact interval) has a 2-horseshoe (or in other words, f is turbulent), if there exist two closed-subintervals of I, A 1 and A 2 , with disjoint interiors, such that While Φ itself cannot feature a horseshoe, its iterates Φ k for k ≥ 2 may, resulting in chaotic dynamics of the map.
To carry out the proof of the forthcoming Theorem 4.7 that establishes the chaotic nature of the adaptation map, we need to introduce one more definition: Definition 4.6 (see e.g. [5]). A trajectory {Φ n (w)} n≥0 of some point w ∈ R will be called alternating if Φ k (w) < Φ j (w) for all even integers k and all odd j, or Φ k (w) > Φ j (w) for all even integers k and all odd j.
Statements 3 and 4 are general implications of 1 and 2. Indeed, statement 3 follows from 1 since the existence of a periodic point whose period is not a power of two is in this case equivalent to positive topological entropy (see Theorem 3.22 in [10] and Corollary 3 in [48]). Statement 4 can be derived as a general consequence of 2 based on Proposition 3.3, Theorem 4.1 and Theorem 4.2 in [1]. Remark. Theorem 4.7 could read "topological chaos occurs almost all the time", since the condition Φ 2 (w * ) < Φ 3 (w * ) < w * < Φ(w * ) is satisfied, for example, for any choice of vR such that F (vR) + I > p0 + 2d and any ε small enough (see Lemma 4.1 and Remark following it). Note also that a condition for existence of periodic orbits of all periods was given in [72,Theorem 3.4]. The assumption of Theorem 4.7 is weaker than the previous result and thus covers more cases.
While the above results prove the existence of an infinite 3 set Y on which the map is chaotic in the sense of the above Theorem 4.7, they do not ensure that the chaotic behavior is generic for arbitrary trajectories. In particular, the set Y may be small, even of zero Lebesgue measure. Similarly, the existing period-3 orbit might be stable and attracting for almost all initial conditions. The following section focuses on a stronger notion of chaos, namely metric chaos.

Metric chaos.
In this section we discuss the theoretical justification for the emergence of visible (metric) chaos occurring in our bifurcation diagrams (see for example the ε = 0.4 panel of Figure 3.2). For this purpose, we come back to the theory of unimodal maps (see e.g. [66]) and show that the standard theory may, under technical conditions on the adaptation map, directly apply to our system and justify the existence of metric chaos at the transitions. We recall the following definitions: Definition 4.8 (Metric chaos). We say that Φ is chaotic if it admits an absolutely continuous invariant probability measure ( acip) µ, i.e. an invariant measure that is finite, normalized and has density with respect to the Lebesgue measure.
Definition 4.9. A map Φ is called a Misiurewicz map if it has no periodic attractors and if critical orbits (i.e. forward orbit of the critical points) do not accumulate on critical points, that is, if where C denotes the set of critical points of Φ and ω(C) is its ω-limit set.
For some bounded interval of parameter values v R ∈ [v R1 , v R2 ] =: V consider the one-parameter family of adaptation maps {Φ v R } v R ∈V . We assume, as previously, . Under this condition it is not hard to construct an associated family of unimodal mapsΦ v R whose orbits are in a one-to-one correspondence with those of Φ v R (at least after a transient period). To this end, we start with a change of coordinates that makes the singular point w * independent of parameter v R , and define: These maps have their critical point atw * = 0 for all v R . Let J = [A, B] be a closed interval strictly containing all dynamical cores [Φ 2 v R (w * ),Φ v R (w * )] of the familyΦ v R . We can define a family of unimodal mapsΦ v R on J that are equal to corresponding . After a finite number of iterates, the orbits ofΦ v R are identical to those ofΦ v R , which are simple translations of those of Φ v R . Therefore, the asymptotic dynamics ofΦ v R and the original map Φ v R are the same. By studying the unimodal mapsΦ v R , we can bring the well-developed theory of unimodal maps [14,66] to bear to characterize the non-transient properties of the orbits of Φ v R .
It is clear that one can perform the above construction of the family {Φ v R } such that the following standard conditions are satisfied: [I ]Φ v R is a one-parameter family of C 3 unimodal maps of an interval J.
[II ] EachΦ v R has a unique and nondegenerate critical point w * (independent of v R ).
[III ] EachΦ v R has a repelling fixed point on the boundary of J.
Such maps display metric chaos as soon as: Unfortunately, it is complex to establish condition [V]. The main difficulty in verifying [V] is in assuring that there are no periodic attractors for the map Φv R . Showing the absence of attractive periodic orbits is a complex task for our general model, and even simpler sufficient conditions such as those relying on negative Schwarzian derivatives atv R present deep difficulties to demonstrate. Another issue is that the accumulation of the orbit of w * in a neighbourhood of w * must not occur; this can be avoided in particular when the critical point w * is mapped in a few iterations onto the unstable fixed point w f v R (see [66,Section 6.2]). This condition is actually satisfied for some intermediate parameter valuev R in the period-incrementing transition. Indeed, it is relatively easy to show, using the intermediate value theorem, that: Proposition 4.10. For sufficiently small ε, the family {Φ v R } v R ∈V of adaptation maps undergoes period-incrementing transitions such that between any two intervals J k = [a k , b k ] and J k+1 = [a k+1 , b k+1 ] of v R values, corresponding, respectively, to admissible k and k + 1 attracting periodic orbits, there exists a parameter valuev R ∈ (b k , a k+1 ) such that the critical point is mapped into a few steps onto the fixed point. The correspondence of the above result is obviously also satisfied for the family {Φ v R }. Moreover, sufficient conditions for the fixed point w f v R to be unstable are provided in Lemma 4.1. For completeness, one would also require some additional non-degeneracy condition on how the unstable fixed point and the point that is eventually mapped onto it evolve with the change of the parameter v R in the neighbourhood ofv R . This condition is technical and therefore for its precise statement we refer to [66] (see condition H6 therein). It is difficult to demonstrate rigorously at full generality, yet it is completely generic.
Our analysis and extensive numerical simulations suggest that these technical conditions are generally met. When they hold, the theory of unimodal maps [66,Theorem 18 and Corollary 19] yields the conclusion that there exist constants γ > 0 and C > 0 and a positive measure set E ⊂ V withv R ∈ E as a Lebesgue density point, such that ∀v R ∈ E, ∀n ∈ N * , |(Φ n v R ) (Φ v R (w * ))| ≥ Ce γn . Furthermore, under the above conditions and if the map has a negative Schwarzian derivative in a neighborhoodV ⊂ V of v R , then E ⊂V can be chosen so that for all v R ∈ E Φ v R exhibits metric chaos with an acip µ v R , describing the asymptotic behavior of almost all orbits, and that has a positive Lyapunov exponent almost everywhere: lim n→∞ 1 n log |(Φ n v R ) (w)| = κ > 0 for a.a. w ∈ R (4.13) We illustrate the application of this theory in Figure 4.2 within the range of values associated with the transition from bursts with 4 to bursts with 5 spikes. We evaluated the critical valuev R numerically and simulated the orbits of the system for distinct initial conditions. These orbits collapse on the unique measure depicted in the figure.  5. Discussion. The study of the bursting patterns in nonlinear adaptive integrateand-fire neuron models led us to identify and elucidate the underlying period-incrementing structure that organizes spike patterns in a regular progression as parameters (particularly, the reset value of the membrane potential) are varied. The structure exhibited here is general and does not depend on the specific adaptive integrate-and-fire model considered. We have shown that this structure is present in particular when the adaptation variable is much slower than the voltage variable. Indeed, we have established that in the limit of perfect separation of timescales, the sequence of adaptation values approach orbits of a discontinuous piecewise linear discrete dynamical system having, for any set of parameters, a unique globally attractive periodic orbit, with a period that is incremented instantaneously as parameters are varied. In the full system, we have shown that while the period-incrementing structure was globally conserved, transitions are no longer instantaneous, and the periodic orbits bifurcate and lose stability. We have investigated in more detail the presence of chaos during these transitions and shown that within these regions, the system exhibits topological chaos and may display metric chaos as well, under the additional assumption that the adaptation map has a negative Schwarzian derivative, at least for a specific value of the voltage reset.
From the mathematical viewpoint, the relative simplicity of the model allowed us to make significant progress in characterizing the complex period-incrementing transition structure. In the class of models investigated, however, the fact that the adaptation map is not known in closed form raises some difficulties and in particular imposes limitations to the results shown. With the same methodology, one may thus achieve stronger results for a specific choice of model in which the adaptation map satisfies a few additional properties. In particular, establishing that a particular map features a negative Schwarzian derivative away from the critical point would allow an extension of the results on metric chaos developed in Section 4.3 as well as a more complete characterization of the structures of the topological and metric attractors and their inter-relationships (see [46,6,8] as well as the survey article [66] and references therein). For instance, considering our result on the non-degeneracy of w * , we know that the Schwarzian derivative is negative in a neighborhood of the critical point, which allows the application of van Strien's theory [73] showing the uniform boundedness of the periods of periodic attractors and non-hyperbolic periodic orbits of Φv R . Using the same theory, if all periodic points of Φv R are hyperbolic and repelling, then Φv R admits an acip of positive entropy.
A number of further questions arise at this stage. For instance, since the adaptation map Φ is continuous and unimodal, a natural question is to investigate to what extent the dynamics of Φ may be similar to that of the canonical logistic map F µ : x ∈ [0, 1] → µx(1 − x) with µ ∈ (0, 4]. While Milnor-Thurston kneading theory [47,14] ensures that our unimodal extension of the adaptation map is semiconjugated to a map {F µ } through a continuous surjective and non-decreasing map, the conclusions that may be drawn remain limited. In particular, one would need again to show the uniqueness of periodic attractors to ensure that the semi-conjugacy is actually a conjugacy [13], but typically this is not the case and the dynamics of Φ is not completely equivalent to that of the logistic map. From the application viewpoint, these results yield a deeper understanding of the dynamics and its parameter-dependence for a widely used class of hybrid neuronal models, in the case where the subthreshold dynamics has no fixed point. This configuration corresponds to settings where the input to the neuron is sufficiently large. Interestingly, one scenario where neurons receive unusually high input is during seizures within the epileptic brain (e.g. [45] and references therein). Chaotic dynamics has been associated with epileptic brain dynamics [62]; it is also possible that seizures represent transitions from chaotic to more regular dynamics within high input regimes [22]. Hence, our work has possible relevance for the use of bidimensional hybrid neuronal models to study epileptic dynamics. Our analysis is based on an adaptation map, which can be defined on the whole real line in this situation. In the companion paper [57], we will switch gears and investigate the case where the subthreshold system has two unstable fixed points, a spiral and a saddle, with a heteroclinic orbit from the former to the latter. In that setting, we obtain and study the corresponding, distinctive forms of dynamics that arise, which feature alternations of small oscillations and spikes or bursts, and are known as mixed-mode oscillations. Again, our analysis will rely heavily on the relatively simple geometric structure of the hybrid system.