Complex oscillations in the delayed Fitzhugh-Nagumo equation
CComplex oscillations in the delayed FitzHugh-Nagumo equation
Maciej Krupa †§ Jonathan D. Touboul ∗§ June 25, 2015
Abstract
Motivated by the dynamics of neuronal responses, we analyze the dynamics of the FitzHugh-Nagumoslow-fast system with delayed self-coupling. This system provides a canonical example of a canard explosion forsufficiently small delays. Beyond this regime, delays significantly enrich the dynamics, leading to mixed-modeoscillations, bursting and chaos. These behaviors emerge from a delay-induced subcritical Bogdanov-Takensinstability arising at the fold points of the S-shaped critical manifold. Underlying the transition from canard-induced to delay-induced dynamics is an abrupt switch in the nature of the Hopf bifurcation.
Keywords
Delayed Differential Equations; Slow-Fast systems; Mixed-Mode Oscillations; Bursting; Chaos
Contents ∗ The Mathematical Neurosciences Team, Center for Interdisciplinary Research in Biology (CNRS UMR 7241, INSERM U1050,UPMC ED 158, MEMOLIFE PSL*), [email protected] § MYCENAE Laboratory, Inria Paris-Rocquencourt † NEUROMATHCOMP Laboratory, Inria Sophia Antipolis-M´editeran´ee a r X i v : . [ m a t h . D S ] J un .2 Stability of fixed points on the critical manifold . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173.3 Bogdanov-Takens singularity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 γ (cid:54) = 0 Dynamical systems with multiple timescales and delays are used as models in a number of applicationsincluding the dynamics of excitable cells in biology [16, 11, 41, 43, 42, 40, 46], mechanical systems [6], chemicalreactions [30] and physical systems [25]. Ordinary differential equations with multiple timescales are knownto display a variety of complex oscillations, including canard explosions [2], relaxation oscillation, mixed-modeoscillations (MMOs) [4] and bursting [39]. A canard explosion is a very fast transition, upon variation of aparameter, from a small amplitude limit cycle to a relaxation oscillator, a type of periodic solution consistingof long periods of quasi static behavior interspersed with short periods of rapid transitions. The topic of thiswork and the companion investigation [32] is to understand the effect of introducing delays in a system whosedynamics, in the absence of delays, is characterized by a canard explosion and relaxation oscillations. In thepresent paper we focus on the FitzHugh Nagumo (FhN) system with delays: x (cid:48) t = x t − x t + y t + J ( x t − x t − τ ) y (cid:48) t = ε ( a − x t ) , (1)for which we find a variety of complex oscillations as delays are varied. System (1) is slow/fast, with the ratioof timescales given by a small parameter ε and a delay τ ≥ ε = 0, and then constructing thedynamics of the full system for ε > ε > ε = 0 approximations. To this purpose, one needs to be able to characterizethe behavior of slow manifolds, which classically relies heavily on the use of Fenichel theory [19]. Recent resultsof Hupkes and Sandstede [26] imply the existence of Fenichel slow manifolds in the system we study, and arefundamental to the generalization of geometric singular perturbation theory developed in this paper as well asin the companion investigation [32].Qualitatively, we will see that as delays are increased, system (1) can undergo two types of transition fromsteady-state to oscillatory dynamics: one through a delay-induced canard phenomenon [2, 1, 12], and one duea destabilization of the slow manifold due to the delays in the fast dynamics. The canard explosions found aregeneric: the qualitative behavior of the system in the vicinity of this transition is similar to that found in systemswith no delays. Their characterization in a general context is the subject of our companion study [32]. Here, weconcentrate on one specific system (1) and investigate the region of the parameter space complementary to thecanard explosion, where delay induced instabilities play a significant role in the dynamics. While some resultsabout Hopf bifurcations and reduction to normal form are proved, our study encompasses a part of numericalsimulations and continuation. This leads us to predict and confirm by means of simulations the existence ofa variety of complex oscillatory behaviors, both periodic and chaotic, including bursting orbits correspondingto recurrent dynamics composed of phases of fast oscillations interspersed with periods of quiescence. Othertypes of dynamics relevant to our findings are MMOs, which are recurrent trajectories of a dynamical systemcharacterized by an alternation between oscillations of very distinct large and small amplitude (see [10] for arecent review on the subject). In our work we show the presence of a new type of MMO, due to the presenceof a Bogdanov-Takens point in the fast dynamics.It is interesting to note that the type of complex oscillations (relaxation, MMO and bursting) found inour simple system may be relevant to the underlying biological problem, as similar solutions are reported inneuronal recordings and models [29, 27].The interest of this work is also mathematical, as our study should be seen as a pilot investigation in theeffort to extend geometric singular perturbation theory to systems with delays in the fast dynamics and finitedimensional slow dynamics.The paper is organized as follows. We introduce in section 1 our model and the differences between our3pproach and previous work. We start by investigating, using bifurcation theory and numerical simulations, theoscillatory dynamics of the full system in section 2. We will observe a number of complex oscillatory patterns,which we investigate using the GSPT. We therefore start by investigating the bifurcations in the fast systemin section 3, and complement this analytical study with extensive numerical simulations leading to a morecomprehensive bifurcation diagram presented in section 4. Section 5 presents how these elements combine toproduce complex oscillatory patterns including MMOs, bursting and chaotic behaviors. The study of the dynamics of delay differential equations with multiple timescales has been the subject of anumber of important works. We briefly review these contributions and summarize the main results of the presentpaper in section 1.1. In this study we are interested in one specific system (1) motivated by the modeling of theelectrical activity of nerve cells, which we present in more detail in section 1.2.
Our results are concerned with the dynamics of a planar delay dynamical system (1) with timescale separationand a delay in the fast equation. A similar system has been investigated in a mechanical context [6], where thebehavior of a machine is described by a delay differential equation with a slow variable. The authors focus onthe behavior of the system in the limit of small delays. In that regime they are able to reduce their equationto an ODE and use finite-dimensional methods to show existence of canard explosions. In the present paper,we shall be interested in the behavior of the system beyond the small delay regime, where the reduction toODEs no longer holds. This is where we will need to use the singular perturbation methods to characterize thebehavior of the system.Another class of models including fast dynamics and delays was introduced in the 1980s. First works in thecontext of delay differential equations with multiple timescales were carried out by Mallet-Paret and Nussbaumin [34] and focused on scalar equations with delays of type: ε ˙ x = − x + f ( x ( t − . (2)In that paper, they showed existence of periodic solutions and investigated their behavior in the limit ε → ε = 0 leads to a discrete dynamical system, which constitutes a significant difference withour model, where one obtains a continuous time delayed differential equation in one dimension. In that sense,our system (1) is more directly related to finite dimensional slow-fast systems. We rely heavily on this analogy4nd develop the GSPT in that context. A number of results were obtained in the context of equation (2):the authors of [51, 50] investigated singular Hopf bifurcation, which is classically closely related to canardexplosions [1]. These authors used a method based on asymptotic expansions, which is relatively different andsomehow complementary to GSPT. Finally, in the present paper we do not consider canard explosions (singularHopf bifurcations) and focus on the parameter region where they do not occur. Our companion paper [32] doesfocus on canard explosions and is similar to [51, 50] in this aspect, but rather than asymptotic expansions usesthe complementary GSPT method. Finally we mention that asymptotic expansions were extended to a systemthat has more similarities to our system in [21] in that the equation is no more scalar and the scaling of thedelay is as in our equation (1). The analysis of [21] focuses mainly on the vicinity of Hopf bifurcation and doesnot make a connection between the fast bifurcation diagram and the patterns of the flow, which is fundamentalto our work.In the present paper, we shall demonstrate a few analytical results that we summarize below. These areconcerned with transition to oscillation through Hopf bifurcations. The proofs are classical, yet the conclusionsare very interesting, namely the lack of smoothness of of the Hopf bifurcation curve in the limit ε = 0 and theexistence of a Bautin (generalized Hopf) point, where the sign of the first Lyapunov coefficient changes. Weprove a result on a Hopf and Bogdanov-Takens bifurcation in the fast system and a result on a Hopf bifurcationin the full system. The statement of the theorem relevant to the fast system is as follows: Theorem 1.
Consider the fast system of (1) with
Jτ > .1. There exists a family of disconnected Hopf bifurcation curves indexed by k ∈ N given by: τ kf ( J, a ) := 1 √ a − √ J + 1 − a (cid:18) cos − (cid:18) − a + JJ (cid:19) + 2 kπ (cid:19) , y = a − a / . (3) with a ∈ (1 , √ J + 1) and τ ∈ (1 /J, ∞ ) .2. There are two Bogdanov-Takens points at Jτ = 1 , y = ± . On the full system, we prove the following:
Theorem 2.
Consider system (1) with
J > .1. For ε > there exists a unique Hopf bifurcation curve in the plane ( a, τ ) . The curve is smooth, goes from a = 1 to a = √ J back and forth infinitely many times and reaches arbitrarily large values of τ . Itcan be decomposed into the following branches: τ k ( J, ε, a ) := A + √ A +4 ε (cid:16) cos − (cid:16) − a J (cid:17) + 2 kπ (cid:17) τ k ( J, ε, a ) := A −√ A +4 ε (cid:16) cos − (cid:16) − a J (cid:17) − k + 1) π (cid:17) (4)5 ith A = (cid:112) ( a − J − a ) , a ∈ [1 , √ J ] and k ∈ N . Qualitatively, we have:2. If τ and J are fixed so that Jτ < then, for ε > sufficiently small, there exists a unique Hopf bifurcation τ , it is generic and supercritical.3. If τ and J are fixed so that Jτ > then, for ε > sufficiently small, the Hopf bifurcation correspondingto the change of stability of the fixed point (at τ ) is generic and subcritical.4. As ε → , the curve τ associated to the change of stability of the fixed point converges to the union ofthe line segment a = 1 , Jτ ≤ and the Hopf bifurcation curve of the fast system (see Theorem 1). Theprimary Hopf bifurcation τ in that limit is continuous but not differentiable, and the sub- or super-criticalnature of the bifurcation changes at Jτ = 1 . Note that Theorem 2 covers in part the parameter region
Jτ <
1, so there is some overlap with [32],where existence of Hopf bifurcation is proved using a different approach (center manifold theorem). In [32] theexistence of Hopf bifurcation is just a minor point in establishing the entire picture of canard explosion. Here wehave a very simple and direct method of analysis which allows us to explain the remarkable shape of the locusof Hopf bifurcation in system (1). We note that as a consequence of 2, there exists along the Hopf bifurcationfor ε >
A central model in mathematical neuroscience, encompassing the excitable nature of the cells and simpleenough to allow analytical developments, is given by the FitzHugh-Nagumo model [20, 36]. This equation,initially introduced as a simplification of the very versatile but much more complex Hodgkin-Huxley neuron6odel, describes the excitable dynamics of the membrane potential of the neuron x coupled to a slow recoveryvariable y , through the slow/fast ordinary differential equations: x (cid:48) = x − x + yy (cid:48) = ε ( a + b x − γy ) (5)where a represents the input current the neuron is subjected to and b is the interaction strength between thevoltage and the recovery variable. The small parameter ε represents the timescale ratio between the voltage andthe recovery variables. Note that the classical van der Pol oscillator corresponds simply to the case γ = 0 [38].In [37], a linear feedback was introduced to model recurrent self-coupling. Here, we shall rather motivateour model by the description of the electrical activity of large networks. Indeed, cortical behaviors generallyarise at macroscopic scales involving a large number of neurons. We consider here that neurons are electricallycoupled through gap junctions, resulting in an input current proportional to the voltage difference between thecommunicating cells [15]. We also incorporate the important fact that neurons interact after a delay related tothe transmission of information through the axon and dendrites. The corresponding N -neurons network modewith stochastic input is given by the equations: dx it = ( x it − ( x it ) + y it + JN (cid:80) Nj =1 ( x it − x jt − τ )) dt + σdW it dy it = ε ( a + b x it − γy it ) dt (6)where W it are independent Brownian motions. The presence of noise allows to prove (see [47, 35]) that in thelimit N → ∞ , the behavior of a given neuron in the network satisfies the implicit stochastic equation: dx t = ( x t − x t + y t + J ( x t − E [ x t − τ ])) dt + σdW t dy t = ε ( a + b x t − γy t ) (7)where E [ x t ] denotes the statistical expectation of the process x t . Note that the presence of noise is essentialto derive this limit. In specific cases [48, 49], such mean-field equations can be exactly reduced to ordinarydifferential equations in aggregate variables (e.g., mean and standard deviation). Here, the nonlinear dynamicsof the cells prevents from such a reduction. However, in the zero noise limit (a ‘viscosity solution’, that maybe seen as a fully synchronized solution of the network equation in the sense of [44]), one obtains the followingself-coupled delayed FitzHugh-Nagumo equation: x (cid:48) t = x t − x t + y t + J ( x t − x t − τ ) y (cid:48) t = ε ( a + b x t − γy t ) . (8)This is precisely the equation we choose to analyze in the present manuscript. When τ = 0, this equation issimply the original FhN model since the interaction term vanishes. In the limit ε = 0 the parameters ( a, b, γ )7nly affect the dynamics of the slow (reduced) equation, while the fast (layer) equation is independent of theseparameters. Among these models, we focus our attention on the simple case γ = 0, and consider withoutloss of generality b = −
1, which is equation (1). This system is particularly appealing because its equilibriahave a very simple expression as a function of the parameters. This simplicity will allow to perform analyticalcharacterization of its dynamics in the main text, and the FhN model can be understood qualitatively fromthe analysis of the fast equation. Numerical analysis of the system for γ (cid:54) = 0 is performed in appendix B. Weeventually note that the system (1), as a delayed differential equation, is classically seen as a dynamical systemin a functional space, and solutions are defined once an initial segment of trajectory is defined on a time intervalof length τ . Here, we will be interested in fixed points, periodic orbits and chaotic attractors: we will notconcentrate on transient behaviors and how the trajectories of the system depends on the choice of the initialconditions. Before undertaking the GSPT analysis of system (1), we start by investigating the fixed points of the systemand their stability for ε >
0. In the regions were the fixed point is unstable, we investigate numerically the typeof dynamics observed.
For any fixed a ∈ R , there exists a unique fixed point to the FhN system (8) given by ( x a , y a ) = ( a, − a + a / Proposition 1.
The fixed point ( a, − a + a / is • stable if | a | > √ J • unstable if | a | < • for a ∈ [1 , √ J ] , there exists one curve of Hopf bifurcations that governs changes on the stability ofthe fixed point. This curve is given by the branches τ k and τ k that together form a smooth curve alongwhich the linearized system has at least one pair of purely imaginary eigenvalues. We will refine this result in section 2.2 by showing that along the branch τ the system undergoes a Hopfbifurcation, whose sub- or super-critical type depends on the value of the parameters ( J, τ ) and changes alonga curve approaching Jτ = 1 when ε →
0. 8 roof.
The stability of the fixed point ( a, a − a /
3) is characterized by the sign of the characteristic roots ξ ofthe system. These are found as solutions of the characteristic equation obtained by linearizing the system atthe fixed point. The linearized equations are: ˙˜ x t = (1 − a + J ) ˜ x t + ˜ y t − J ˜ x t − τ ˙˜ y t = − ε ˜ x t and the characteristic equation is obtained when looking for solutions of the form (˜ x, ˜ y ) e ξt . Substituting thisinto the linearized equation and denoting by we readily obtain: ξ ˜ x = (1 − a + J (1 − e − ξτ )) ˜ x + ˜ yξ ˜ y = − ε ˜ x which has non-trivial solutions when det( ξI − M ) = 0 with M = − a + J (1 − e − ξτ ) 1 − ε i.e.: ξ ( ξ − a − J (1 − e − ξτ )) + ε = 0 . This equation corresponds exactly to equation (14) when ε →
0. Similarly as in the fast case this equation canbe solved by elementary algebra. In contrast with the fast system, ξ = 0 is never a solution when ε >
0. Thisallows to write down the equations: Je − ξτ = − εξ − ξ + (1 − a + J ) . Denoting ξ = α + i β and taking the real part of both sides of the equality, we obtain the equation: Je − ατ cos( βτ ) = − εαα + β − α + (1 − a + J ) . If α >
0, we therefore necessarily have − ≤ e − ατ cos( βτ ) ≤ − a J and there is no possible pairs of real values α > , β ∈ R satisfying this inequality for a > J . Similarly, if α <
0, taking the modulus on both parts of the equality, we can write down the inequality:
J > Je − ατ = (cid:115)(cid:18) − a + J − α ( εα + β + 1) (cid:19) + β (cid:18) εα + β − (cid:19) ≥ (1 − a + J )which is not possible for a <
1. 9he stability of the fixed point necessarily changes within the interval [1 , √ J ], and therefore character-istic roots, which are continuous with respect to the parameters of the system, cross the imaginary axis. Sincewe already noted that ξ = 0 is not a possible root, changes of stability necessarily occur when the system haspurely imaginary solutions ξ = i ζ , that hence satisfy the complex equation: Je − i ζτ = i ( εζ − ζ ) + (1 − a + J ) . Taking the modulus, one finds, as expected, that solutions exist for any | a | ∈ [1 , √ J ] (and only in thatinterval), and in that case we find two possible characteristic roots (together with their complex conjugate). Wechoose by convention the two roots corresponding to ε/ζ − ζ = A : ζ = ( − A − √ A + 4 ε ) ζ = ( − A + √ A + 4 ε ) . With these expressions of the eigenvalues, we can now solve the dispersion relationship and find the parametersfor which Hopf bifurcations occur. Taking the real part, we find that the delay can take the two following values: τ ± = ± ζ cos − (cid:18) − a J (cid:19) + 2 kπζ with ζ = ζ or ζ = ζ and k ∈ Z . The imaginary part of the dispersion relationship provides an additionalequation that is only satisfied when τ = τ − . This leads to the formulae (4) where the k have been fixed to ensurepositivity of the curves. It remains to show that these two families of curves actually form together a uniquesmooth curve in the parameter plane ( a, τ ) that tends to infinity as visible in Figure 1(b). We demonstrate thisfact by showing that the expansion of the different branches of the curve match at second order at the turningpoints a ∈ { , √ J } . At a = 1, we have: τ k ( J, ε, a ) = kπ √ ε + K k √ ε √ a − O ( a − τ k ( J, ε, a ) = k +1) π √ ε − K k +11 √ ε √ a − O ( a − K k = 2 / √ J − Jkπ/ √ ε . This implies that the union of τ k +11 and τ k , seen as a function of τ , is the graphof a map at least twice continuously differentiable at τ = k +1) π √ ε . Similarly, at a ∗ = √ J we have: τ k ( J, ε, a ) = (2 k +1) π √ ε + K k √ ε √ a − a ∗ + O ( a − a ∗ ) τ k ( J, ε, a ) = k +1) π √ ε − K k √ ε √ a − a ∗ + O ( a − a ∗ )ensuring a smooth connection between τ k and τ k . We eventually notice that τ k belongs to the interval [2 k +1 , k + 2] π/ √ ε and therefore disappear at infinity in the limit ε →
0. In contrast, τ k ≤ τ kf the k th branch of10opf bifurcations of the fast system given by equation (3) and announced in theorem 2: as ε →
0, the twofamilies become infinitely close in any compact contained in (1 , √ J ). We conclude that τ is the primaryHopf bifurcation curve for τ < π/ √ ε , i.e. for fixed a , the first value of the delay leading to a change of stabilityof the fixed point. It smoothly connects τ = 0 at a = 1 to τ = π/ √ ε at a = √ J , where it has a turningpoint and merges smoothly with τ at a = a ∗ , a branch ending at a = 1 and τ = 2 π/ √ ε , and the matchingpersists for increasing values of k : • τ k connect 2 kπ/ √ ε at a = 1 to (2 k + 1) π/ √ ε at a = √ J • τ k that connects (2 k + 1) π/ √ ε at a = √ J to (2 k + 2) π/ √ ε at a = 1 • τ k and τ k smoothly connect at a = √ J and τ k smoothly connects to τ k +11 at a = 1.The figure 1(b) moreover shows that the curve forms a sequence of loops, visible for ε large enough. Sufficientlyclose to the boundary a = √ J , increasing τ leads to a series of Hopf bifurcations changing the stability ofthe fixed point. But as ε →
0, all this complexity escape to infinity, and the curve τ essentially controls thestability of the fixed points for finite delays. In this section we compute the first Lyapunov coefficient λ of the normal form at the primary Hopf bifurcationat τ = τ ( J, ε, a ). We show the following transition:
Proposition 2.
In the limit ε → , the Hopf bifurcation is supercritical for Jτ < and subcritical otherwise.Proof. We first focus on the case when
Jτ <
1. Note that in this case the Hopf bifurcation is found very closeto the line a = 1. Consider the extended system obtained from (1) by adding the equation ε (cid:48) = 0. For thissystem there a center manifold reduction, based on the classical theory developed, for example, in [17]. We omitthe straightforward details of this computation, that lead to show that on the center manifold, the reductionfollowed by a reflection in the ˜ x variable, yields the following system on the center manifold: d˜ x d θ = − ˜ y + ˜ x − ˜ x + a ε ˜ x + O (˜ y , ˜ x ˜ y, ˜ x ) d˜ y d θ = ε (˜ x − ˜ a + O ( ε ˜ y , ε ˜ x ˜ y, ε ˜ x )) , (9)where ˜ a = 1 − a , a = Jτ − Jτ ) , (10)If the higher order terms are omitted (9) differs from the classical van der Pol system by the term a ε ˜ x . Notethat a is positive and blows up as τ approaches 1 /J . 11 SN f H f (a) ε = 0 . (b) ε = 2 Figure 1: Hopf bifurcations for system (1) as a function of ( a, τ ) with J = 2 and (left) ε = 0 .
01: the red curverepresents τ which is the only visible Hopf bifurcation in that range of parameters. The Hopf bifurcation issupercritical below the Bautin point B (solid line) and subcritical above (dotted line). Superimposed is thebifurcations of the fast system (12) with J = 2 (blue curves). SNf: saddle node, Hf: subcritical Hopf, closelyapproaching τ on Jτ <
Jτ > ε = 2 shows how the different branches of Hopfbifurcation form a single smooth curve. In the range of values plotted, we see τ k and τ k for k ∈ { · · · } .12ystem (9) has a canard point (singular Hopf point) at an associated Hopf bifurcation, which is the sameas computed in Section 2.1. A straightforward application of formula (3.4.11) in [22] yields the first Lyapunovcoefficient for this Hopf bifurcation: λ = −√ ε/ O ( ε ) . (11)Hence the bifurcation is supercritical and yields a stable limit cycle.Concerning the case Jτ > ε >
0, yet the Lyapunov coefficient of this Hopf bifurcation may not converge to the Lyapunov coefficient of thefast Hopf bifurcation as ε →
0. This discontinuity relies on the presence of quadratic terms in the slow equation.Note that for (1) no such terms are present and that when ε = 0 the slow variable is a parameter, so that thecenter manifold reduction does not introduce quadratic terms without a factor of ε (see the remainder terms in(9)). Hence it is easy to check using the approach [52] that in this case the Lyapunov coefficient is continuousin the limit ε →
0. We conclude that Hopf bifurcation for
Jτ > ε > λ directly using the normal form reduction, butthe expressions become complicated. Aided by formal calculation software (Maple R (cid:13) ), we did the calculationof the first Lyapunov coefficient along the Hopf bifurcation line τ ( J, ε, a ) following the classical theory [18, 7].Details are provided in appendix A. The very complicated expression of the Lyapunov coefficient obtained inthis manner is evaluated numerically and confirms that there exists τ s ( ε, J ) >
0, approaching 1 /J for small ε ,such that for τ < τ s ( ε, J ), the bifurcation is supercritical and for τ > τ s ( ε, J ), it is subcritical. A plot of theLyapunov coefficient for ε = 0 . J and τ appears in Figure 2 (A). The change of sign of theLyapunov coefficient occurs along a line τ s closely following τ = 1 /J . For τ < /J , the Lyapunov coefficientis negative and of small absolute value. It then sharply jumps towards higher values before slowly returningto smaller positive values. This is visible in the zoom of the Lyapunov coefficient for J = 2. This is perfectlyconsistent with the properties that: (i) the Lyapunov coefficient is of order √ ε for Jτ <
1, and (ii) the remarkthat the Lyapunov coefficient of the fast system tends to infinity at Jτ →
1. Eventually, we note that at τ = τ s ,the the first Lypaunov coefficient changes sign in a non-degenerate manner, and the map that associates to τ the eigenvalues of the Jacobian matrix and the first Lyapunov coefficient is regular. We conjecture that thesystem undergoes a Bautin bifurcation at this point. In order to rigorously demonstrate the presence of thisbifurcation in the system, one needs to compute the second Lyapunov coefficient at this point and show that it13s not zero: this amounts to extending the analysis performed in appendix A to access terms of order 5 of thereduction of the system on the center manifold, leading to tedious computations. / J J (A) (B) < 0 Figure 2: First Lyapunov coefficient of the Hopf bifurcation at τ = τ ( J, ε, a ) with ε = 0 .
01. (A) Lyapunovcoefficient as a function of J and τ . Negative values are plotted in white: the Lyapunov coefficient changes signalong a curve closely approximating τ = 1 /J . (B) We zoomed for J = 2 on the negative part of the curve ofLyapunov coefficients. We observe the sharp change occurring close to 0 .
503 and disappears to positive values(red arrow).
The presence of Hopf bifurcations marks the transition between regions of parameters where the fixed pointis stable or unstable: in particular, in the parameter region depicted in yellow in Figure 1, the fixed point isunstable, and oscillations or chaotic solutions arise. In this section we give a preview of what can happen as τ grows for fixed a > ε > J >
0. The topic of section 5 will be to use the GSPT to uncover the dynamicalorigin of these phenomena. Distinct qualitative behaviors as τ is increased, summarized in Figure 3.For values of the delay smaller than the value corresponding to the Hopf bifurcation ( τ = τ ( J, ε, a ) (cid:39) . v a , w a ) remains stable, and precisely at crossing the Hopf bifurcation line, a small cycle appears.In a very small range of delay values, the cycle can become more complex, indicating the possible presence ofperiod doubling bifurcations. A cascade of such bifurcations seem to lead to a chaotic behavior as delays arefurther increased (although in a very small neighborhood of the Hopf bifurcation), as depicted in Figure 3- SmallOscillations .As the delay is further increased a sudden switch to large oscillations takes place, a hallmark of thepresence of canard transitions in the system. These oscillations are more complex than relaxation oscillationsemerging in two-dimensional slow-fast systems. They show the typical shape of Mixed-Mode oscillations (see14igure 3-
MMOs ), with three distinct phases: one small oscillation related to the ghost of the small periodicorbit (highlighted in blue boxes in the figure, ‘classical’ MMO small oscillations), small oscillations related to theproperties of the convergence of the delayed fast equation towards the critical manifold (red boxes), and largerelaxation oscillations corresponding to switches between the two attractive parts of the critical slow manifold.We will study this phenomenon more in details in section 5.1.As the delay is further increased, these MMOs suddenly disappear in favor of bursting periodic orbits (seeFigure 3- (Bursts) ), characterized by the presence of very fast oscillations interspersed with periods of slowbehavior. The analysis of these orbits will be performed in section 5.2, of which we display time series andtrajectories in the extended phase space ( x t , y t , x t − τ ).We discuss also the behavior found at the transition from MMOs to bursts, where the system displaysextremely wild chaotic behaviors characterized by irregular switches between MMOs and bursts (Figure 3- (Chaos) ). These phenomena will be accounted for in section 5.3. In order to illustrate the chaotic nature of thetrajectory, it is generally convenient to display a Ruelle plot, as we did for the chaotic small oscillation. Buthere, Ruelle plots on Poincar´e section in the axes ( x, y ) are scattered over very complex attractors that do notseem to be one-dimensional, and do not provide a concise information of the chaotic nature of the trajectory.Rather, we found that the value of x at time t − τ at crossing a Poincar´e section in y shows a lower-dimensionalbehavior that we display in Figure 4. In order to explain the complex dynamics presented in the previous section, the singular perturbation theorytakes into the account properties of both the fast and the slow flow in the limit ε = 0. In the case of (1) theslow dynamics is very simple and the challenge is to understand the fast dynamics given by the solution ofequation (8) in the limit ε = 0: x (cid:48) t = x t − x t y + J ( x t − x t − τ ) (12)The slow variable y is now looked at as a parameter of the dynamics. An important object are families ofequilibria of (12), known also as critical manifolds, as they give rise to slow manifolds of (1). In this sectionwe identify the critical manifolds of (12), characterize their stability and find non-hyperbolic bifurcation points(Hopf) induced by the presence of delays. 15 .0510.951.051.10 50 100 -0.67 -0.66-0.65 -2 -0.67 -0.65 -2 -0.69 -0.63 -2 time -3030 2000 4000 -2
500 600 S m a ll O s c ill a t i o n s MM O s C h a o s B u r s t s time -2 time -0.50.5 -3030 2000 4000 timetime time time time timetime -2 time -3 -3 Figure 3: Trajectories of the full system (1) with a = 1 . J = 2 and ε = 0 .
05 as delays are increased. (SmallOscillations) from left to right, τ = 0 . τ = 0 .
401 and τ = 0 . τ = 0 .
41, as illustrated by the Ruelle plot on a Poincar´esection (blue line at y = − . τ = 0 .
51 (left) τ = 0 .
55 (right). (Chaos) τ = 0 .
7: time series showsirregular alternation of bursts, MMOs and mixed behaviors between MMO and burst, (Bursts) τ = 1.16 Figure 4: Sequence of values of x t n − τ where t n are the consecutive times the trajectory crosses the Poincar´esection y = − .
4. Parameters J = 2, ε = 0 .
05 and τ = 0 . The first step in analyzing the fast equation is to characterize its fixed points (critical manifold). These aregiven by the solutions to the algebraic equation: x − x y = 0 . (13)We therefore have three fixed points when | y | < and one fixed point otherwise. A closed from expression isobtained using Cardano’s method: • for | y | > /
3, the unique solution is given by: x = (cid:32) y + (cid:112) y − (cid:33) / + (cid:32) y − (cid:112) y − (cid:33) / • for | y | < /
3, ∆ <
0, the three solutions are given by x k = 2 cos (cid:18)
13 cos − (cid:18) y (cid:19) + 2 kπ (cid:19) , k = 0 , , , that we order and denote x − ( y ) < x ( y ) < x + ( y ). • Eventually, for y = ± /
3, we have a double root x = ∓ x = ± x + ( y ) ≥ y ≥ − / x − ( y ) ≤ − y ≤ /
3, and x ( y ) ∈ [ − ,
1] for y ∈ [ − / , / In the absence of delays, x ( y ) is always a saddle fixed point while x ± ( y ) are stable fixed points, and the fastequation has saddle-node bifurcations at y = ± . The presence of the delay term, which did not affect theshape of the critical manifold, may affect its stability. The stability of the fixed points on the critical manifoldwith respect to the dynamics of (12) is characterized in the following:17 /J BTBT SN SNH H y (a) Bifurcations fast system in ( y, τ ) SNBT fff (b) Bifurcations fast system in ( a, τ ) with a = x + ( y ) Figure 5: Hopf bifurcation for the fast system (12) for J = 2. (a) The branches of solutions x ± ( y ) undergosaddle-node bifurcations at y = ± / τ = τ f with verticalasymptotes at y = ± / y = ± √ J ( J −
1) (black dashed line). (b) Primary and secondary Hopfbifurcations in the plane ( a, τ ) with a = x + ( y ). The visible curves in the range plotted are τ kf for k ∈ { , , } .The BT star identify the generic subcritical Bogdanov-Takens bifurcations characterized in theorem 3. Proposition 3.
System (12) undergoes the following bifurcations: • saddle-node bifurcations at y = ± that correspond to parameters for which the fixed points x ± ( y ) collidewith x ( y ) . These bifurcations are independent of τ . • For any y ∈ [ − , ( J − √ J ] , there exists a unique value of τ given by τ f ( x ± ( y )) (recall the definitionof τ kf in equation (3) ) above which the fixed point x ± ( y ) loses stability through a generic subcritical Hopfbifurcation. Secondary Hopf bifurcations arise for τ = τ kf ( a ) for any k ∈ N \ { } .These bifurcations organize the stability of the three fixed points. For J > , their stability is as follows • For any | y | < / fixed point x ( y ) is a saddle • For | y | > ( J − √ J , the fixed points x ± ( y ) are stable regardless of the value of τ • For | y | < ( J − √ J , the fixed points x ± ( y ) are stable for τ < τ ( x ± ( y )) , and unstable otherwise. The Hopf and saddle-node bifurcation curves are depicted for J = 2 in Figure 5.18 roof. We start by characterizing the stability of the fixed points away from the bifurcation points announcedin the proposition. To this end, we investigate the sign of the real part of the characteristic roots of the systemat a given fixed point x ∗ . These are the solutions ξ of the equations: ξ = 1 − ( x ∗ ) + J − Je − ξτ (14)This equation may be solved using special functions , but elementary calculus can also be used to characterizethe bifurcation curves in the parameter space.Indeed, from the complex characteristic roots equation (14), we can show that x ( y ) is unstable whateverthe parameters ( J, τ ) of the system, and x ± ( y ) are stable when | y | > ( J − √ J . Indeed, if ξ = α + i β ,we can deduce from (14) that: Je − ατ = (cid:113) ( − α + (1 − ( x ∗ ) + J )) + β , which is not possible with α < | x ∗ | <
1, since in that case we would have the contradiction: J ≥ Je − ατ ≥ − α + (1 − ( x ∗ ) + J ) ≥ − ( x ∗ ) + J, i.e. ( x ∗ ) ≥ . Similarly, we can show that no unstable solution exist with | y | > ( J − √ J . Indeed, in that case wehave | x ∗ | > √ J . Assuming α >
0, equation (14) implies that necessarily: − J ≤ Je − ατ cos( βτ ) ≤ − ( x ∗ ) + J, i.e. ( x ∗ ) ≤ J In the parameter region not covered by the above characterization, the system may display changes ofstability. We now use the characteristic equation (14) to determine parameters corresponding to such changesof stability. We therefore only need to access to the set of parameters for which the characteristic roots of oneof the equilibria cross the imaginary axis. There are two such situations: the characteristic root crossing theimaginary axis can either be (i) a real characteristic root (the saddle-node case) or (ii) a pair of purely imaginarycharacteristic roots (Hopf case).Saddle-node bifurcations (corresponding to ξ = 0) occur if and only if: x ∗ = ± Indeed, the solutions to the characteristic equation are given by the Lambert functions W k (the different branches of the inverseof x (cid:55)→ xe x , see e.g. [9]): ξ = A + 1 τ W k (cid:16) − τJe − τA (cid:17) with A = 1 − ( x ∗ ) + J . The stability of x ∗ is hence governed by the sign of the real part of the rightmost eigenvalue, given by thereal branch W of the Lambert function, and if the argument of the Lambert function has a real part greater than − e − the rootis unique. If not, we have two eigenvalues with the same real part corresponding to k = 0 or − x ∗ as a fixed point of the system, to the point in theparameter space: y = ± τ and the coupling strength J .Hopf bifurcations occur when ξ = i ζ for ζ >
0, thus necessarily when J = (1 − ( x ∗ ) + J ) + ζ i.e. ζ = (cid:112) ( x ∗ ) − (cid:112) J + 1 − ( x ∗ ) . (15)The first observation is that Hopf bifurcations may only arise for | x ∗ | ∈ (1 , √ J ), i.e. the only fixed pointsthat may undergo a Hopf bifurcation are x ± ( y ). These fixed points lose stability when τ exceeds τ kf ( x ∗ ) givenby: τ kf ( x ∗ ) = 1 ζ (cid:18) cos − (cid:18) − ( x ∗ ) + JJ (cid:19) + 2 kπ (cid:19) . (16)In order to show that indeed, the system undergoes a generic subcritical Hopf bifurcation at this point, oneneeds to reduce the system to normal form at this point and compute the first Lyapunov coefficient. Similarlyto what we did in order to characterize the criticality of the Hopf bifurcation of the full system in section 2and appendix A, we compute a complex but exact formula for the first Lyapunov coefficient of the fast systemaided by a formal computation software (Maple R (cid:13) and Mathematica R (cid:13) ). Numerical evaluation of this coefficientshows that it is always strictly positive (see Figure 6), ensuring that the system undergoes a subcritical Hopfbifurcation.We now characterize in more detail the shape of the Hopf bifurcation curves. When y → ± √ J ( J − x + ( y ) approaches √ J and all branches of Hopf bifurcations escape to infinity, since theHopf bifurcation curves τ kf ( x ∗ ) behave, for x ∗ approaching √ J + 1 from the left, as: τ kf ( x ∗ ) ∼ (2 k + 1) π √ J (cid:112) (2 J + 1) − ( x ∗ ) . For x ∗ → + , the behavior of τ f differs from that of τ kf for k ∈ N \ { } . For k (cid:54) = 0, the Hopf bifurcation curves τ kf ( x ) diverge to infinity when x → + . In contrast, the curve of Hopf bifurcations τ f approaches the curve ofsaddle node bifurcations tangentially at this point, since using the series expansion of cos − (1 + z ), we have: τ f ( x ∗ ) = 1 J + 112 J (1 − ( x ∗ ) ) + O ((1 − ( x ∗ ) ) ) . This is also the case of the graph of the Hopf bifurcation curve parameterized by y : for instance looking at thebifurcation around x − ( y ) at y → − shows that the curve of Hopf bifurcation meets the line of saddle nodebifurcations y = 2 / τ = 1 J − J (cid:114) − y + O ( (cid:112) / − y ) . .51 20.5 -2/3 y Figure 6: First Lyapunov coefficient of the normal form of the Hopf bifurcation of the fast system. For thewhole range of parameters (
J, y ) considered, the Lyapunov coefficient is strictly positive identifying a genericsubcritical Hopf bifurcation. The Lyapunov coefficient diverges when y → ± / τ → /J and goes to zeroat the other end of the Hopf bifurcation curve ( y → ± √ J ( J −
1) and τ → ∞ ). (Bottom) Lyapunovcoefficient along the Hopf bifurcation curve for J = 2 (red curve on the top figure), as a function of y (left) or τ (right). 21his is consistent with the hypothesis that the system undergoes a generic subcritical Bogdanov-Takens (BT)bifurcation, which we demonstrate in the next section. This BT bifurcation point plays a central role in thedynamics of the full system, organizing the shape and criticality of the Hopf bifurcation of the system asdiscussed in section 2. The curve of Hopf bifurcations ends at a codimension two point corresponding to y = 2 / τ = 1 /J . Atthis point, the saddle-node bifurcation line y = 2 / y close to one, and we shall now reduce the systemto normal form around this point. Theorem 3.
System (12) undergoes generic subcritical Bogdanov-Takens bifurcation at y = ± / and τ = 1 /J .Proof. Let us start by noting that at this point the linearized operator corresponds to a the Bogdanov-Takenssingularity. Indeed, the linearized system at this point is given by the simple equations: u (cid:48) t = J ( u t − u t − τ ) . (17)It is very easy to see that the two dimensional space of affine functions t (cid:55)→ α + tβ is contained in the solutionspace of (17) with the action of the RHS of (17) given by α + βt → β . In the basis given by the functions 1and t this gives the matrix . Using the dispersion relationship (14) we can draw some conclusions on the bifurcation diagram in the plane( τ, y ). First of all the saddle node bifurcation occurs for y = − / y = 2 /
3, case treated in an identicalmanner). Near such a bifurcation point, the fixed point relationship (13) ensures that y = − / x ∗ − ) + O (( x ∗ − ), and using this expansion and the dispersion relationship (14), one can easily show that the curveof Hopf bifurcations is tangent to the y = − / z = y + 2 / u ( t ) = x ( τ t ), so that equation (12) writes˙ u = τ J ( u t − u t − ) − τ u t + τ z = δ u + νδ u − τ u + τ z δ u = ( u t − u t − ) and ν = ( Jτ − Jτ = 1. Normal form reduction for such delayed equations was analyzed in [17]. Applyingdirectly the methodology presented in that article, we obtain: ˙ x = νx + x − τ x + τ y ˙ x = 2 νx − τ x + 2 τ y. From these equations, we pursue the reduction to bring the system into canonical form of the Bogdanov-Takensbifurcation as for instance given by [33, 22], and obtain, with our notations, ˙ x = z ˙ z = 2 τ y + 2 νz − τ x − τ x z. This is the canonical normal form of the subcritical Bogdanov-Takens bifurcation.This theorem ensures in particular that beyond the Hopf and saddle-node bifurcations already identified,the system presents also a global bifurcation: a curve of saddle-homoclinic bifurcations which, near the point Jτ = 1, has the expansion { y = 2 / − / Jτ − /τ, Jτ > } . (18) Analytical methods have allowed characterizing all local codimension one bifurcation curves (saddle-node andHopf bifurcations) and a codimension two BT bifurcation. From these bifurcations emerge families of periodicorbits and homoclinic cycles that we characterize numerically in the present section. Universal unfolding ofthe BT bifurcation discussed in section 3.3 provides us with a local equivalent of the curve close to y = 2 / Jτ = 1, given by the quadratic curve (18). In order to understand the global fast dynamics away fromthe BT bifurcation, we need to continue curves of homoclinic loops, which is very complex mathematicallyand a delicate numerical task. Here, aided by numerical simulations performed using XPPAut, we obtaineda relatively simple bifurcation picture (represented schematically in Figure 7), that may be explained throughsimple, two-dimensional representations of a virtual phase plane given by ( x t , x t − τ ) (Figure 8). These phaseportraits suggest that the dynamics of (1) is essentially two-dimensional, which would suggest the existence ofan inertial manifold possibly described by the two variables x ( t ) and x ( t − τ ).The numerically evidenced bifurcations of cycles of the system, obtained through extensive simulations, arethe following: The proof of existence of such a manifold is beyond the scope of this paper. BTBT SN H Sh1
FLC Dh A B C D EFG
Sh2
Figure 7: Bifurcation diagram of the fast system in the parameter space (
J, y ). The diagram is symmetric withrespect to the reflection y (cid:55)→ − y . The unfolding of the Bogdanov-Takens points (BT) at τ = 1 /J yields aline of saddle-node bifurcations (SN, blue) at y = 2 /
3, a curve of Hopf bifurcations (H, red), together with asaddle-homoclinic bifurcation (Sh1, green) enclosing one fixed point. The fold of limit cycles (FLC) is indicatedin orange, and the second line of saddle homoclinics (Sh2) enclosing two fixed points is depicted in dashedgreen. Saddle homoclinic lines intersect along the line y = 0, giving rise to a double homoclinic loop (Dh).Seven specific locations are indicated with blue stars, with extended phase portraits provided in Figure 8. • The saddle-homoclinic bifurcation curves emerge from the BT points and end on the opposite manifoldof saddle-node bifurcations (solid green curves Sh1). These homoclinic loops enclose only one fixed point.Because of the symmetry of the system ( x, y ) ↔ ( − x, − y ), there exists a double-homoclinic loop at y = 0for a given value of τ > /J (Dh point). • We numerically establish the presence of additional periodic orbits of relatively large amplitude, one stableand enclosing the other one which is unstable. As parameters are varied, we find that the two cycles areconnected through a codimension two fold of limit cycles (orange curve in Figure 7). • The smaller and unstable limit cycle described above terminates, in a certain range of parameters, on themanifold of saddle fixed point in a saddle-homoclinic loop (dotted green line Sh2 in Figure 7). Contrastingwith the curve Sh1 of homoclinic loops, the Sh2 curve correspoonds to homoclinic loops enclosing the twofixed points. Moreover, the curve also intersects the line y = 0 at the point Dh where the system displaysthe double-homoclinic loop already mentioned above. • for the values of the parameters for which the fold of limit cycles exists but the Sh2 curve no longerexists (for values of τ larger than the value corresponding to the double homoclinic loop Dh), it actually24ontinues the family of unstable limit cycles emerging from the Hopf bifurcation.Although this description may appear relatively complex, it seems to be constrained to a planar dynamicalsystem and therefore has a fairly simple geometrical interpretation in the virtual extended phase plane ( x t , x t − τ )mentioned above. We now describe, in this space, the changes of the dynamics as a function of τ and y . -6-4-2024-6 -4 -2 0 2 4 -6-4-2024-6 -4 -2 0 2 4 (E) (F) -6-4-2024-6 -4 -2 0 2 4 (G) Figure 8: Phase space representation of the fast system in the axes ( x t − τ , x t ) for y = 0 and different values of τ , corresponding to the stars in Figure 7 (see text). The trajectories were obtained with XPP Aut. As before, J = 2. (A-D): y = 0 and distinct values of τ . (A) τ = 0 .
4, (B) τ = 0 .
6, (Dh) τ = 0 . τ = 0 .
7, (D) τ = 1. (E-G) y = 0 .
8, (E): τ = 1 .
5, (F) τ = 0 .
8, (G): τ = 0 .
5. All curves are actual solutions of the delayedfast system, except the dashed lines corresponding to the unstable trajectories, namely the unstable limit cyclesand the stable manifold of the saddle fixed point. See text for a description of the phase portraits.Let us first discuss the symmetric case corresponding to y = 0, in which case the system has the reflectionsymmetry x ↔ − x (see Figure 8). For small delay, the system has a phase plane similar to the case τ = 0:it features two stable foci, a saddle equilibrium and no cycle (Figure 8(A)). The first bifurcation occurring as25elays are increased is the fold of limit cycles (orange curve in Figure 7). At this bifurcation point appear a twocycles, one attractive and one repulsive (Figure 8 (B), orange and blue loops). In this regime, the stable manifoldof the saddle fixed point winds around the unstable cycle, while the unstable manifolds converge towards thestable fixed points. As delays are increased, the unstable cycle progressively shrinks and approaches the stableand unstable manifolds of the saddle fixed points. It eventually becomes identical to the stable and unstablemanifolds precisely when τ reaches the value corresponding to the double-homoclinic bifucation (Dh point inFigure 7). This loop is also the trace of the symmetric continuation of the saddle-homoclinic loops emergingfrom the two BT bifurcations, and at this point, stable and unstable manifolds of the saddle fixed point coincide(Figure 8 (Dh)). For slightly larger values of the delay, the system shows two unstable periodic orbits aroundthe stable fixed points, towards which the stable manifold of the saddle fixed point converges, while the unstablemanifold converges towards the large stable cycle (Figure 8 (C)). The unstable orbits cycling around the stablefixed points shrink, until τ reaches the value of the Hopf bifurcation (arising simultaneously at two fixed pointsdue to symmetry). At the Hopf point the stable fixed points become unstable foci and the unstable cyclesdisappear. The unstable manifold of the saddle fixed point keeps converging towards the stable relaxationcycle, and the stable manifold winds around the unstable foci (Figure 8 (D)). Asymmetric situations arisewhen y (cid:54) = 0. When | y | > /
3, the system has a unique fixed point, which cannot be unstable if | y | > √ J .The system can feature an additional pair of limit cycles (one stable and one unstable). A typical case for2 / < y < √ J is depicted in Figure 8(E-G). Let us finally discuss the what happens if 0 < y < /
3. Inthis case homoclinic bifurcations do not arise simultaneously at a double homoclinic loop, but sequentially. Anexample of phase portraits is provided in Figure 9 for y = 0 .
1. In this case, the first homoclinic bifurcation (Sh1)already arises before the fold of limit cycles. A similar situation as in Figure 8(A) occurs for small delays, and ahomoclinic loop arises (Figure 9(a)), giving rise to an unstable periodic orbit enclosing the smallest fixed point(Figure 9(b)). After the fold of limit cycles, a pair of stable and unstable cycles emerges (Figure 9(b)). Theunstable cycle progressively shrinks towards the saddle fixed point, until reaching the homoclinic bifurcationSh2 (Figure 9(c)). At this point, as delays are further increased, the unstable loop disappears, and the systemis left with an unstable orbit enclosing the smallest fixed point, a stable fixed point, a saddle and a stable limitcycle (Figure 9(d)). Further increasing the delays leads to crossing the saddle homoclinic line that emergesfrom the BT point at y = − /
3: a homoclinic loop encloses the the fixed point with the largest value of x (Figure 9(e)). For larger delays, the loop becomes an unstable periodic orbit enclosing the the fixed point withthe largest value of x (Figure 9(f)), and we are in a situation similar to Figure 8(C). As delays are furtherincreased, the the fixed point with the least value of x will undergo a Hopf bifurcation and lose stability asthe loop enclosing it collapses (not shown). The phase portrait is similar to Figure 9(e), except that the loop26igure 9: Phase space representation of the fast system in the axes ( x t − τ , x t ) for y = 0 . τ . The trajectories were obtained with XPP Aut. As before, J = 2. (a) τ = 0 . τ = 0 . τ = 0 . τ = 0 .
67, (e) τ = 0 . τ = 0 .
7. See text for a description.around the the fixed point with the least value of x disappeared and that fixed point lost stability. The stablemanifold of the saddle fixed point that converged towards this cycle now winds around this stable fixed point.Based on these descriptions of the fast dynamics, we are now in a position to account for the complexoscillatory patterns observed in our delayed slow-fast equation.27 Complex dynamics of the delayed FitzHugh-Nagumo system
Now that we characterized the dynamics of the fast system, we investigate the dynamics of the full system.The origin of small oscillations arising in the regime
Jτ < y = ± / Jτ < /
2, and we conjectured that it persists for
Jτ <
1. We refer to that paper for the explanations of the emergence of small oscillations for a close to 1, andtheir instability as delays are increased (Figure 3- Small Oscillations ).Here, we focus on the regimes
Jτ >
1, in which the dynamics of the fast system is nontrivially governedby the delays. The combination of the slow (quasi-static) evolution of the variable y as the fast system evolvesaccording to the above descriptions will provide an explanation for the different behaviors observed. Thedescriptions provided below are all valid for ε sufficiently small.For any J and τ , when a ≤ √ J , the system converges towards the stable fixed point ( a, a − a / a < √ J when τ < τ f ( J, a ) (these provide two boundaries, for any fixed τ , delineatingthe region where solutions are fixed points, or ‘stationary’ yellow region, in Figures 10 and 11).But if a < < a < √ J and τ > τ f ( J, a ), the full system does not have any stable equilibrium,and non-stationary behaviors therefore emerge. These can be of distinct types depending on the shape of thebifurcation diagram of the fast system as y varies. We describe these in more detail below. For τ slightly larger than 1 /J , we have observed that the fast system shows a subcritical Hopf bifurcation,associated with a family of limit cycles undergoing a saddle-homoclinic bifurcation, and no stable cycle ispresent in the fast system (as long as we are below the value corresponding to the fold of limit cycles). Wedisplay in Figure 10 the typical bifurcation diagram of the fast system in that regime as a function of y . For ε > y varies slowly in time and the trajectories of the full system can be understood as slowmotions along the bifurcation diagram of the fast system.In that regime, depending on a , the system either shows fixed-point behaviors (for a within the yellow regionof Figure 10) or MMOs (for a within the pink region of Figure 10).Let us consider a within the pink region. Since the fast system does not have any periodic orbit, for almostall initial conditions of the full system the fast variable converges rapidly towards a stable branch of the slowmanifold (blue line in Figure 10). From that point, as time goes by, the value of the slow variable startsevolving according to the slow dynamics, along the arrows plotted on the slow manifold. Assume for instance28 − − SNSNH H
StationaryStationaryMMO
SNSNH H − − − Figure 10: The MMO cycle: J = 2 and τ = 0 .
65. The bifurcation diagram of the fast system is superimposedwith a trajectory of the system in the plane ( x, y ). Blue lines are fixed points (solid: stable, dashed: unstable).SN: saddle-node, H: Hopf, red dashed lines: unstable periodic orbits, and orange cycle is the saddle-homocliniccycle. Dashed horizontal line corresponds to the value of a . (a) ε = 0 .
05, (b) ε = 0 .
01. Yellow and pink regionscorrespond to choices for a leading respectively to stationary or periodic (relaxation oscillations or MMO)behaviors.that it converges to the branch of the slow manifold corresponding to x < a . Then y slowly increases and thevariable x drifts accordingly along the slow manifold, and will necessarily reach the subcritical Hopf bifurcationpoint where the branch of the stable manifold loses stability. At this point, the fast system cannot stay inthe vicinity of this branch of the stable manifold and eventually jumps towards the upper branch of the slowmanifold. Because of the delays, the fast variable makes an overshoot, and subsequently displays fast dampedoscillations while converging towards the critical manifold (negative delayed feedback). On this branch x > a ,hence the variable y slowly decreases, crosses the Hopf bifurcation on the upper branch and jumps towards thelower branch of fixed points. This dynamics creates a periodic solution showing small oscillations interspersedwith large amplitude relaxation oscillations, the Mixed-Mode Oscillations shown in Figure 3. In the phaseportraits of the fast system pictures, the system slowly switches back and forth between the situation describedin Figure 8(G) and the symmetric situation. Note that for a (cid:54) = 0, the oscillation is asymmetric: for a > y is much faster on the lower branch than on the upper branch. Therefore, on theupper branch, the trajectory gets closer of the slow manifold than on the lower branch, which is highly visiblefor larger values of ε (left panel of Fig 10 for ε = 0 .
05) compared to cases where ε is smaller (right panel, ε = 0 . Remark 1.
The MMOs described in this section are different than canard mediated MMOs arising near folded ingularities [4, 10], for which transitions between regions of different numbers of small oscillations occur bymeans of a passage through a canard solution (similar to canard explosion). MMOs of this type require thepresence of two slow variables and a fold singularity. Interestingly, canard transitions seem to also play a rolein the BT mediated MMOs discussed in this paper. Typical trajectories pointing towards the presence of thisphenomenon can be computed in the delayed FhN system. Such a trajectory is depicted in Figure 3- (MMOs) (left side). As τ is increased, the bifurcation diagram of Figure 7 indicates that stable periodic orbits arise. These will beassociated to faster oscillations of the full system, either in the form of a fast oscillation or a burst. We focushere on the case where the fast system shows two folds of limit cycles connecting the branches of stable andunstable periodic orbits. This situation corresponds to the case where τ is large enough so that the homoclinicbifurcations have disappeared (see Figure 7). This is the case for instance for J = 2 and τ = 1, as we depicted inthe bifurcation diagram of Figure 11 computed with the Matlab package DDE Biftool [13, 14]. In this diagram,shown in in Figure 11(b) we display the minimum and maximum amplitude of the periodic orbits of the system(red), while in Figure 11(a) we replaced the amplitude of the cycle by the average value of the fast variable overa period of the cycle (green line). This quantity is relevant in order to understand the behavior of slow variablewhen ε > x w ( t ) the fast system dynamics for a fixed value of the slow variable y = w (whichis a parameter of the fast system), and let us virtually uncouple the slow and fast equation. Integrating theslow dynamics for x ( t ) = x w ( t ) yields: y ( t ) = y (0) + ε t ( a − t (cid:90) t x w ( s ) ds )which in the slow timescale θ = εt , reads: y ( θ ) = y (0) + θ (cid:16) a − εθ (cid:90) θε x w ( s ) ds (cid:17) (19)We have shown that x w ( t ) is one of two types: either a fixed point or a periodic orbit. In the limit ε → m ( w, τ ) = lim t →∞ t (cid:90) t x w ( s ) ds. In the case of a stationary behavior of the fast subsystem, this quantity is precisely equal to the value of thefixed point, and in the case of a periodic activity of period T , this quantity is equal to the average of the cycle30 LCHSN SNFLC Bursting Orbit H Fast OscBurstingStationary − − − (a) Abstract view SN HFLC − − (b) Actual trajectory, ε = 0 . a =1 . SN HFLC -4 time − − (c) Actual trajectory, ε = 0 . a = 0 SN HFLC -4 -303 800 time − − time (d) Actual trajectory, ε = 0 . a = 0 . Figure 11: FhN system, τ = 1. (a):Bifurcation diagram of the fast system as a function of y . Fixed points arein blue, average amplitude of cycles in green (plain line: stable, dashed: unstable). Regions of fast oscillations,bursting and stationary solutions are depicted for a < a > a = 1 . ε case. For legibility, we chose ε relatively large ε = 0 .
1, and this explains why the system slightly departs from bursting cycle depicted in subfigure. (a). Smaller ε more closely match the bursting cycle, but the number of oscillations in the burst dramatically increase (as1 /ε ) affecting legibility. (c) a fast oscillation case: a = 0. T (cid:82) T x w ( s ) ds . Therefore, in the slow timescale, introducing the delay has the effect of modifying the slowmanifold by adding to the stationary solution manifold a branch related to stable periodic orbits of the fastsystem.From this point of view, the slow manifold can be seen as the union of the critical manifold and the averagevalue of x along the cycles. It is depicted in Figure 11. It present folds that are associated to the saddle-nodebifurcations and to the fold of limit cycles of the fast system (see section 3). Along the branch of fixed points,31e also find changes of stability due to the presence of Hopf bifurcations.On that manifold, the slow variable y ( θ ) evolves depending on the relative position of the manifold comparedto the parameter a : y is increasing when m ( y ( θ ) , τ ) > a and decreasing otherwise. When the slow dynamicsreaches one of these extremal values of the slow attractive manifold, a fast switch will occur, depicted in thisfigure by the black vertical line with double arrows.This analysis allows characterizing the behavior of the system as a function of a , and reveals new behaviorsthat were not observed in Figure 3. If the value of a intersects a stable branch of the limit cycles manifold (lightblue region of Figure 11), the value of y will stabilize at a and the system will display fast oscillations (see 11(c)).For larger a (in absolute value), within the yellow region, the line x = a intersects the critical manifold at a stablefixed point of the fast system, and the full system will stabilize at a fixed point. In contrast, in an intermediateregion of values of a (pink region), the line x = a neither intersects the stable branch of average value of limitcycle nor a stable branch of the critical manifold. In that case, as plotted schematically in Figure 11(a), thesystem shows a relaxation-like cycle including a branch of the fast oscillations manifold: the system switchesperiodically between the periodic orbits manifold the fast system and the branch of fixed point of the slowmanifold, thus producing bursts: • when the fast system stabilizes on the branch of stable limit cycles, the system presents fast oscillationsin the fast variable, and the slow variable slowly increases (with a speed bounded below) • this persists until the slow variable y reaches the value associated to the fold of limit cycles of the fastsubsystem. Subsequently, the system jumps on to the slow manifold associated to stable fixed point ofthe fast dynamics • On that manifold, the slow variable y decreases with a lowerbounded speed, and therefore will reach theHopf bifurcation, subsequently leading the trajectory to leave the stable manifold and jump onto themanifold of stable oscillations.This scenario hence occurs periodically in time, defining a periodic orbit composed of a phase of fast oscillationsfollowed by a silent period, in other words a burst. A trajectory of the dynamical system superimposed tothe slow manifold is plotted in Figure Figure 11(b). In the phase portraits of the fast system pictures, thesystem slowly switches back and forth between the situation described in Figure 9(a) to that of Figure 8(E).Eventually, we note that for values of a close to the transition from fast spiking to bursting, complex oscillatorypatterns made of fast oscillations with modulated amplitude are found (see 11(d)) evocative of torus canardphenomena [5]. 32 .3 Chaotic transition We now investigate the origin of the chaotic orbits displayed in Figure 3 (Chaos) . These solutions, arising for arange of values of the delay between those related to MMOs and those related to the presence of stable burstingcycles, display irregular alternations of MMO-like or bursting orbits, together with mixed trajectories. Moreprecisely, the chaotic orbits observed (Figure 3 (Chaos) ) show the presence of relaxation cycles of relativelyregular period, but during the increase of the phase of the slow variable, the fast variable shows chaotic alterna-tions between at least 3 behaviors: (i) damped oscillations convergence towards the stable manifold (MMO-liketrajectory with no fast oscillation), (ii) large periodic oscillations (bursting-like trajectory) or (iii) a mixturebetween the two types of orbits.In order to understand the origin of this chaotic switching, one needs to understand the fate of the trajectoriesleaving the upper branch of the stable manifold (in the case a > x ( t ) , x ( t − τ ))) for the value of τ corresponding tothe chaotic region in Figure 3 (Chaos) . We observe indeed that in this case, the system can either switch towardsthe stable fixed point or to the stable periodic orbit: one branch the unstable manifold of the saddle fixed pointconverges on one side to the fixed point and the other branch to the periodic orbit, while both branches of thestable manifold wind around the unstable fixed point. The system is therefore extremely sensitive to the preciselocation of the trajectory as it leaves the neighborhood of the now unstable fixed point: in the two-dimensionalextended phase space metaphor of the system, the two branches of stable manifold of the saddle do separatethose trajectories converging towards the fixed point and those converging towards the periodic orbit. Theattraction bassins of the fixed point and cycles are interwoven in a very intricate fashion . We interpret thismechanism as the origin of the emergence of the irregular alternation of MMO-like and burst-like trajectories In the vicinity of the fixed point after the Hopf bifurcation, it is clear that the expansion rate in the delayed system becomesvery large because of the separation of trajectories between those converging to fixed points and to the periodic orbit
33n this parameter regime.In Figure 12, we represent different branches of the trajectories on the phase plane ( x, y ), together withthe bifurcation diagram of the fast system. In this diagram, we represented the regions of y for which the fastsystem has a topology similar to that represented in Figure 12(a), i.e. the presence of a stable periodic orbitand a stable fixed point, two unstable fixed point and no unstable periodic orbit. This region overlaps withthe Hopf bifurcation point, and makes the system very sensitive to the specific shape of the trajectory. Werepresented on this diagram an actual trajectory of the system during one period of the relaxation oscillation,in all three cases (MMO, burst and mixed trajectory). We have considered here a simple model of delay differential equation based on the FitzHugh-Nagumo neuronmodel, that provides a canonical example of canard explosion. This system shows a surprisingly rich phe-nomenology. For small delays, a classical canard explosion occurs similarly to what happens in the non-delayedFitzHugh-Nagumo equation, which is the topic of [32]. However, as delays are increased, they destabilizebranches of the critical manifold that are necessary to support canard cycles, and complex dynamics occur. Asdelays are increased, the first phenomenon observed is the emergence of instabilities of the small cycles. Asshown in [32], these are related to an increase of the expansion rate along the critical manifold. Indeed, althoughthe stability of the branches of the critical manifold remains unchanged, increasing the delay induces a decreaseof the contraction rate of the trajectories along the stable branch of the critical manifold passing below theincreasing expansion rate along the unstable manifold. This induces a sequence of bifurcations of canards cycles,and numerical evidence shows the emergence of period doubling bifurcations and chaotic small oscillations. Asdelays are further increased, stable branches of the critical manifold may lose stability. In that case, there is nomore room for canard cycles. In the case of the FhN system, the branches of the critical manifold lose stabilitythrough a subcritical Hopf bifurcation emerging from one of the folds through a Bogdanov-Takens bifurcations.This results, in particular, in the emergence of a new type of MMOs. For even larger delays, the fast systemshows the presence of a stable limit cycle, which yields fast oscillations and bursting in the full system. Chaotictrajectories showing irregular alternations of MMO- and burst-like trajectories were also found within a specificparameter region, in which stable manifolds in the extended phase plane ( x t − τ , x t ) show a very specific shape.Delay differential equations are infinite dimensional dynamical systems, and as such may have very complexdynamics. However, we conjecture that the phenomena observed in this system may be observed in a threedimensional system with two fast and one slow variables. Such a toy model would allow going deeper into a34 a) Phase Portrait SNSN H H FLCFLC Sh Sh − − − (b) Fast system bifurcations SNSN HH
FLCFLC − − − (c) MMO cycle SNSN HH
FLCFLC − − − (d) Bursting cycle SNSN HH
FLCFLC − − − (e) Mixed cycle Figure 12: Chaotic orbits. τ = 0 . ε = 0 .
05. (a) Phase plane in a chaotic situation ( y = − .
4) the stablemanifold of the saddle fixed point winds around the unstable orbit. (b) Bifurcation diagram of the system asa function of y . The saddle homoclinic bifurcations (Sh2) and the Hopf bifurcations (H) are such that thereexists a range of values of y (gray region) in which the phase portrait is of type (a): no unstable cycle actsas a separatrix between the fast periodic orbit and the fixed point. (c,d,e) represent three trajectories for ε = 0 .
05 (large for legibility). (c) MMO-like trajectory (delayed switch favors such trajectories), (d) bursting-like trajectory (favored by early switches), (e) is a mixed trajectory: switch from the critical manifold atan intermediate location, yielding loose convergence towards the lower branch of the critical manifold beforeswitching to a periodic orbit. 35umber of phenomena that we identified in our infinite-dimensional setting and that are still not completelyunderstood from a theoretical viewpoint. These include accounting for chaotic small oscillations when expansionovercomes contraction, the development of a single Hopf bifurcation curve with a Bautin point out of the mergingof singular Hopf bifurcations and Bogdanov-Takens bifurcation in the fast system. Eventually, such a finite-dimensional version of the system would allow understanding more precisely the nature of the MMOs emerginghere in the absence of two slow variables, as well as the chaotic orbits mixing MMOs- and burst-like trajectories.From the application viewpoint, the model analyzed is a simple representation of the activity of neuronalnetworks. This relative simplicity allowed us to uncover a number of interesting phenomena that are indeedobserved in collective firings of cells, such as transitions from subthreshold oscillations to MMOs and bursting,that are for instance observed in the inferior olive [3]. We expect that the phenomena observed here persistin more complex systems and with different kinds of synaptic interactions, as long as these have an effectivenegative delayed feedback and the same kind of local bifurcation as the delayed FhN system. Our study wasbased on assuming in particular γ = 0 in equation (8). Simulations of the delayed FhN model with γ > A Normal form reduction at Hopf bifurcation
In sections 2.1 and 3, we have identified parameters for which the linearized equation at a fixed point, respectivelyof the full system (1) or of its fast subsystem, have a pair of purely imaginary eigenvalues. This provides uswith the locus of possible Hopf bifurcations, and we show here aided by numerical evaluations that the systemundergoes generic Hopf bifurcations. To this end one needs to reduce the system to normal form on the centermanifold. This appendix is devoted to providing the outline of the method, in the context of the full system(slightly more complex than the fast one), and along the branch of Hopf bifurcations given by the first branch τ = τ ( J, ε, a ) for a ∈ (1 , √ J ), defined in equation (4). This curve correspond to a change of stability ofthe fixed point ( a, a − a / ε →
0, the curve collapses to a non-smooth curve made of asegment { Jτ < , a = 1 } and the curve of Hopf bifurcations of the fast system τ ( a ) given by equation (16).We further showed that for when ε →
0, the Hopf bifurcation corresponding to
Jτ < ε > R (cid:13) , before presenting the obtained coefficients and their shapeas a function of the parameters.This method is based on considering the delay differential equation (1) as a dynamical system in the Banachspace X of continuous functions from [ − τ,
0] to R endowed with the uniform norm. (cid:107) z (cid:107) = sup θ ∈ [ − τ, | z ( t ) | where | z ( t ) | is the Euclidian norm on R of the value of the function z ∈ X at time t ∈ [ − τ, z t ∈ X as the portion of solution ( x ( t ) , y ( t ) , t ∈ [ t − τ, t ]), with the definition z t ( θ ) = z ( t + θ ) , − τ ≤ θ ≤ , we can rewrite equation (1) as: ddt x t ( θ ) = ddθ x t ( θ ) − τ ≤ θ < (cid:104) y t (0) + x t (0) + J ( x t (0) − x t ( − τ )) (cid:105) − x t (0) θ = 0and similarly, ddt y t ( θ ) = ddθ y t ( θ ) − τ ≤ θ < (cid:104) ε ( a − x t (0)) (cid:105) θ = 0The terms within the brackets is the affine part of the flow, separated from the nonlinear cubic term in theequation on the first variable, which does not involve delays. Calculations are much simplified when changingvariables so that the equilibrium is at the origin. This change of variable simply amounts changing the origin,i.e. consider the equation in terms of ˜ x = x − x ∗ . The equations now read: ddt ˜ x t ( θ ) = ddθ ˜ x t ( θ ) − τ ≤ θ < (cid:104) (1 − a )˜ x t (0) + J (˜ x t (0) − ˜ x t ( − τ )) (cid:105) − a ˜ x t (0) − x t (0) θ = 0 ,ddt ˜ y t ( θ ) = ddθ ˜ y t ( θ ) − τ ≤ θ < − ε ˜ x t (0) θ = 0 . The linear operator decomposes into a part only depending on ˜ z t (0): A = − a + J − ε , x t ( − τ ): A = − J
00 0 . At Hopf bifurcation points, i.e. when the characteristic equation det ( λI − A − e − λτ A ) = 0has a complex solution λ = ± i ω , a complex right eigenvector is: v = i εω , which corresponds to a two-dimensional center eigenspace N : cos( ωθ ) sin( ωθ ) − sin( ωθ ) εω cos( ωθ ) εω and an infinite-dimensional stable eigenspace S . The corresponding center manifold is given by: M = { φ ∈ X, φ = Φ u + h ( u ) } where u = ( u , u ) t are the coordinates on the nullspace N and h ( u ) ∈ S . Classically, we project the solutionsto the delay differential equation on M and obtain: z t ( σ ) = cos( ωσ ) u ( t ) + sin( ωσ ) u ( t ) − sin( ωθ ) εω u ( t ) + cos( ωθ ) εω u ( t ) + h ( σ ) u ( t ) + h ( σ ) u ( t ) u ( t ) + h ( σ ) u ( t ) h ( σ ) u ( t ) + h ( σ ) u ( t ) u ( t ) + h ( σ ) u ( t ) + O ( (cid:107) x (cid:107) )where the h ijk characterize the Taylor expansion of the solution on the center manifold. These coefficients satisfya system of affine ODEs, in which the affine term depends on the orthogonal basis of the nullspace, denoted ψ ( θ ) (which is straightforward to calculate) arising from the projection of the original equation on N . Theseare now classical methods, introduced and used in [7, 45, 18]. In our case, the solution to the linear ordinarydifferential equation of h ijk yield relatively complex expressions, in terms of six constants C · · · C , that arethen solved in order to match boundary values (and solve the original DDE on the center manifold). One finds: h = a ω (cid:0) cos( ωθ )Ψ − sin( ωθ )Ψ (0) (cid:1) + C − C sin( ωθ ) cos( ωθ ) + C (cos( ωθ ) −
12 )38nd similar expressions for the other terms, and it is not hard to determine the constants. From these expressions,we obtain a system of ODEs describing the evolution in time of u : ˙ u = ωu + Ψ (0)( G ( u , u )) + O ( (cid:107) u (cid:107) )˙ u = ωu + Ψ (0)( G ( u , u )) + O ( (cid:107) u (cid:107) )with G a cubic polynomial, whose coefficients can be deduced from the evaluation of h ijk and C i , and the firstLyapunov coefficient, characterizing the type of the Hopf bifurcation, is a simple function of these coefficients.With the help of the formal calculations software Maple R (cid:13) , we compute all these coefficients analytically. Ithappens that, as of quadratic terms, only the coefficients of u in G and G are non zero, and cubic coefficientsinvolved in the computation of the Lyapunov coefficients enjoy a relatively simple form. These coefficients aregiven by: G ( u , u ) = − a Ψ (0) u + u Ψ (0) (cid:0) − − h (0) a (cid:1) − (0) ah (0) + · · · G ( u , u ) = − a Ψ (0) u − ah (0) u u + · · · where the dots correspond to terms that are not necessary to compute the first Lypaunov coefficient. From thisexpression, classical formulae (see e.g. [33]) yield a very complex expression for the first Lypaunov coefficientas a function of the parameters. This expression is however exact, and allows direct numerical computation ofthe Lyapunov coefficient, that we presented in section 2.2. B The FitzHugh-Nagumo system with γ (cid:54) = 0 In this appendix, we provide one numerical example of delayed self-coupled FitzHugh Nagumo system (8) with γ (cid:54) = 0 that has similar dynamics as the ones studied in the main text, i.e. for γ small enough. ˙ x = x − x + y + J ( x ( t ) − x ( t − τ ))˙ y = ε ( a − x + γy ) (20)The fixed points of the system are given by y = ( x − a ) /γ , where x is the solution of the cubic polynomial: γ + 1 γ x − x − aγ = 0 . This equation can be solved using Cardano method, exactly as we solve the cubic equation of the fast system.In detail, the Cardano discriminant of the equation writes:∆ = 27 b (9 a − b + 1 b ) . >
0, the system has a unique solution: x = (cid:18) − a b + √ ∆ (cid:19) / − (cid:18) a b + √ ∆ (cid:19) / , and when ∆ <
0, we find 3 roots: x k = 2 (cid:18) b + 1 b (cid:19) / cos (cid:32)
13 arccos (cid:32) − a √ b (cid:112) ( b + 1) (cid:33) + 2 kπ (cid:33) . The stability of the fixed points can be analyzed in a similar fashion as done before. In that case, denoting by x ∗ one of the equilibria found, the dispersion relationship corresponds to the determinant of the matrix: ξ − A − ε εγ with A = 1 − ( x ∗ ) + J (1 − e − ξτ ), and the dispersion relationship therefore reads: ξ − ( A + γε ) ξ + ε = 0 , or: e − ξτ = ξ − ξ (1 − ( x ∗ ) + J + γε ) + ε ( γ (1 − ( x ∗ ) + J ) + 1) Jξ + γεJ . This appears much more complex to solve in closed form.The fast dynamics is identical as that of the delayed FhN system, and therefore the analysis of section 3applies. This allows to account for most of the behaviors of the FhN system. The only phenomenon requiringsome more work is the presence of canard explosion and their type. Similarly to the non-delayed case, thedetails of the slow dynamics matter for the demonstration of the presence of canard explosions. Here, weprovide numerical evidences for canard explosion, and show that similar phenomena of emergence of Mixed-Mode oscillations and bursting appear, for the same values of parameters, in the parameter regime where FhNsystem has a unique fixed point. Let us for instance consider the FhN model for b = − γ = − .
3. Weset a = 0 .
85 (which leads the system to a similar situation as studied in the FhN system with a = 1 . τ > /J , we have seen that geometric analysis of the fast system accounted for most phenomena arising in theFhN system. This property persists for the FitzHugh-Nagumo system, and we therefore observe the very samebehaviors and transitions at similar delays. Conflict of Interest:
The authors declare that they have no conflict of interest and that they have complied with ethical standards.40igure 13: Self-coupled delayed FitzHugh-Nagumo system for τ > /J and ε = 0 .
05. (a) τ = 0 .
6, (b) τ = 0 . τ = 0 .
9. The system transitions from MMOs to bursting through a period of chaos, arising for similarparameter values as in the delayed FhN system. 41 eferences [1] Baer, S.M., Erneux, T.: Singular Hopf bifurcation to relaxation oscillations. SIAM Journal on AppliedMathematics (5), 721–739 (1986)[2] Benoit, E., Callot, J., Diener, F., Diener, M., et al.: Chasse au canard (premi`ere partie). CollectaneaMathematica (1), 37–76 (1981)[3] Bernardo, L., Foster, R.: Oscillatory behavior in inferior olive neurons: mechanism, modulation, cellagregates. Brain research Bulletin , 773–784 (1986)[4] Brøns, M., Krupa, M., Wechselberger, M.: Mixed mode oscillations due to the generalized canard phe-nomenon. Fields Inst. Comm. , 39–63 (2006)[5] Burke, J., Desroches, M., Barry, A.M., Kaper, T.J., Kramer, M.A., et al.: A showcase of torus canards inneuronal bursters. The Journal of Mathematical Neuroscience (3) (2012)[6] Campbell, S., Stone, E., Erneux, T.: Delay induced canards in high speed machining. Dynamical Systems (3), 373–392 (2009)[7] Campbell, S.A.: Calculating centre manifolds for delay differential equations using maple. In: DelayDifferential Equations, pp. 1–24. Springer (2009)[8] Chiba, H.: Periodic orbits and chaos in fast-slow systems with bogdanov-takens type fold points. Journalof differential equations , 112–160 (2011)[9] Corless, R., Gonnet, G., Hare, D., Jeffrey, D., Knuth, D.: On the lambertw function. Advances in Com-putational mathematics (1), 329–359 1019–7168 (1996)[10] Desroches, M., Guckenheimer, J., Krauskopf, B., Kuehn, C., Osinga, H., Wechselberger, M.: Mixed-modeoscillations with multiple time scales. SIAM Review (2), 211–288 (2012)[11] Drover, J., Rubin, J., Su, J., Ermentrout, B.: Analysis of a canard mechanism by which excitatory synapticcoupling can synchronize neurons at low firing frequencies. SIAM journal on applied mathematics (1),69–92 (2004)[12] Dumortier, F., Roussarie, R.H.: Canard cycles and center manifolds, vol. 577. American MathematicalSoc. (1996)[13] Engelborghs, K., Luzyanina, T., Roose, D.: Numerical bifurcation analysis of delay differential equationsusing dde-biftool. ACM Transactions on Mathematical Software (TOMS) (1), 1–214214] Engelborghs, K., Luzyanina, T., Samaey, G.: Dde-biftool v. 2.00: a matlab package for bifurcation analysisof delay differential equations. Technical Report TW-330, Department of Computer Science, K.U.Leuven,Leuven, Belgium (2001)[15] Ermentrout, G., Terman, D.: Mathematical foundations of neuroscience (2010)[16] Ermentrout, G.B., Wechselberger, M.: Canards, clusters, and synchronization in a weakly coupled in-terneuron model. SIAM Journal on Applied Dynamical Systems (1), 253–278 (2009)[17] Faria, T., Magalh˜aes, L.: Normal forms for retarded functional differential equations and applications tothe bogdanov-takens singularity. Journal of Differential Equations , 201–224 (1995)[18] Faria, T., Magalh˜aes, L.T.: Normal forms for retarded functional differential equations with parametersand applications to hopf bifurcation. Journal of differential equations (2), 181–200 (1995)[19] Fenichel, N.: Geometric singular perturbation theory for ordinary differential equations. J. DifferentialEquations (1), 53–98 (1979)[20] FitzHugh, R.: Mathematical models of threshold phenomena in the nerve membrane. Bulletin of Mathe-matical Biology (4), 257–278 0092–8240 (1955)[21] Friart, G., Weicker, L., Danckaert, J., Erneux, T.: Relaxation and square-wave oscillations in a semicon-ductor laser with polarization rotated optical feedback. Optics express (6), 6905–6918 (2014)[22] Guckenheimer, J., Holmes, P.J.: Nonlinear Oscillations, Dynamical Systems and Bifurcations of VectorFields, Applied mathematical sciences , vol. 42. Springer (1983)[23] Guckenheimer, J., Osinga, H.M.: The singular limit of a hopf bifurcation. Discrete and Continuous Dy-namical Systems series A (8), 2805–2823 (2012)[24] Hale, J., Lunel, S.: Introduction to functional differential equations. Springer Verlag (1993)[25] Higuera, M., Knobloch, E., Vega, J.: Dynamics of nearly inviscid faraday waves in almost circular contain-ers. Physica D: Nonlinear Phenomena (1), 83–120 (2005)[26] Hupkes, H.J., Sandstede, B.: Traveling pulse solutions for the discrete FitzHugh-nagumo system. SIAMJournal on Applied Dynamical Systems (3), 827–882 (2010)[27] Izhikevich, E.: Neural excitability, spiking, and bursting. International Journal of Bifurcation and Chaos , 1171–1266 (2000) 4328] Jones, C.K.: Geometric singular perturbation theory. In: Dynamical systems, pp. 44–118. Springer (1995)[29] Kandel, E.R., Schwartz, J.H., Jessell, T.M., et al.: Principles of neural science, vol. 4. McGraw-Hill NewYork (2000)[30] Koper, M.: Bifurcations of mixed-mode oscillations in a three-variable autonomous van der pol-duffingmodel with a cross-shaped phase diagram. Physica D: Nonlinear Phenomena (1), 72–94 (1995)[31] Kozyreff, G., Erneux, T.: Singular hopf bifurcation in a differential equation with large state-dependentdelay. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science (2162),20130,596 (2014)[32] Krupa, M., Touboul, J.: Canard explosion in delayed equations with multiple timescales. arXiv preprintarXiv:1407.7703[33] Kuznetsov, Y.: Elements of applied bifurcation theory, 1995[34] Mallet-Paret, J., Nussbaum, R.D.: Global continuation and asymptotic behaviour for periodic solutions ofa differential-delay equation. Annali di Matematica Pura ed Applicata (1), 33–128 (1986)[35] Mischler, S., Quininao, C., Touboul, J.: On a kinetic FitzHugh-nagumo model of neuronal network. hal-01108872 (2015)[36] Nagumo, J., Arimoto, S., Yoshizawa, S.: An active pulse transmission line simulating nerve axon. Pro-ceedings of the IRE (10), 2061–2070 (1962)[37] Plant, R.E.: A FitzHugh differential-difference equation modeling recurrent neural feedback. SIAM Journalon Applied Mathematics (1), 150–162 (1981)[38] Van der Pol, B.: On relaxation-oscillations. The London, Edinburgh, and Dublin Philosophical Magazineand Journal of Science (11), 978–992 (1926)[39] Rinzel, J.: A formal classification of bursting mechanisms in excitable systems. In: Mathematical topics inpopulation biology, morphogenesis and neurosciences, pp. 267–281. Springer (1987)[40] Rinzel, J., Ermentrout, G.B.: Analysis of neural excitability and oscillations. In: Methods of NeuronalModeling, pp. 251–292. MIT Press (1998)[41] Rubin, J., Wechselberger, M.: Giant squid-hidden canard: the 3D geometry of the Hodgkin-Huxley model.Biological Cybernetics (1), 5–32 (2007) 4442] Sherman, A., Rinzel, J., Keizer, J.: Emergence of organized bursting in clusters of pancreatic beta-cells bychannel sharing. Biophysical journal (3), 411–425 (1988)[43] Simpson, D.J., Kuske, R.: Mixed-mode oscillations in a stochastic, piecewise-linear system. Physica D:Nonlinear Phenomena (14), 1189–1198 (2011)[44] Stewart, I., Golubitsky, M., Pivato, M.: Symmetry groupoids and patterns of synchrony in coupled cellnetworks. SIAM Journal on Applied Dynamical Systems (4), 609–646 (2003)[45] Stone, E., Campbell, S.A.: Stability and bifurcation analysis of a nonlinear dde model for drilling. Journalof Nonlinear Science (1), 27–57 (2004)[46] Terman, D.: Chaotic spikes arising from a model of bursting in excitable membranes. SIAM Journal onApplied Mathematics (5), 1418–1450 (1991)[47] Touboul, J.: Limits and dynamics of stochastic neuronal networks with random delays. J. Stat. Phys 149(4), 569–597 (2012)[48] Touboul, J., Hermann, G., Faugeras, O.: Noise-induced behaviors in neural mean field dynamics. SIAM J.on Dynamical Systems (1), 49–81 (2012)[49] Touboul, J., Krupa, M., Desroches, M.: Noise-induced canard and mixed-mode oscillations in large stochas-tic networks with multiple timescales. arXiv preprint arXiv:1302.7159 (2013)[50] Weicker, L., Erneux, T., D’Huys, O., Danckaert, J., Jacquot, M., Chembo, Y., Larger, L.: Slow–fastdynamics of a time-delayed electro-optic oscillator. Philosophical Transactions of the Royal Society A:Mathematical, Physical and Engineering Sciences (1999), 20120,459 (2013)[51] Weicker, L., Erneux, T., D’Huys, O., Danckaert, J., Jacquot, M., Chembo, Y., Larger, L.: Stronglyasymmetric square waves in a time-delayed system. Physical Review E (5), 055,201 (2012)[52] Zhang, W., Kirk, V., Sneyd, J., Wechselberger, M.: Changes in the criticality of hopf bifurcations dueto certain model reduction techniques in systems with multiple timescales. The Journal of MathematicalNeuroscience (JMN)1