Wild oscillations in a nonlinear neuron model with resets: (II) Mixed-mode oscillations

This work continues the analysis of complex dynamics in a class of bidimensional nonlinear hybrid dynamical systems with resets modeling neuronal voltage dynamics with adaptation and spike emission. We show that these models can generically display a form of mixed-mode oscillations (MMOs), which are trajectories featuring an alternation of small oscillations with spikes or bursts (multiple consecutive spikes). The mechanism by which these are generated relies fundamentally on the hybrid structure of the flow: invariant manifolds of the continuous dynamics govern small oscillations, while discrete resets govern the emission of spikes or bursts, contrasting with classical MMO mechanisms in ordinary differential equations involving more than three dimensions and generally relying on a timescale separation. The decomposition of mechanisms reveals the geometrical origin of MMOs, allowing a relatively simple classification of points on the reset manifold associated to specific numbers of small oscillations. We show that the MMO pattern can be described through the study of orbits of a discrete adaptation map, which is singular as it features discrete discontinuities with unbounded left- and right-derivatives. We study orbits of the map via rotation theory for discontinuous circle maps and elucidate in detail complex behaviors arising in the case where MMOs display at most one small oscillation between each consecutive pair of spikes.

1. Introduction. In this paper, we continue our study of hybrid integrate-andfire neuronal models from [49], turning our attention to the analysis of mixed-mode oscillations (MMOs). MMOs are trajectories exhibiting small or subthreshold oscillations alternating with one or more large amplitude oscillations or spikes. These appear in a variety of cell types and brain areas including inferior olive nucleus neurons [5,34,35], stellate cells of the entorhinal cortex [1,2,25,61], and neurons in the dorsal root ganglia [4,32,33], as well as in thalamocortical spindle waves [36]. Neurons transmit information through the timing of spikes and the pattern of spikes fired, and subthreshold oscillations and associated MMO patterns may contribute to the precision, timing, and robustness of neuronal spiking [29,61,45] as well as to spatial navigation [13]. The major goal of this work is to provide a detailed mathematical analysis of MMOs in a class of neuronal models. Specifically, in this article, we (i) consider a class of planar hybrid models widely used to model the electrical activity of neurons, (ii) show that models in this class are able to generate a wide range of MMO patterns, (iii) introduce a general mathematical framework for studying the dynamical structure involved and the orbits that result, and (iv) describe the geometric mechanism underlying these patterns.
From the biological viewpoint, neuronal activity patterns, including MMOs, rely on ionic and biochemical mechanisms that are accurately described by nonlinear dynamical systems of relatively high complexity, such as variants on the celebrated Hodgkin-Huxley model [17,50]. As described in our companion paper [49], in contrast with detailed biophysical models, integrate-and-fire models are abstractions of the voltage dynamics in which differential equations describing the dynamics of membrane depolarization are combined with a discrete reset corresponding to the emission of an action potential (a spike) and subsequent hyperpolarization. These models, first introduced more than a century ago [30], have evolved to incorporate nonlinearities to model the fast dynamics of spike initiation [8,10] and additional variables modeling adaptation [19], synaptic dynamics [38] or resonant properties [22]. Among these models, nonlinear bidimensional integrate-and-fire models with blow-up and resets are widely used in computational neuroscience, owing to their relative simplicity yet very rich dynamical phenomenology [7,19,21,57,60]. However, none of these studies reported the presence of MMOs in these systems.
More generally and despite their importance for applications, MMOs have so far received little attention in hybrid systems. A notable exception is the work of Rotstein and collaborators on linear bidimensional resonate-and-fire neuron models [22]. These models are organized around an unstable focus and naturally exhibit a variety of MMOs. Numerical simulations guided by characterizations of the trajectories and timescale analysis were used to explore associated subthreshold dynamics [47,48] and to show how they can give an abrupt increase in firing frequency [44]. MMOs in this class of models result from a combination of subthreshold oscillations together with subsequent threshold crossings corresponding to spikes. The multi-timescale MMO scenario in these models does not necessarily represent the general mechanism for MMOs in hybrid models, however, and to date, there has not been a thorough analytical investigation of the detailed mechanisms underlying MMOs, incorporating both subthreshold and spiking components, in these models in the absence of timescale separation.
In the current manuscript we present a rigorous study of MMOs in nonlinear bidimensional integrate-and-fire neuron models. We show that these can exhibit a wide variety of MMO patterns when the subthreshold dynamics features two unstable fixed points, a saddle and an unstable focus. We investigate spike patterns through iterates of a discrete map, the so-called adaptation map introduced in [60] (see the companion paper [49] for more details on the construction and use of this map). While previous works have considered settings in which the adaptation map is continuous, here, in the presence of an unstable focus, we will show that the adaptation map is singular: it may be undefined on a countable set of values at which the map has well-defined and finite left and right limits and infinite one-sided derivatives. A number of difficulties emerge from the irregular nature of the map; since associated circle maps may also feature analogous singularities, classical theories of Poincaré and Denjoy of circle homeomorphisms or their extensions to continuous non-invertible maps ( [39]) do not apply, and because of the unbounded derivative, neither do theories of discontinuous contractive maps [14,27]. This contrasts with previous detailed studies of interspike intervals for periodically driven one-dimensional integrate-andfire models [9,12,28,37,53,56]. Here, we will demonstrate a fundamental relationship between the type of MMO pattern arising and the rotation number of the adaptation map. With the aim of characterizing rotation numbers of these maps, we build upon a number of theoretical results on circle maps that may have discontinuities [6,14,27,39,42,43] and sometimes extend these to singular maps with unbounded derivative. In this way, we describe a new mechanism underlying robust MMOs, not requiring multiple timescales, in hybrid dynamical systems constituting an important class of neuron models.
Our presentation of these results is organized as follows. In section 2, we introduce the model studied, review a few results on its dynamics, and describe the geometric mechanisms underlying the generation of MMOs. We detail the properties of the adaptation map in section 3, with a particular focus on discontinuity points and divergence of the derivative, which is proved to be a general result based on a Poincaré section encompassing the stable manifold of a saddle. We further show that the particular structure of the map ensures that any type of transient MMO can be generated by these neuron models. In section 4, we use discontinuous rotation theory to develop a precise description of the dynamics in the case where the adaptation map admits one discontinuity in its invariant interval. Implications and perspectives in dynamical systems and neuroscience, as well as some extensions and prospects for analyzes of cases with more discontinuities, are discussed in sections 5 and 6.
2. Hybrid neuron model and the geometry of the MMO mechanism.
In this work, we study the class of integrate-and-fire neuron models introduced in [57], described in detail in the companion paper [49]: where ε, b > 0 and I are real parameters. Following [57,60], we will assume that the real function F is regular (at least three times continuously differentiable), strictly convex, superquadratic at infinity, with its derivative having a negative limit at −∞ and an infinite limit at +∞ (see Assumption (A1) in the companion paper [49]). These assumptions imply in particular that the membrane potential blows up in finite time and at this explosion time, say t * , the adaptation variable converges to a finite value w(t * − ) [58]. At time t * , it is considered that the neuron has fired a spike; the voltage is instantaneously reset to the fixed reset value v R and the adaptation variable is updated as follows: with γ ≤ 1 and d ≥ 0; in this work, we allow γ < 1, which arises in accounting for spike duration (see [49] for details). Numerical simulations performed in the present article correspond to the case of the quartic model F (v) = v 4 + av with, unless otherwise stated, The values of the parameters associated with resets, d and γ, are left free and will be used as bifurcation parameters.
As discussed in the companion paper [49], the 1-dimensional adaptation map Φ can be defined based on the orbits of the system in the phase plane. Specifically, if is the associated value of the adaptation variable after spike and reset. In [49,60], we detailed the mathematical analysis of (2.1)-(2.2) in the absence of fixed points (yellow region of Fig. 1.1 of [49]). In that regime, the system fires an action potential for any initial condition in the phase plane. The adaptation map is thus well-defined on the whole real line and smooth. The orbits of the map can be used to analyze spike patterns and transitions between them.
In the present manuscript, we analyze the dynamics in regions in which the system features an unstable spiral point and a saddle (pink region of Fig. 1.1 of [49]). The stable manifold of the saddle is a one-dimensional heteroclinic orbit spiraling out from the unstable focus (see Fig. 1.1 of [49] and Fig. 2.1). This geometry of the phase space constrains trajectories reset within the spiral to proceed to a prescribed number of rotations around the unstable fixed point before firing ( Fig. 2.1, inset and bottom). The rotations around the unstable fixed point provide small oscillations used to define mixed-mode oscillations as follows.
Definition 2.1. Mixed-mode oscillations (MMOs) for the system (2.1)-(2.2) are spiking orbits consisting of an alternation of small oscillations and spikes. MMO patterns formed by a sequence of L k ∈ N * := {1, 2, ...} spikes followed by s k ∈ N * small oscillations are characterized by their signature L s1 1 L s2 2 L s3 3 · · · . Periodic signatures with period k are only denoted by finite sequence of length k, L s1

Remark.
• MMOs featuring bursts of two or more consecutive spikes not separated by periods of small oscillations (i.e., L k ≥ 2 for some k) are referred to as mixed-mode bursting oscillations (MMBOs). We use the term MMO as a generic term to describe any combination of spikes and small oscillations, and the term MMBO is applied specifically to distinguish those trajectories featuring bursts and small oscillations.
• In the present paper, we will be able to distinguish small oscillations at half-rotation precision, and thus will extend the definition above to signatures with half-integer number of small oscillations s k ∈ 1 2 N * .
We henceforth assume that the reset line {v = v R } intersects the spiraling stable manifold of the saddle, as in Fig. 2.1. The adaptation map is undefined at each intersection of the reset line with the stable manifold of the saddle, since the orbit of (2.1) starting from such a point converges to the saddle, and thus no spike follows. For any initial condition (v R , w) not on the stable manifold, the associated orbit performs a specific number of small oscillations before firing, resulting in an MMO pattern. As indicated in Fig. 2.1, the present framework allows us to perform a detailed analysis of this scenario, since: • the fact that the stable manifold is bounded in the v variable implies that the amplitudes of small oscillations, similarly to biological MMOs, are considerably smaller than the spike amplitude, and • the intersections of the stable manifold with the reset line partition the values of w associated with a specific number of small oscillations (with half-rotation precision). The signature of the MMO patterns can be deduced from a dynamical analysis of the adaptation map. The main objective of the manuscript is to characterize these patterns, and the main results are summarized below.  Phase plane with v and w nullclines (dashed black) and stable (red) and unstable (blue) manifolds of the saddle; the stable manifold winds around the repulsive singular point. The reset line {v = v R } (solid vertical line) intersects the stable manifold, separating out regions such that trajectories emanating from each undergo a specific number of small oscillations (colored segments, here from 0 to 3 below the w-nullcline and from 3.5 to 0.5 above). (Lower rows) The solution for one given initial condition in each segment. Note that the time interval varies in the different plots (indicated on the x-axis). Simulations had initial conditions v = v R = 0.012 and w chosen within the different intervals on the reset line.

Stable manifold
We establish that, as a transient behavior, the system can feature MMOs with all possible finite signatures (Proposition 3.3 and Corollary 3.4). Non-transient behaviors are deduced from the iterates of the adaptation map, which may feature several discontinuities and therefore support a very wide range of possible dynamics. We concentrate in section 4 on the case where the adaptation map features a single dis-continuity within its invariant interval. As in the seminal study of Keener on maps with one discontinuity [27], we distinguish two cases depending on the monotonicity of the lift, called overlapping or non-overlapping cases. In the non-overlapping case (see subsection 4.1), we characterize the rotation number of the associated adaptation map and show that it characterizes the MMO signature (Theorem 4.3) or the chaotic nature of the spike pattern fired. In the overlapping case (see subsection 4.2), the adaptation map yields rotation intervals with rational numbers corresponding to periodic orbits with MMOs (Proposition 4.6). To go beyond this description, we provide conditions that guarantee existence of periodic orbits with arbitrary periods, all displaying MMBOs (Proposition 4.8). Eventually, we discuss how the methods used here could be extended to cases with multiple discontinuity points within the invariant interval of the adaptation map (Theorem 5.1).
3. The adaptation map. We start by characterizing the properties of the adaptation map Φ given in (2.4). In the absence of singular points of the subthreshold dynamics, it is defined and continuous on R, and the nature of its orbits distinguishes regular spiking (fixed point of the map), bursting (periodic orbit of the map) or chaotic spiking [49,60]. In the present case with two singular points (unstable focus and saddle), we show that Φ is undefined at specific points, no longer continuous and has unbounded derivative, but its orbits still provide all the information necessary to characterize the associated MMO patterns.
3.1. Properties of the adaptation map. Throughout the manuscript, we assume that the vector field (2.1) has two unstable singular points, the repulsive singular point (v − , w − ) = (v − , F (v − ) + I) and the saddle singular point (v + , w + ) = (v + , F (v + ) + I), with v − < v + (see [57,60] providing detailed bifurcation analysis of the subthreshold dynamics). We denote by W s and W u the stable and unstable manifolds of the latter singular point; each of these is decomposed into two branches ( Fig. 3.1) with W s − (W s + ) extending towards w < w + (w > w + ) and W u − (W u + ) extending towards v < v + (v > v + ). The shape of the map Φ is organized around a few important points (see Fig. 3.1): • We denote by w * = F (v R ) + I the intersection of the reset line v = v R with the v-nullcline. • We denote by w * * = bv R the intersection of the reset line with the w-nullcline. • We denote by {w i } p i=1 the sequence of intersections of the reset line with W s , labeled in increasing order with respect to the value of w. As long as v R = v − , there exists a finite number of such points or none depending on the parameter values: an even number of intersections for v R < v − and an odd number for v R > v − . We denote by p 1 the index such that {w i } i≤p1 are below the v-nullcline and {w i } i>p1 are above; it is easy to see that p 1 is the smallest integer larger than or equal to p/2, i.e. p 1 = p/2 . The points {w i } split the real line into p + 1 intervals that we denote {I i } p i=0 , with I 0 = (−∞, w 1 ), I i = (w i , w i+1 ), i = 1, 2, ..., p − 1, and I p = (w p , ∞). Remark that these intervals precisely correspond to those in which the number of small oscillations occurring between two consecutive spikes is constant, except the interval I p1 , which is split into two subintervals by w * (see Fig. 2.1). The number of small oscillations for trajectories starting from I i is if i = p 1 and w < w * , p 1 + 1/2 if i = p 1 and w > w * and p is even if i = p 1 and w > w * and p is odd. (3.1) • We denote by w − lim < w + lim < ∞ the limit of the adaptation variable when v → +∞ along W u − and W u + respectively. In addition, we introduce the corresponding values obtained through the reset mechanism:  With these points defined, we can characterize the shape of the adaptation map. When we refer to the adaptation map, we abuse notation and use w i to denote the w-coordinates of the points of intersection of W s with {v = v R }.
Theorem 3.1. The adaptation map Φ has the following properties.
In any given interval I i with i ∈ {0, · · · , p}, the map is increasing for w < w * and decreasing for w > w * . 4. At the boundaries of the definition domain D, {w i } p i=1 , the map has welldefined and distinct left and right limits: 5. The derivative Φ (w) diverges at the discontinuity points 1 : 6. Φ has a horizontal plateau for w → +∞ provided that Moreover, for any w taken between the two branches of the unstable manifold of the saddle, hence in particular for w ∈ (w 1 , w p ), Φ(w) > β. In comparison to the case without singular points [60, Theorem 3.1], the map loses continuity, convexity, and uniqueness of the fixed point, but the monotonicity property (point 3), the presence of a plateau (point 6) and the comparison with identity (point 7) remain true. The presence of discontinuities and divergence of the map derivative substantially change the nature of the dynamics as we will see below.
It is worth noting that this divergence is a general property of correspondence maps in the vicinity of saddles (see Fig. 3.2), as we show in the following: Lemma 3.2. Consider a two-dimensional smooth vector field (at least C 2 ) with a hyperbolic saddle x 0 associated with the eigenvalues −µ < 0 < ν of the linearized flow. We denote by W s and W u the stable and unstable manifolds of the saddle and consider two transverse sections S s and S u intersecting W s and W u at x s and x u , respectively. There exists Ω s a one-side neighborhood of x s on S s that maps onto a one-side neighborhood Ω u of x u on S u . The correspondence map Ψ : Ω s → Ω u is differentiable in Ω s and we denote by Ψ the one-sided differential of Ψ at x s . We have: When µ = ν, the differential is finite and its value depends on the precise location of the sections.
Proof. Let us start by considering the linearized system in the vicinity of the saddle singular point. In the basis that diagonalizes the Jacobian, we can write the system in the simple form: ẋ = −µẋ y = νy and considering a section S s corresponding to y = y 0 and a section S u corresponding to x = x 0 > 0, simple calculus leads to the formula that the correspondence map ϕ of the linearized system between S s and S u is defined for x ≥ 0 by ϕ(x) = (y 0 /x ξ 0 )x ξ with ξ = µ/ν. Hence, the differential of ϕ at 0 + is such that:

Linearized system
Smooth conjugacy Typical topology of manifolds and sections in Lemma 3.2: we consider correspondence map between sections Ss (red) and Su (orange) transverse, respectively, to the stable and unstable manifolds (black lines) of the saddle (orange circle). Typical trajectories are plotted in blue.
The key arguments are the characterization of correspondence maps associated with the linearized system (upper left inset) between two transverse sections S s and S u , and the smooth conjugacy between the nonlinear flow and its linearization.
• it diverges if ξ < 1, hence for ν − µ > 0 (i.e. if the dilation along the unstable direction is stronger than the contraction along the stable direction); • it vanishes if ξ > 1, hence for ν − µ < 0 (i.e. if the contraction along the stable direction is stronger than the dilation along the unstable direction); • when ξ = 1 (i.e. contraction and dilatation are of the same intensity), we find ϕ (0 + ) = y 0 /x 0 which depends on the precise location of the sections. To demonstrate the lemma, we thus need to show that the same result holds for the nonlinear system. The Hartman-Grobman Theorem [15], which ensures that the nonlinear system is conjugated to its linearization through a homeomorphism in the vicinity of the (hyperbolic) saddle, will not be sufficient; we need to ensure that the nonlinear and linear flows are locally conjugated via smooth diffeomorphisms (at least C 1 ). Finding smooth conjugacy is a subtle question for a general dynamical system that has been the object of significant research and generally requires avoiding resonances in the eigenvalues, which may lead to a relatively complex relationship [51,54]. In two dimensions, the problem is simpler and it was proved in [55] that any C 2 planar dynamical system in the neighbourhood of a saddle is smoothly (with at least C 1 regularity) conjugated with its linearization, and the derivative of this conjugacy is bounded away from 0 in a sufficiently small neighborhood of the saddle (since this conjugacy converges in a C 1 -sense towards the identity close from the saddle). Completing the proof thus only amounts to showing that the correspondence maps from a neighborhood of S s to S s and from S u to a subset of S u are smooth with differential bounded away from zero and infinity. This is a classical consequence of the flow box theorem and regularity with respect to the initial condition.
Now that this general result is proved, we proceed to establish the properties of the adaptation map by proving Theorem 3.1.
Proof. of Theorem 3.1: First, note that generalization of the reset mechanism by introducing a parameter γ ∈ (0, 1] does not substantially impact the shape of the adaptation map. Indeed, all properties rely on the map associating with a point on the reset line (v R , w) the value ϕ(w) := W (t * − ; v R , w) of the adaptation variable at the time t * of the subsequent spike, since ϕ(w) = (Φ(w) − d)/γ. In other words, the generalization of the reset mechanism does not introduce any new mathematical difficulty. Hence, the proofs for items 1. to 3. and 6. to 8. are straightforward extensions of the analogous proof in [57] or simple algebra. Similarly, the proof of property 4. follows reasoning analogous to that of 3., with left and right limits found by finely characterizing the shape of the trajectories as w approaches one of the discontinuity points w i . In all cases, the trajectory will initially remain very close to the stable manifold, before leaving the vicinity of the stable manifold near the saddle and following the unstable manifold. Depending on whether the trajectory approaches the saddle from the right or from the left, it will either follow the left or right branch of the unstable manifold, hence either converge towards w + lim or w − lim . We focus on the proof of property 5., which requires specific analysis. We use Lemma 3.2 and prove that the conditions on the contraction and dilation near the saddle are satisfied. We consider the specific sections that define Φ, namely S s = {v = v R } (which is valid as long as the stable manifold is not tangent to the reset line) and a section corresponding to spiking, denoted with a slight abuse of notation as S u = {v = +∞}. The use of a section at ∞, however, does not exactly fit the statement of Lemma 3.2, and requires us to show that the differential of the correspondence map does not vanish as v → ∞.
First, note that since (v − , F (v − ) + I) is an unstable focus, the linearized flow there has two complex conjugate eigenvalues with positive real part and therefore the trace of the Jacobian, given by F (v − ) − ε, is strictly positive. Since F is convex, the trace of the Jacobian at the saddle equals Hence the dilation at the saddle is always stronger than the contraction; in the notation of (3.3), we have ν − µ > 0.
To show that the infinite derivative persists when one considers S u = {v = +∞}, we express the map Φ formally in the region below the stable manifold of the saddle, which all spiking trajectories cross. In this region, any trajectory has a monotonically increasing voltage (that blows up in finite time), and the orbit with initial condition (v R , w 0 ) can be expressed as the parametric curve (v, W (v)) with The expression of the differential of W with respect to w 0 at v is given by: with solution given, as a function of the trajectory W , by (see Peano's Theorem in [16]) 2 : Hence, for any section S u = {v = θ} (with θ < ∞), the derivative of the map Φ cannot vanish. Furthermore, for u large, we know that W (u) remains finite and thus the integrand in (3.6) behaves as −ε/F (u) which is integrable at infinity (cf. Assumption (A1), [49]). Consequently, the integral within the exponential term does not diverge towards −∞ as v → ∞ and the derivative (3.6) does not vanish at S u = {v = ∞}. We further note that all correspondence maps away from singularities (v − , F (v − )+I) and (v + , F (v + ) + I) are well-defined and with finite derivative bounded away from zero for the same reason. The intervals (−∞, w 1 ), {I i } p i=1 , and (w p , ∞) on the line {v = v R } are transverse sections of the flow and correspondence maps from I i to (−∞, w 1 ) are increasing for i ≤ p 1 (hence the left and right differentials of Φ at w ± i for i ≤ p 1 are equal to +∞) and decreasing otherwise (hence the left and right differentials at w ± i for i > p 1 are equal to −∞).

Transient MMO behaviors.
We recall that at each discontinuity point w i , the right and left limits of the adaptation map are always equal to either α or β. This property, related to the fact that all discontinuities correspond to intersections of the reset line with the stable manifold of the saddle, is a very important property that endows the system with a rich phenomenology, ensuring that it can generate MMOs of any signature.
We start by treating the case where the adaptation map has an infinite number of discontinuity points, which occurs in particular 3 when the reset line We denote by {m i } i∈N * the w values of the discontinuity points below the intersection w * of the reset line with the vnullcline, with m i < m i+1 for any i. Similarly, we denote by {M i } i∈N * the w values of the discontinuity points satisfying M i > w * and M i > M i+1 . We note that the left and right limits of Φ at m i (M i ) are well defined, equal to α and β, respectively (β and α, respectively) 4 . Proposition 3.3. Assume that the reset line has an infinite number of intersections with the stable manifold of the saddle. If moreover all the discontinuity points , then for every n ∈ N * and every finite sequence with non-empty interior such that for any w 0 ∈ J, the orbit with initial condition (v R , w 0 ) has a transient signature Proof. We recall that for w ∈ (m k , m k+1 ) (resp. w ∈ (M k , M k+1 )), the orbit passing through (v R , w) performs exactly k (resp. k+1/2) small oscillations before spiking. Thus, proving the proposition amounts to finding a set of initial conditions with a prescribed topological dynamics. In detail, given an MMO pattern 1 s1 1 s2 . . . 1 sn , where the s i are as above, we are searching for sequences of iterates of Φ falling sequentially in the intervals The set of initial conditions corresponding to this prescribed signature is therefore , and proving the theorem amounts to showing that this set is not empty, which relies on the particular shape of the map Φ and specifically on the fact that Φ(I s k ) = (β, α) for any admissible s k . This property implies that for any interval J ⊂ (β, α) with non-empty interior and any admissible s k , the intersection of the pre-image Φ −1 (J ) with I s k is an interval with non-empty interior. In turn, this fact allows us to establish the proposition by recursion on the length of the transient signature n. Indeed, for n = 1, the set of initial conditions associated to a transient signature 1 s1 is the set J 1 = Φ −1 (I s1 ), which has a nonempty interval of intersection with both (m k , m k+1 ) and (M k , M k+1 ). Let us now assume that the same property is true for some n ∈ N * , namely that for any sequence {s i } i=1···n , there exists a set of initial conditions with a non-empty interval of intersection with both (m k , m k+1 ) and (M k , M k+1 ) from which trajectories have the transient signature 1 s1 · · · 1 sn . Let us now fix a sequence {s i } i=1···n+1 . By the recursion assumption, the set J n associated to the transient signature 1 s2 1 s3 ...1 sn+1 is such that J n ∩ I s1 is an interval with non-empty interior. Consequently, the set J n+1 = Φ −1 (J n ∩ I s1 ) contains a non-empty interval of intersection with all I s k and any trajectory with initial condition within that interval has the transient signature
If the reset map has a finite number of discontinuities, exactly the same proof applies to show that any MMO pattern with an accessible number of small oscillations exists. This extension is precisely summarized in the following result.
Corollary 3.4. Suppose that the reset line v = v R has a finite number p of intersections with W s , with w-coordinates ordered as w 1 < w 2 < ... < w l < w l+1 < ... < w l+k < .. < w p , where l ∈ N * denotes the largest index i such that w i < β and exactly k ≥ 2 intersections lie in [β, α]: β ≤ w l+1 , ..., w l+k ≤ α. Let S denote the set of numbers of small oscillations performed by the trajectories with initial condition in I l+1 = (w l+1 , w l+2 ), I l+2 = (w l+2 , w l+3 ), ..., I l+k−1 = (w l+k−1 , w l+k ) according to the formula (3.1). Then for any n ∈ N * and any sequence If the number of intersections w i of the reset line v = v R with W s is infinite but only k of them lie in the interval [β, α], then we have exactly two possibilities: • all the points w i in [β, α] are not greater than w * and for every N ∈ N * and every sequence {s i } N i=1 with s i ∈ {l + 1, l + 2, ..., l + k − 1} there exists an interval J ⊂ [β, α] such that every initial condition w ∈ J yields an MMO with the pattern 1 s1 1 s2 ...1 s N , where the index l is obtained from the ordering • all the points w i in [β, α] are greater than w * and for every N ∈ N * and every sequence {s i } N i=1 with s i ∈ {l + 3/2, l + 5/2, ..., l + (k − 1) + 1/2} there exists an interval J ⊂ [β, α] of initial conditions yielding MMOs with the pattern 1 s1 1 s2 ...1 s N , where the index l is obtained from the ordering This corollary covers all cases studied in this paper, including finite or infinite number of intersections of the stable manifold with the reset line. Only the number of these points in the interval [β, α] determines the possible MMO patterns.
4. Adaptation map with one discontinuity point in the invariant interval. The general description developed above does not yield a precise specification of the dynamics of the system. For clarity of exposition, from now on, we shall assume that v R < v − and we focus chiefly on the case where the adaptation map has exactly one discontinuity point w 1 in the interval [β, α] (although there may be arbitrarily many outside of that interval). One of the main limitations of this situation is that the resulting MMOs have at most one small oscillation between spikes. Nonetheless, this case is advantageous in that the number of possible configurations of the map and identity line remains relatively limited, while there is a combinatorial explosion in cases with more intersections. It will be clear that most of our techniques extend beyond these situations under suitable technical assumptions 5 .
The shape of the map depends on the relation of certain points, as represented in Fig. 4.1, and we list several relevant conditions that we will consider as we proceed: (C1) There exists a unique discontinuity point w 1 in the interval [β, α]: is invariant, which is set by the conditions: Non-transient regimes only depend on the properties of the map Φ in a bounded invariant interval. Indeed, we have seen in Theorem 3.1 that Φ is bounded above and that for w small enough Φ(w) > w, implying the existence of an invariant compact set I in which any trajectory is trapped after a finite number of iterations. This remark opens the way to consider Φ as a circle map (after identifying the endpoints of I) and thus to use rotation theory in order to rigorously discriminate (i) whether the firing is regular, bursting or chaotic, corresponding respectively to fixed points, periodic orbits, or chaotic (non-periodic) orbits of Φ (see [60]), as well as (ii) the number of  the non-overlapping case (assumption (C4)) is decomposed into regions A (yellow), B (pink) and C (orange) according to the position of the discontinuity point w 1 with respect to Φ(α) and Φ(β). The overlapping case (assumption (C4')) is decomposed into regions D (assumption (C2)) and E (assumption (C2')) according to the position of the critical point w * . Note that, for the given value of v R , in the range of (d, γ) values, assumption (C3) is always satisfied. For each subregion, a prototypical scheme of adaptation map Φ is displayed. The grey regions in the map plots D and E highlight the overlap. small oscillations occurring before a spike, according to the partition of Fig. 2.1, i.  and extended on R through the relationship that for any x ∈ R and k ∈ Z, The rotation number of Ψ at w ∈ R is defined as: provided that the limit exists.
An example of the lift is given in Fig. 4.5. Note that the lift is continuous on the interior of the invariant interval I. Indeed, it is continuous on (β, w 1 ) and (w 1 , α) since Φ is continuous therein, and at x = w 1 , both its left limit Ψ(w − 1 ) = Φ(w − 1 ) and right limit Ψ(w + 1 ) = Φ(w + 1 ) + (α − β) are equal to α. However, the map Ψ is generally not continuous on R and displays discontinuities at the points x k = α + k(α − β) for k ∈ Z when Φ(β) = Φ(α) (which is generally the case). By convention, the above definition introduces Ψ as a left-continuous map. As will be emphasized at relevant places, this choice does not impact our developments, and in particular does not affect possible values of rotation numbers. We finally note that the maps Ψ and Φ restricted to [β, α] induce the same circle map ϕ : S |I| → S |I| on the circle of length |I| = α − β, and the orbits of Ψ coincide modulo |I| with the orbits of Φ, except at w 1 where the map Φ is not defined 6 . Therefore Ψ captures well the general dynamical properties of Φ.
The sign of the jump of Ψ at its discontinuity points x k will be particularly important in our developments. We will distinguish the following cases: (C4) We say that the map is non-overlapping if (C1), (C2) and (C3) are satisfied, and moreover: When the inequality (4.7) does not hold, we identify another case: (C4') We say that the map is overlapping if conditions (C1), (C3) and either (C2) or (C2') are satisfied and Ψ has a negative jump: The terminology follows [27] and refers to the property that Φ is injective in [β, α] under assumption (C4), while the images of (β, w 1 ) and (w 1 , α) under Φ have non-empty intersections (overlap) under assumptions (C4'). We also emphasize that in the non-overlapping (resp., overlapping) case, the map Ψ has non-negative (resp., negative) jumps at its discontinuity points (x k ) k∈Z , i.e., Ψ( ). These conditions may seem complex to check theoretically since they involve relative values for the adaptation map Φ, the discontinuity points, and α and β. However, they are very easy to check numerically for a specific set of parameters. In Fig. 4.1 we illustrate these different situations for a quartic model with a particular choice of the subthreshold parameters and for a fixed value of the reset voltage v R , and we identify the regions with respect to the reset parameters γ and d where the above conditions are satisfied.
We will provide an exhaustive description of the MMO patterns produced by the adaptation map when it has exactly one discontinuity in the interval I. In our framework, we can classify MMOs with half-oscillation precision. However, in this section, we choose for the sake of simplicity in the formulation of the results to consider integer numbers of small oscillations; that is, the points in (β, w 1 ) correspond to no small oscillations whereas the points in (w 1 , α) result in one small oscillation. Thus, referring to the signature of MMOs, we have s i = 0 or s i = 1 and by grouping together in the signature consecutive spikes followed by no small oscillations, we can assume that s i = 1 for any i.
We start with a simple remark stating, roughly speaking, that MMOs occur frequently in our system: Proof. The first two statements about MMOs follow from the monotonicity of Φ in [β, w 1 ) and its limits at w 1 ; indeed, it is easy to see that because of these properties the considered orbits recurrently visit the set (w 1 , α], whereas any point of the orbit in this set undergoes one small oscillation before firing a spike. Hence persistent MMOs result, which are regular if these orbits are periodic. Similarly under all assumptions (C1), (C2) and (C3), Φ is monotone increasing on (w 1 , α] in addition to [β, w 1 ), with Φ(w − 1 ) > w 1 > Φ(w + 1 ). Hence, if an orbit of Φ| I is not trapped by a fixed point in one of these intervals, then it necessarily escapes to the other. In particular, any non-trivial periodic orbit thus features small oscillations as well as consecutive spikes with no small oscillations in between, leading to MMBOs.
Note that under (C1), (C2'), (C3) it is possible that there are periodic orbits fully contained in (w 1 , α]. Such periodic orbits always display one small oscillation before each spike. Hence these are MMOs but not MMBOs. When conditions (C1), (C2) and (C3) are satisfied, the singular case Φ(β) = Φ(α) can be treated using the classical Poincaré theory of orientation preserving circle homeomorphisms. In all other cases the corresponding lift is discontinuous and possibly non-monotonic. Our study will build upon previous works of Keener [27], Misiurewicz [39], Rhodes and Thompson [42,43] and Brette [6]. We link their general results to MMOs in our system, as well as extend and strengthen some of them to more specific subcases arising in our study, allowing for more refined characterization of the dynamics of Φ.
4.1. Non-overlapping case. We start by investigating the non-overlapping case (C4). In that situation, the lift Ψ is discontinuous (unless Φ(β) = Φ(α)) but conserves the orientation-preserving property since it only admits positive jumps. It is well-known that monotone circle maps conserve the properties of continuous orientation-preserving maps: they have a unique rotation number, and rational rotation numbers imply asymptotically periodic behaviors.
To ensure convergence towards periodic orbits, one needs to take special care about the presence of discontinuities. Indeed, when Φ has a periodic orbit with period q, then necessarily there exists x 0 ∈ R such that Ψ q (x 0 ) = x 0 + p(α − β) for some p ∈ N * , p, q relatively prime, i.e. x 0 is periodic mod(α − β) for the lift Ψ. However, since map Φ is discontinuous at w 1 , it might happen that, although the rotation number is rational, no truly periodic orbit of Φ exists but point w 1 acts as a periodic point. This means that one of the two following properties is necessarily fulfilled, with x 0 mod |I| = w 1 (see [42]): • for all t ∈ R, Ψ q (t) > t + p|I| and • for all t ∈ R, Ψ q (t) < t + p|I| and Remark. By allowing the lift to be bivalued at the discontinuity points, Brette [6] and Granados et al [14] avoid the distinction of the three cases for rational rotation numbers (i.e., the existence of the actual periodic orbit and the two cases listed above). That formalism indeed ensures that the periodic orbit always exists, since the two situations above can happen only if x0 mod |I| = w1, i.e. when the periodic orbit bifurcates.
For simplicity, with a little abuse of terminology, in both above cases, we will refer to the orbit of x 0 (mod|I|) under Φ as the periodic orbit. Bearing that in mind we now relate the orbits of Φ to the dynamics of the neuron model and show that the rotation number in the non-overlapping case fully characterizes the signature of the resulting MMO. Theorem 4.3. We assume that the adaptation map Φ satisfies condition (C4) and consider its lift Ψ : R → R. Then the rotation number := (Ψ, w) of Ψ exists and does not depend on w ∈ R.
Moreover, the rotation number is rational, = p/q ∈ Q with p ∈ N 0 := {0, 1, 2, ...} and q ∈ N * relatively prime, if and only if Φ has a periodic orbit, which is related to the MMO pattern fired in the following way: (i) If = 0, then the model generates tonic asymptotically regular spiking for every initial condition w 0 ∈ [β, α] \ {w 1 } (see Figure 4.2, top). (ii) If = 1, then the model generates asymptotically regular MMOs for every initial condition w 0 ∈ [β, α] \ {w 1 }, with periodic signature 1 1 1 1 1 1 ... = (1 1 ). (iii) If = p/q ∈ Q \ Z (p, q relatively prime, q > 1 and 1 ≤ p < q), then the model generates asymptotically regular MMBOs for every initial condition w 0 ∈ [β, α] \ {w 1 } (see e.g., Figure 4.2, bottom). Defining 0 < l 1 < · · · < l p ≤ q − 1 as the unique integers such that l i p/q mod 1 ≥ (q − p)/q and L i = l i+1 − l i for i = 1 · · · p (with the convention l p+1 = q + 1), the MMBO signature is L 1 1 · · · L 1 p . (iv) If ∈ R \ Q, then there are no fixed point and no periodic orbit, and the system fires chaotic MMOs. Remark. This result establishes that in the non-overlapping case, the MMO signatures of orbits are determined by the rotation number and provides a constructive algorithm to compute the MMO signature associated to a given rotation number. We illustrate this construction on two examples: • Orbits of the adaptation map with rotation number = 1/q have signature q 1 . When = (q − 1)/q the signature is 2 1 1 1 · · · 1 1 (with q − 2 repetitions of the pattern 1 1 ).
• For the rotation number = 3/5, up to cyclic ordering, the periodic orbits are ordered as those of the corresponding rotation by 3/5 on the unit circle, i.e. {0, 3 5 , 2·3 5 = Proof. Since the induced lift Ψ : R → R is strictly increasing, we can apply the theory of monotone circle maps theory developed by Rhodes and Thompson [42,43] and Brette [6]. The existence and uniqueness of the rotation number is shown in [42, Theorem 1] and [6], and the proof for orientation preserving homeomorphisms applies 7 . The characterization of the orbits in the case of rational rotation numbers results from [42, Theorem 2] and the fact that Ψ is strictly increasing.
Moreover, if = p/q, then it can be shown that every non-periodic point w ∈ [β, α] of Φ tends under Φ q to some periodic pointw ∈ [β, α]: lim n→∞ Φ nq (w) =w. This is a consequence of [6, Proposition 5] since the monotonicity of Ψ ensures that the underlying circle map is strictly orientation preserving. From the proof therein it also follows that the asymptotic behavior is consistent for all the points of a given orbit, i.e. that if w tends under Φ q tow, then Φ k (w), k = 0, 1, ..., q − 1, tends to its corresponding point Φ k (w) on the periodic orbit ofw. This provides the classification of orbits for the adaptation map, analogous to the one for a circle homeomorphism with rational rotation number (cf. [26,Proposition 11.2.2]). Next, we consider the subcases of firing patterns.
(i-ii) When (Ψ) = 0 mod 1, the adaptation map admits a fixed point. Moreover, under the current assumptions and the way we have defined the lift Ψ we either have (Ψ) = 0 if the fixed point belongs to (β, w 1 ), in which case there is no (full) small oscillation between spikes, or (Ψ) = 1 if the fixed point belongs to (w 1 , α), in which case the orbit displays one small oscillation between every two consecutive spikes.
(iii) As mentioned in Proposition 4.2, periodic orbits necessarily correspond to MMBO. Moreover, it is not hard to show that q-periodic orbits with rotation numbers p/q have exactly p points to the right of w 1 . These points split the periodic orbit into firing events consisting of either one spike or a burst, separated by a small oscillation. Since the lift preserves the orientation, the consecutive points of a periodic orbit {w, Φ(w), . . . , Φ q−1 (w)} with rotation number p/q are ordered as the sequence of numbers (0, p/q, 2p/q, ..., (q − 1)/q) in [0, 1] (up to cyclic permutation, see e.g. [26, Proposition 11.2.1]). The signature of the MMBO is directly related to the indexes l ∈ {0, 1, ..., q −1} such that Φ l (w) > w 1 , and hence such that lp/q ≥ (q −p)/q mod 1. We easily conclude that the signature of the MMBO indeed is L 1 1 L 1 2 · · · L 1 p . (iv) If the rotation number is irrational, then Φ admits no periodic orbit, and all orbits under Φ have the same limit set Ω, which is either the circle or a Cantor-type set as in the continuous case (Φ(β) = Φ(α)), as proved in [6,Proposition 6].
When w 1 is periodic mod (α − β), the corresponding forward attracting periodic orbit is unique. Otherwise, several attracting periodic orbits may exist with the same rational rotation number, and hence with the same period and the same ordering. In [14], the authors have proved the uniqueness of the periodic orbit of maps such as Φ in the non-overlapping case with the assumption that Φ is contractive on both (β, w 1 ) and (w 1 , α). Here, because of the divergence of the differential at the discontinuity points, we cannot use the contraction assumption.
We emphasize that since Ψ is a strictly increasing lift of a degree-one circle map, changing its value at a discontinuity point (while conserving monotonicity) does not change the value of the rotation number (see e.g., [42]). The above remark means that for the characterization of the dynamics of Φ, it does not matter whether we define the lift Ψ to be left-or right-continous at its discontinuity points β + k(α − β), nor that Φ is formally not defined at w 1 , since lim We now provide a simple sufficient condition for the existence of 2 1 MMBOs. This result is analogous to [27, Lemma 3.2] but does not necessitate the boundedness assumption on the derivative of the map made in [27], which our map Φ obviously does not satisfy.
When the second assumption of Proposition 4.4 is not valid and Φ(β) > Φ(α) > w 1 , the dynamics may generate complex orbits of higher period or even chaos. Different MMBO patterns may therefore be observed in the non-overlapping case, depending sensitively on the parameters. We now focus on this dependence on the reset parameters (d, γ) and show that the rotation number varies as a devil's staircase (in the sense of Theorem 4.5 below). This result is based on a theorem in [6]. However it does not follow from [6] immediately, since varying reset parameters changes the invariant interval [β, α] and one needs to add some technical assumptions to ensure that the lifts display an increasing relation. The detailed proof can be found in the Appendix.  A similar result holds for the dependence of the rotation number on the parameter γ in the regime where we can ensure the suitable monotonicity of γ → γ . Fig. 4.4 illustrates a case where this theorem applies.

4.2.
Overlapping case. Now let Φ satisfy the properties of the overlapping case (C4'). In this case, the lift Ψ is no longer increasing (it has negative jumps at the points x k = α + k(α − β) for k ∈ Z), and a number of important properties inherited from the well-behaved dynamics of orientation-preserving circle homeomorphisms that persist in the non-overlapping case [6,42,43] are now lost, leaving room for still richer dynamics.
In the overlapping case, it is easy to see that our map restricted to its invariant interval [β, α] falls in the framework of the so-called old heavy maps [39], since it is a lift of a degree-one circle map with only negative jumps. These maps have interesting dynamics with non-unique rotation numbers. More precisely, we can define a rotation As noted in [39], these two quantities are the (unique) rotation numbers of the continuous orientation preserving maps: We can now conclude after [39]: Proposition 4.6. Under assumption (C4'), 1. if Φ admits a q-periodic point w with rotation number (Ψ, w) = p/q, then a(Ψ) ≤ p/q ≤ b(Ψ); 2. if a(Ψ) < p/q < b(Ψ), then Φ admits a periodic point w of period q and rotation number (Ψ, w) = p/q. In both cases, the orbit of w displays MMOs, unless p = 0, in which case w is a fixed point in [β, w 1 ). Moreover, for any 1 and 2 such that a(Ψ) This result implies in particular that the rotation set in the overlapping case is closed, meaning that every number ∈ [a(Ψ), b(Ψ)] is the rotation number (Ψ, w) of an orbit with initial condition w ∈ [β, α], and the rational numbers in its interior correspond inevitably to periodic orbits. The property of having a non-trivial rotation interval implies coexistence of infinitely many periodic orbits of distinct periods; this situation is sometimes referred to as 'chaos' (see [27]), although this notion differs from the chaos associated with non-regular ('chaotic') behavior of orbits with non-rational rotation numbers. We next consider (i) the variation of the rotation interval as a function of the reset parameters and (ii) the relationship between rotation intervals and MMO patterns.
A result from [39,Theorem B] ensures continuous dependence of the boundaries of the rotation interval a(Ψ) and b(Ψ) on map parameters when Ψ l and Ψ r , regarded as elements of the space C 0 (R) occupied with the uniform topology, depend continuously on these parameters. The following proposition makes this dependence more precise in our case by showing that these vary as a devil's staircase under mild assumptions.
Proposition 4.7. Consider fixed parameters v R , a, b, γ and I and vary d ∈ [λ 1 , λ 2 ] such that, for each d ∈ [λ 1 , λ 2 ], the corresponding adaptation map Φ d satisfies the assumptions of the overlapping case (C4'). Then the maps d → a(Ψ d ) and d → b(Ψ d ) assigning to d the endpoints of the rotation interval of Φ d are continuous.
If we further assume that, for any pair then the maps d → Ψ d (w), d → Ψ d,r (w) and d → Ψ d,l (w) are increasing for each w (Ψ d,r and Ψ d,l denote, respectively, upper and lower enveloping maps of the lift We note that the sufficient condition (4.17) is equivalent to where Ψ d2 (β + d2 ) denotes the right limit of Ψ d2 at β d2 . This latter condition is satisfied for instance when, for every d ∈ [λ 1 , λ 2 ], Φ d < 1 in the whole interval [β d , β d + (λ 2 − λ 1 )]. The proposition is proved in Appendix A and illustrated in Fig. 4.6. In the overlapping case, the pointwise rotation numbers (Ψ, w) may not exist for some initial conditions and generally depend upon w, implying that we have non-trivial (i.e., non-singleton) rotation intervals. Moreover, despite Proposition 4.6, even knowing the rotation interval [a(Ψ), β(Ψ)] does not fully determine yet the structure of the set of all (minimal) periods of orbits of Φ. Specific cases were fully characterized, however, notably degree-one continuous non-injective circle maps [3,40]. For maps with discontinuities, the issue is very complex and, to our knowledge, periods of periodic orbits are completely described only for lifts of monotonic modulo 1 transformations (see [18]), which corresponds to the overlapping case with the additional monotonicity assumption (C2). In addition to these difficulties, the overlap prevents systematic deduction of the MMO signature from knowledge of the rotation number, since the rotation number (Ψ, w) does not determine the ordering of points on the orbit of w.
Specific analysis on the maps considered here, however, provides some information about this characterization. First, when (C2) holds, the map Φ is piecewise increasing and thus Proposition 4.2 applies and ensures that periodic orbits are associated to MMBOs. When (C2) is not valid, we still know that periodic orbits fully contained in (w 1 , α) correspond to MMOs with signature 1 1 . Beyond these particular cases, we now demonstrate more general results under milder assumptions.
Proof. First of all, we note that since w f ∈ [β, w 1 ) is a fixed point of Φ, it is also a fixed point for Ψ and thus the associated rotation number is equal to 0. Moreover, for w f ∈ (w 1 , α) we have Ψ(ŵ f ) =ŵ f + (α − β) and thus the associated rotation number under Ψ is equal to 1. We thus conclude that the rotation interval of Ψ contains the full interval [0, 1], which concludes the proof. Proposition 4.9. Assume that Φ satisfies (C4') and that Φ admits at least one fixed point in (w 1 , α), the smallest of which we denote by w f ∈ (w 1 , α). Assume moreover that there is no fixed point of Φ in [β, w 1 ). Then • if Φ(β) < w f , then there existsq > 1 such that for all q >q, Φ admits a periodic orbit of period q, and the associated trajectories display MMOs; and periodic orbits correspond to MMOs with signature 1 1 . If additionally α ≤ w * , then Φ admits no periodic orbit of period q > 1, every orbit converges towards a fixed point in (w 1 , α], and associated trajectories display asymptotically regular MMOs with signature 1 1 . Proof. We first assume that Φ(β) < w f . In this case, the lower envelope Ψ l intersects neither the identity (Id) line nor the Id + (α − β) line (and, obviously, none of the lines Id + k(α − β) for k ∈ Z). Thus the graph of Ψ l is fully contained between the lines Id and Id + (α − β) and since the functions Ψ l (w) − w + k(α − β), k ∈ Z, are continuous and α − β periodic, there exists δ > 0 such that for every w and n ∈ N * we have Ψ n l (w) < w + n(α − β) − nδ and .
On the other hand, Ψ r (w f ) = w f + (α − β) and thus b(Ψ) = (Ψ r ) = 1. Therefore, the rotation interval is not trivial and For every q > 1 large enough, we have that is, there exists a periodic orbit of Φ with period q and rotation number q−1 q . We now assume Φ(β) ≥ w f . Then w f is also a fixed point mod (α − β) of Ψ l , i.e.
Thus, if there was some periodic orbit of period q > 1, all points of such an orbit would lie in (w 1 , α) and would have rotation number 1, yielding MMOs with signature 1 1 . However, if additionally α ≤ w * , then the map Φ is increasing in (w 1 , α) and no periodic orbit can be fully contained in this interval.
Assuming that α ≤ w * , we notice that the interval [w f , α] is invariant for Φ and that Φ is continuous and increasing therein. Consequently, every point w ∈ [w f , α] tends under Φ to one of the fixed points in [w f , α]. But as every point in [β, w f ) in mapped finally into [w f , α], this holds for all the points in [β, α] and the proof is completed.
Later, we shall complement the above result in a slightly more general situation, in Theorem 4.12, by treating maps admitting fixed points in [β, w 1 ) and lacking a fixed point in (w 1 , α]. We can also easily justify the following: Corollary 4.10. In the overlapping case, the existence of a fixed point of Φ and of a periodic orbit with rotation number p/q, for some period q > 1 and p = q, implies the existence of periodic orbits with all arbitrary periods greater than q, each yielding MMOs. In In contrast to the non-overlapping case, in the overlapping case we have dropped the assumption (C2) that the map is piecewise increasing. However, under this assumption we can describe the chaotic behavior of the map's iterates more precisely: Corollary 4.11. Assume that Φ satisfies (C4') with (C2), that Φ(α) < w 1 , that Φ has at least two periodic orbits with periods q 1 = q 2 and that exactly one point of each of these periodic orbits is greater than w 1 . Then the mapping w → Φ(w) is a shift on a sequence space.
To obtain the above proposition it suffices, for example, to look at the proof of [27, Theorem 2.4] and notice that the piecewise contraction assumption made in [27] does not interfere in the proof of this particular theorem and therefore it extends to our class of discontinuous maps with unbounded derivative.

A general result for adaptation maps with one discontinuity in the invariant interval [β, α].
In previous sections we have classified the dynamics of the adaptation map and the associated spiking patterns in terms of rotation numbers and rotation intervals. However, for particular values of the parameters, we lack an explicit analytical expression with which to characterize the corresponding rotation properties. Below we address this problem under the assumptions (C3), that [β, α] is an invariant interval, and (C1), that the map has a unique discontinuity point within this interval, regardless of whether the map is in the overlapping-or non-overlapping case or neither of these (e.g. when the jumps at β + k(α − β) are positive but there is an overlap in values of Φ| [β,w1) and Φ| (w1,α] ).
Theorem 4.12. Assume that conditions (C1) and (C3) hold and that Φ has a fixed point in [β, w 1 ). By w f denote the largest fixed point in [β, w 1 ). Then 1. if max{Φ(w) : w ∈ (w 1 , α]} < w f , then the rotation number (Ψ, w) = 0 is unique and the system displays no MMOs; 2. if max{Φ(w) : w ∈ (w 1 , α]} ≥ w f > β, then there are subintervals of [β, α] of points with rotation number 0, corresponding to orbits with no MMOs. However, if simultaneously Φ(α) ≥ Φ(β), then there existsq ∈ N * such that for every q ≥q, Φ admits also a periodic pointŵ ∈ (w f , α) of period q, displaying MMOs. Proof. The first result follows from the fact that every point w ∈ [β, α] \ {w 1 } is mapped into [β, w f ] after at most a few iterates, and, since Φ([β, w f ]) ⊂ [β, w f ] and Φ is increasing therein, it is eventually attracted to one of the fixed points located in For the second result, the same argument applies to show the existence of subintervals of [β, α] with rotation number 0 under the assumptions made. To establish that there isq ∈ N * such that periodic points of every period greater thanq exist under the additional assumption that Φ(α) ≥ Φ(β), it suffices to show that the rotation interval is of the form [0, δ] for some δ > 0. Since w f is a fixed point for the lower enveloping map Ψ l , we clearly have a(Ψ) = (Ψ l ) = 0. In contrast, the upper enveloppe Ψ r has no fixed points, and using a similar argument as in the proof of Proposition 4.9, we show that b(Ψ) = (Ψ r ) ≥ δ α−β > 0 for some δ > 0, which completes the proof.    Fig. 4.1 are superimposed on the colormap. Regions A, B and C comprise the non-overlapping case, i.e. assumption (C4) is fulfilled, and general Theorem 4.3 applies for (d, γ) values in these regions. In particular, the rotation number of Φ is unique, i.e. does not depend on the initial condition.
• In region A, Φ(α) < Φ(β) < w 1 . Along certain paths in this region, Theorem 4.5 applies and the rotation number varies as illustrated in Fig. 4.4. • In region B, Φ(α) < w 1 < Φ(β), hence Proposition 4.4 applies and ensures the existence of a period-2 orbit of Φ, with rotation number equal to 1/2. • In region C, w 1 < Φ(α) < Φ(β). This region may feature a variety of different dynamics including all types of behavior arising in the other regions. In the example in the right panel of Fig. 4.4, the unique rotation number is 1/2, but this value depends on the choice of (d, γ), as can be seen in the left panel. Regions D and E comprise the overlapping case and Φ may admit different rotation numbers depending on the initial condition. For (d, γ) in these regions, the lift Ψ associated with the adaptation map exhibits only negative jumps. The general Proposition 4.6 applies, which ensures the existence of a rotation interval. Using the left and right lifts Ψ l and Ψ r associated with Φ, one computes the endpoints of the rotation interval and their evolution according to parameter d (Proposition 4.7 and Fig. 4.6).
• In region D, α < w * and Φ is piecewise increasing. The rotation number is not uniquely defined in the general case. Nevertheless, along the particular chosen path in the parameter space (d, γ) shown in the right panel of Fig. 4.4, Ψ l and Ψ r present the same rotation number 1/2 and the rotation number of Φ does not depend on the initial condition. This particular simulation illustrates a way to demonstrate that the rotation number of the adaptation map is unique by showing that the rotation interval is reduced to a singleton. • In region E, w * < α. The rotation numbers of Ψ l and Ψ r differ and the rotation interval of the adaptation map varies with changes in (d, γ) within the region bounded by the black and red lines in the right panel of Fig. 4.4. Note that the global Theorem 4.12 applies in all regions A to E. One may track the appearance and disappearance of the fixed points according to the values of d and γ together with the evolution of the rotation number or rotation interval. Outside of regions A to E, the structure of the lift is more complex due to the presence of additional discontinuity points. Yet, the numerical calculation of the rotation number can be performed for a given initial condition.

5.
A note on the case of two or more discontinuities. One challenge in this study is related to the fact that the map under scrutiny, the adaptation map, is not known analytically. Our mathematical analysis has covered in detail the cases of overlapping and non-overlapping maps with one discontinuity in the invariant interval. These situations do not cover all possible shapes of adaptation maps that can induce lifts with more discontinuity points; indeed, multiple discontinuities can yield a combinatorial explosion of cases with different combinations of possible jumps as well as maps that are non-monotone but with only positive jumps. While in these cases it is still possible to obtain upper and lower bounds for the rotation set by computing the rotation numbers of the non-decreasing maps Ψ l and Ψ r , defined in the same way as in the overlapping case (C4'), it remains an open question to determine when every value within this interval corresponds to the rotation number of a given orbit, and it is not hard to find elementary examples for which this is false 8 . Thus the general, complete and precise characterization of the dynamics of the system is a complex and rich mathematical problem that raises several deep questions of iterates of interval maps with discontinuities. In particular, we have seen that in the non-overlapping case (C4), the rotation number allowed us to completely decode the MMO signature. A natural extension of this work is thus to define for maps with more discontinuities a mathematical invariant (perhaps some vector of numbers) that would either provide the exact signature of each supported MMO pattern or allow calculation of how many points from a random orbit would be expected to fall into each continuity interval and hence how frequently a given number of small oscillations occurs between two consecutive spikes.
Let us conclude with the following exemplary result, which allows for multiple and even infinitely many intersections w i of the reset line {v = v R } with W s , assuming that only finitely many of them lie in the interval (β, α): Theorem 5.1. Suppose that Φ(β) > β and that there are finitely many discontinuity points of the map Φ in (β, α), all located in (β, w * ). Then the adaptation map Φ induces the rotation interval with the same properties as in Corollary 4.6. In particular, to every rational rotation number in the interior of this interval, there corresponds a periodic orbit, displaying regular MMOs.
The above theorem is straightforward once it is noted that the suitably defined lift Ψ for the adaption map under the given assumptions is an old heavy map, as the maps studied in [39]. Therefore, in particular, one can also derive conditions e.g. for periodic orbits of all possible periods exhibiting MMOs (with richer structure than what we considered earlier due to the additional discontinuities), and the corresponding regions in the space of reset parameters for specific models can be computed numerically, in the same way as in the previous subsection. We emphasize that due to the properties of the adaptation map (in particular, the fact that Φ(w i , w i+1 ) = (β, α) for the consecutive discontinuity points w i , w i+1 < w * ), the lift Ψ of Φ is very likely to be an old heavy map. Typically, in case of multiple discontinuities one can expect the rotation interval to cover the whole interval [0, 1] and the occurrence of periodic orbits of all periods and rich MMO structure.
6. Discussion. Nonlinear bidimensional hybrid neuron models, which combine continuous subthreshold dynamics with a spike-related jump or reset condition, are easily defined and show an astonishingly rich mathematical phenomenology. A number of studies have already revealed their subthreshold dynamical properties [57], investigated their spike patterns in the absence of any equilibrium state of the subthreshold dynamics [60], and highlighted their versatility [7,20,52] and capacity to reproduce neuronal dynamics [21,23,41,59]. The present paper and its companion [49] add to this body of works by studying (i) chaotic dynamics and period-incrementing structures, and (ii) oscillating solutions associated with multiple unstable equilibria. The latter led us to investigate the dynamics of a particular class of interval maps that feature both discontinuities and divergence of the derivative. Interestingly, in the presence of an unstable focus of the subthreshold dynamics, we have shown that the spike patterns fired may correspond to complex oscillations that combine action potentials (or bursts of action potentials) and subthreshold oscillations, trajectories known as MMOs or MMBOs in continuous dynamical system.
In contrast to continuous dynamical systems, these forms of complex oscillations can occur in hybrid systems with only two variables. Moreover, the mechanism of generation of these trajectories differs between these two models; in the hybrid case, MMOs result entirely from the topology of the invariant manifolds of the continuoustime dynamics. As such, these trajectories can occur in systems that lack timescale separation and based on a mapping approach, discrete dynamical systems methods can be used to rigorously establish their existence and properties. One may however wonder if there exists a relationship between the two systems, and particularly it is tempting to interpret the hybrid system as the reduction of a differentiable multiple timescale system in a certain singular limit. The wide variety of MMOs (in particular the wild signatures encountered) produced with the reset mechanism indicates that such a differentiable system should be at least four-dimensional and the vector field should induce a highly complex return mechanism within the region of the phase space where small oscillations are generated (funnel). The construction of such a return mechanism for reproducing the same versatility in the MMOs signature in the differentiable case remains a challenging problem from the dynamical viewpoint, involving complex interactions between the different timescales.
To tune the model parameters to attain the regime studied in this work, we introduced a parameter γ, which yields an attenuation of the adaptation variable during the reset. This adjustment to the reset mechanism accounts for the durations of spikes fired (see [49]). With this new parameter, the quartic model (and, we expect, all other models of the class, including the Izhikevich model [20] and the adaptive exponential [7]) can be tuned to achieve any of the cases we have identified. Therefore, our analysis provides useful information for tuning model parameters to achieve outputs fulfilling a list of qualitative and quantitative specifications. In particular, the ability to reproduce fine trajectories of MMOs may be useful when modeling neurons in situations in which synchronization is essential. Indeed, in neuroscience, it has been shown that in the presence of noise, small subthreshold oscillations support the generation of precise and robust rhythmic spike patterns, as recorded in specific rhythmic pattern generators such as the inferior olive nucleus [5,34,35], in the stellate cells of the entorhinal cortex [1,2,25], and in the dorsal root ganglia [4,32,33]. A possible direction for future work would be to go deeper into the analysis of the shape of the adaptation map of the adaptive exponential integrate-and-fire system to relate the presence and possible signature of MMOs to variations in biophysical parameters, following e.g. [59].
Another important direction related to the roles of model parameters would be to characterize the structural stability of trajectories and their possible bifurcations. First works in that direction have been developed in [11]: taking into account the infinite contraction of the trajectories in the voltage variable associated with the reset, the authors proposed to compute expansion or contraction exponents along transverse directions, providing a notion of stability of hybrid orbits that is more explicit than criteria on the shape of the adaptation map. It would be interesting to develop these methods in the cases of non-monotonic spiraling trajectories associated with the presence of MMOs. Alternatively, using models with simpler subthreshold dynamics, for instance linear or piecewise linear [24,46], may allow for a derivation of an explicit expression of the reset maps, thus for fine characterization of the stability of the orbits.
At the level of the adaptation map, a question that is open in the overlapping case is to characterize the stability of orbits when the system has multiple possible rotation numbers. Indeed, even if the rotation interval is not a singleton, one often observes in simulations that only one rotation number is actually realized. There are two typical reasons why this could occur: either there is an attracting periodic orbit that attracts most initial conditions or the system has an invariant measure µ, absolutely continuous with respect to the Lebesgue measure, in which case the observed rotation number is just the average displacement Ψ(w) − w with respect to the measure µ. Nonetheless, rigorously establishing the existence of such a measure is a challenge in most of systems arising from applications. In particular, we cannot use e.g. the classical Lasota-Yorke theorem ( [31]), since the derivative Φ diverges at the discontinuity points. On the other hand, for investigating stability of orbits a possible approach would be to use and develop symbolic dynamics and kneading theory for such discontinuous interval maps. However, we emphasize that in our characterization of the orbits and the patterns of complex oscillations fired, rotation theory turned out to be the most useful tool since we have a unequivocal, bidirectional link between the rotation number and the signature of the MMO (Theorem 4.3), which allows us to characterize situations in which the neuron shows regular spiking, MMO, bursting, MMBO or chaotic behavior.
In these studies, we have made a crucial use of the planar nature of the system. MMOs will of course exist in higher dimensional hybrid dynamical systems, and analysis would require fine characterization of the invariant manifolds. The extension of the theory to higher dimensional systems would be particularly interesting from the computational neuroscience viewpoint for understanding the behavior of neuron networks in which several neurons driven by such dynamics are coupled and communicate at the times of the spikes.
Moreover, we also observe that for any d ∈ [d 1 , d 2 ], the maps Φ d have the same discontinuity point w 1,d , and the lifts Ψ d are continuous at points w 1,d + k(α d − β d ), have positive jumps at α d + k(α d − β d ) and satisfy Ψ d (w + θ) = Ψ d (w) + θ. So in fact all these lifts Ψ d can be seen as lifts of non-continuous invertible circle maps under the same projection p : t → exp( 2πıt θ ). However, even if the map Φ d is increasing with d, this is not necessarily the case for Ψ d , because of the simultaneous fluctuation of the invariant interval. Indeed, when each lift Ψ d is obtained from Φ d [β d ,α d ] the relation Ψ d1 (w) < Ψ d2 (w) for d 1 < d 2 might be violated in the intervals [β d1 , β d2 ], as at the point α d2 we glue the right part of the graph of Φ d2 | [β d 2 ,α d 2 ] to its left part (shifted up by θ). But noticing that under the additional condition Φ d (α d2 ) < Φ d (β d1 ) for any d ∈ [d 1 , d 2 ], the interval [β d1 , α d2 ] constitutes a particular invariant interval in which the adaptation map Φ d is piecewise increasing and non-overlapping, we can build well-behaved liftsΨ d : R → R based on the shape of the map Φ d on this bigger invariant interval [β d1 , α d2 ]. In contrast to Ψ d , these new lifts are discontinuous at the points w 1,d + k(α d2 − β d1 ) (where they have positive jumps of amplitude d 2 − d 1 ), in addition to their discontinuity at α d2 +k(α d2 −β d1 ), k ∈ Z. The latter jump also remains positive under our assumption that Φ d (β d1 ) is strictly greater than Φ d (α d2 ). Constructing lifts Ψ d on an enlarged invariant interval [β d1 , α d2 ] instead of [β d , α d ] has the advantage of ensuring that the mapping (w, d) →Ψ d (w) is increasing in both variables. Moreover, it has no effect on the dynamics, since any orbit of Φ d with an initial condition in [β d1 , α d2 ] enters after a few iterations into the interval [β d , α d ]. Since the orbits {Ψ n d (w)} project mod (α d2 − β d1 ) to the orbits {Φ n d (w)}, we therefore have (Ψ d ) = (Ψ d ). Concluding the proof therefore only amounts to showing that the map d →Ψ d is continuous in the Hausdorff topology, which is very simple once it is noted, as mentioned above, that this property is equivalent to the uniform convergence at all points in the interior of [β d1 , α d2 ] \ {w 1,d } and thatΨ d −Ψ d = d − d on this interval. Thus the mappingρ : d → (Ψ d ) has the properties listed in the theorem (compare with Theorem 2 in [6]) and consequently, the same holds for ρ : d → (Ψ d ).
We have noticed that while continuity of the lifts under the Hausdorff topology was always satisfied in our case, an additional assumption is necessary to ensure that the mapping (s, w) → Ψ s (w) (where s denotes a parameter, here d or γ) is increasing in both variables, which otherwise is not always true. We emphasize that even in situations in which this mapping is not increasing in both variables, the rotation number remains continuous under the H-convergence provided that the limit function Ψ s0 is strictly increasing, see [43,Proposition 5.7].
The plateaus of rotation number observed in the devil's staircase situation are a general property of our system, called locking (see [43] for precise definition of locking).
Remark. We observe that no condition beyond monotonicity of the lifts in w and d is required to show locking of the rational rotation number in the strictly non-overlapping case (i.e. Φ(α) < Φ(β)), unlike the case of continuous circle maps. When Φ(α) = Φ(β), the lift Ψ would be in fact a lift of an orientation preserving circle homeomorphism and thus locking of the rotation number at rational values requires that there is no conjugacy with rational rotation for such a map (see e.g. Propositions 11.1.10 and 11.1.11 in [26]).
Proof of Theorem 4.7 The first part of the proof amounts to showing that the upper and lower envelopes of Ψ d , denoted Ψ d,l and Ψ d,r , are uniformly continuous in d for d ∈ [λ 1 , λ 2 ].

Remark.
To ensure that the mapping t → (Ft) behaves as a devil's staircase for a continuous increasing family {Ft} t∈[T 1 ,T 2 ] of continuous non-decreasing degree-one maps Ft, we also need to make sure that there exists a dense set S ⊂ Q such that, for s ∈ S, no map Ft is conjugated to the rotation Rs by s and that the map t → (Ft) is not constant (see Proposition 11.1.11 in [26]). However, in practice, these two specific cases do not occur for any of the envelopes Ψ l and Ψr of the adaptation map.