Bifurcations of mixed-mode oscillations in three-timescale systems: an extended prototypical example
Panagiotis Kaklamanos, Nikola Popović, Kristian Uldall Kristiansen
BBifurcations of mixed-mode oscillations in three-timescale systems:an extended prototypical example
P. Kaklamanos, N. Popovi´c, and K. U. Kristiansen ∗ Abstract
We introduce a multi-parameter three-dimensional system of ordinary differential equationsthat exhibits dynamics on three distinct timescales. Our system is an extension of both aprototypical example introduced by Krupa et al. [14] and a canonical form suggested by Letsonet al. [18]; in the three-timescale context, it admits a one-dimensional S -shaped supercriticalmanifold that is embedded into a two-dimensional S -shaped critical manifold in a symmetricfashion. We apply geometric singular perturbation theory to explore the dependence of thegeometry of our system on its parameters. Then, we study the implications of that geometry forboth the local and global dynamics. Our focus is on mixed-mode oscillations (MMOs) and theirbifurcations; in particular, we uncover a geometric mechanism that encodes the transition fromMMOs with single epochs of small-amplitude oscillations (SAOs) to those with double epochsof SAOs, and we show that the latter are more robust than in the two-timescale context. Wedemonstrate our results for the Koper model from chemical kinetics [13], which represents oneparticular realisation of our prototypical system. Finally, we illustrate how some of our resultscan be extended to more general systems with similar geometric properties, such as to a three-dimensional reduction of the Hodgkin-Huxley equations derived by Rubin and Wechselberger[20]. We consider the family of three-dimensional singularly perturbed systems of the form ε ˙ x = − y + f x + f x =: f ( x, y, z ) , (1a)˙ y = αx + βy − z =: g ( x, y, z ) , (1b)˙ z = δ ( µ + φ ( x, y, z )) =: δh ( x, y, z ) , (1c)with f , ε , and δ positive and f negative; moreover, φ : R → R is a smooth function in ( x, y, z )that will be specified further in the following. When ε and δ are sufficiently small, Equation (1)exhibits dynamics on three distinct timescales; the variables x , y , and z are then called the fast , intermediate , and slow variables , respectively. Correspondingly, Equations (1a), (1b), and (1c) arecalled the fast , intermediate , and slow equations , respectively.Mixed-mode oscillations (MMOs) are trajectories that are characterised by the alternation ofsmall-amplitude oscillations (SAOs) and large-amplitude excursions (LAOs) in the correspondingtime series. MMO dynamics is frequently observed in singularly perturbed slow-fast systems of ∗ P. Kaklamanos and N. Popovi´c: School of Mathematics, University of Edinburgh, James Clerk Maxwell Building,King’s Buildings, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom ( [email protected] and [email protected] ); K. U. Kristiansen: Department of Applied Mathematics and Computer Science,Technical University of Denmark, Asmussens All´e, Building 303B, 2800 Kgs. Lyngby, Denmark ( [email protected] ) a r X i v : . [ m a t h . D S ] S e p (a) SAOs “above” ( k = − . λ = − . k = − . λ = 2 . k = − . λ = 0 . k = − . λ = 0 . Figure 1: Oscillatory dynamics in the Koper model, Equation (2), for different values of the param-eters k and λ . (a) MMO trajectory with single epochs of SAOs and Farey sequence 2 s s s · · · ; (b)MMO trajectory with single epochs of SAOs and Farey sequence 2 s s s · · · ; (c) MMO trajectorywith double epochs of SAOs and Farey sequence 1 s s s s · · · ; (d) relaxation oscillation.the type of our prototypical example, Equation (1); one commonly accepted mechanism for theemergence of MMOs in that context is the “generalised canard mechanism”, which refers to acombination of local, dynamical passage through a so-called canard point and a suitably definedglobal return. ( Canards are trajectories that connect attracting portions of slow manifolds withrepelling ones [15, 21]; a relatively recent, exhaustive review of canard-induced MMOs can befound in [6].) One particular realisation of (1), modulo an affine transformation, is given by thethree-timescale Koper model from chemical kinetics [13], ε ˙ x = ky + 3 x − x − λ, (2a)˙ y = x − y + z, (2b)˙ z = δ ( y − z ) , (2c)with parameters k, λ ∈ R and ε and δ sufficiently small. Representative MMO trajectories thatare realised in Equation (2) can be seen in Figure 1; each such trajectory can be associated witha sequence of the form { F F . . . } , called the Farey sequence , that describes the succession of largeexcursions and small oscillations, where the segments F j are of the form F j = (cid:40) L s if the segment consists of L LAOs, followed by s SAOs “above”; L s if the segment consists of L LAOs, followed by s SAOs “below” . If a Farey sequence consists of L s -type or L s -type segments only, we say that the correspondingMMO trajectory contains single epochs of SAOs, as seen in panels (a) and (b) of Figure 1, respec-tively; Farey sequences that consist of both L s -type and L s -type segments correspond to MMOtrajectories that contain double epochs of SAOs, as shown in Figure 1(c). Finally, relaxation os-cillation refers to oscillatory trajectories that contain large excursions and no SAO segments, i.e.,trajectories with associated Farey sequence { L } ; cf. Figure 1(d).MMO dynamics has been widely studied in slow-fast systems with two distinct timescales; aparticularly fruitful approach, which is based on dynamical systems theory, combines Fenichel’sgeometric singular perturbation theory (GSPT) [10] with the desingularisation technique known as“blow-up” [15]. Of particular relevance to that approach are localised, non-hyperbolic singularitieson the corresponding critical manifolds which generate SAOs in the resulting MMO trajectories,whereas LAOs arise via a global return mechanism along normally hyperbolic portions of thosemanifolds; the reader is again referred to [6] for details and references.The geometric theory of MMO dynamics in singular perturbation problems with more than twoscales is less well-developed. In [14], a prototypical example of a “slow passage through a canardexplosion” is considered, with ε ˙ x = − y + f x + f x , (3a)˙ y = x − z, (3b)˙ z = ε ( µ + φ ( x, y, z )) , (3c)which refers to the special case with β = 0 and δ = ε in Equation (1). Asymptotic formulae forthe return map induced by the flow of (3) are derived; in particular, the underlying near-integrablestructure allows for a both qualitative and quantitative description of the corresponding MMOdynamics which includes predictions on the associated bifurcation (Farey) sequences and estimatesfor the relevant parameter regimes. In a follow-up article [4], a two-parameter modification ofEquation (3c) is considered, with ε replaced by an independent singular perturbation parameter δ in (3c); crucially, it is shown that the interplay between ε and δ can result in MMO trajectoriesin which SAO epochs of canard type alternate with those of delayed-Hopf type. In [18], the localcanonical form ε ˙ x = y + x , (4a)˙ y = − α x + βy + z, (4b)˙ z = δ (4c)is studied, which is obtained for f = 0 and µ + φ ( x, y, z ) = 1 in (1); the emergence of oscillatorydynamics is related to the interaction of a Hopf bifurcation in the corresponding fast subsystemand a folded node singularity via the resulting “canard-delayed-Hopf singularity”. Finally, variousrealisations of the general three-timescale system ε ˙ x = − y + a + a x + a x + a x , (5a)˙ y = − z + b + b x + b x + b x + cy, (5b)˙ z = δ ( µ + φ ( x, y, z )) , (5c)have been considered in [19, 7].Our principal aim in this article is a classification of the MMO dynamics in Equation (1), aswell as in the realisation thereof that is provided by the Koper model, Equation (2), in the three-timescale scenario where ε and δ are sufficiently small, on the basis of Fenichel’s geometric singularperturbation theory (GSPT) [10]. While the underlying local and global mechanisms that cangenerate SAOs and LAOs, respectively, in singularly perturbed systems of the type in Equation (1)are known, the combination of those in the present context is novel; our analysis is thereforedeliberately qualitative in nature. Our focus here is firmly on local and global bifurcations of theresulting MMO trajectories which encode transitions between the corresponding Farey sequences,as illustrated in Figure 1: we construct families of singular cycles for Equation (1) in the doublesingular limit of ε = 0 = δ ; then, we study the persistence of these families for ε and δ sufficientlysmall, and we show how the MMO dynamics of (1) can be classified in dependence of the underlyingsingular geometries – denoted as “remote”, “aligned”, or “connected”. We showcase our results intwo examples, the Koper model from chemical kinetics and the Hodgkin-Huxley equations frommathematical neuroscience. In particular, the Koper model, Equation (2), has typically beentreated as a two-timescale system with one fast and two slow variables, with δ = 1 in (2c); while thecorresponding MMO trajectories are induced by a folded node, and are highly regular [6, 11, 17],our three-timescale analysis near folded saddle-nodes of type II [21] in (1) uncovers rich MMOdynamics which is not captured by the conventional two-scale approach. One of our main resultsis the classification of that dynamics for Equation (2), as summarised in the bifurcation diagramin Figure 2: we identify subregions in the ( k, λ )-parameter plane which correspond to the variousFarey sequences illustrated in Figure 1. Moreover, we derive asymptotic approximations for theboundaries between those subregions, which are marked by transitions between the above singulargeometries, as well as by bifurcations of the associated fast subsystems; see Section 4 for details.Of particular interest here are MMOs that exhibit double epochs of SAOs, as shown in panel (c)of Figure 1; while such MMOs were reported in [6, 2] in the context of two-timescale systems, theywere found to be highly delicate there [6, Figures 16 and 22], whereas they are relatively robustin the three time-scale Koper model. Our analysis provides a clear geometric explanation for thetransition from MMOs with single epochs of SAOs to those with double epochs, as well as for therobustness of the latter.As will become apparent through our analysis, Equation (1) exhibits properties of both (3) and(4); however, it captures a plethora of phenomena that are not captured by either of those systems,as discussed in detail in Section 3. Specifically, due to the absence of a linear term in y in (3b), theintermediate dynamics therein is regular, which implies that the so-called supercritical manifoldconsidered in [14, 4] admits no degeneracies; the fact that no cubic x -dependence is present in (4a),on the other hand, diminishes the global applicability of results in [18] due to the lack of a returnmechanism. Finally, the realisations of (5) considered in [19, 7] exhibit singular geometries thatdiffer from the one considered here due to the potential for interaction between singularities on thecritical and supercritical manifolds.The article is organised as follows. In Section 2, we describe the geometry of the three-time-scale Equation (1) in the double singular limit of ε = 0 = δ : we define critical and supercriticalmanifolds; then, we construct families of singular cycles which form the basis for MMO trajectoriesof Equation (1). In Section 3, we study the singularly perturbed system in (1) for ε and δ suffi-ciently small; we classify the MMO dynamics of (1), as illustrated in Figure 1, by establishing acorrespondence with the cycles constructed in Section 2. In Section 4, we apply our results to theKoper model from chemical kinetics, Equation (2), and we elucidate in detail the structure of theFigure 2: Two-parameter bifurcation diagram of the three-timescale Koper model, Equation (2),for ε and δ sufficiently small.two-parameter bifurcation diagram in Figure 2. In Section 5, we indicate how our analysis can beextended to to a three-dimensional reduction of the Hodgkin-Huxley equations derived by Rubinand Wechselberger [20] which generalises our extended prototypical example, Equation (1). Finally,we conclude in Section 6 with a discussion, and an outlook to future research. In this section, we study the double singular limit of ε = 0 = δ in Equation (1). To that end, wefirst describe the singular geometry for ε = 0; then, we consider the resulting flow in the limit of δ →
0. Finally, we construct singular cycles which will form the basis of MMO trajectories forEquation (1) when ε and δ are sufficiently small, as considered in Section 3 below. M For ε sufficiently small and δ = O (1) fixed, Equation (1) is singularly perturbed with respect tothe small parameter ε ; in particular, (1) describes the dynamics in terms of the intermediate time t . Rewriting the governing equations in the fast time τ = t/ε , we have x (cid:48) = − y + f x + f x , (6a) y (cid:48) = ε ( αx + βy − z ) , (6b) z (cid:48) = εδ ( µ + φ ( x, y, z )) , (6c)which is a two-timescale system with one fast variable x and two slow variables y and z . Thereduced problem of the above is obtained by setting ε = 0 in (1),0 = − y + f x + f x , (7a)˙ y = αx + βy − z, (7b)˙ z = δ ( µ + φ ( x, y, z )) , (7c)while the layer problem is found for ε = 0 in (6): x (cid:48) = − y + f x + f x , (8a) y (cid:48) = 0 , (8b) z (cid:48) = 0 . (8c)We will refer to the flow that is induced by the one-dimensional vector field in Equation (8) as the fast flow ; the corresponding trajectories will be denoted as the fast fibres . The critical manifold M for (1) is a set of equilibria for (8), and is given by M := (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) f ( x, y, z ) = 0 (cid:9) = (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) y = F ( x ) (cid:9) , (9)where we define F ( x ) = f x + f x . (10)The manifold M can be written as M = S a ∪ S r ∪ F M , where S a = (cid:26) ( x, y, z ) ∈ S (cid:12)(cid:12)(cid:12) ∂f∂x ( x, y, z ) < (cid:27) and S r = (cid:26) ( x, y, z ) ∈ S (cid:12)(cid:12)(cid:12) ∂f∂x ( x, y, z ) > (cid:27) are normally attracting and normally repelling, respectively, whereas F M is degenerate due to aloss of normal hyperbolicity: F M := (cid:26) ( x, y, z ) ∈ M (cid:12)(cid:12)(cid:12) ∂f∂x ( x, y, z ) = 0 (cid:27) = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) f x + 3 f x = 0 (cid:9) . (11)In particular, we may write F M = L − ∪ L + , where L − = (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) x = 0 = y (cid:9) and L + = (cid:26) ( x, y, z ) ∈ R (cid:12)(cid:12)(cid:12) x = − f f , y = 427 f f (cid:27) ; (12)hence, it follows that S a = S a − ∪ S a + , with S a − = { ( x, y, z ) ∈ S | x < } and S a + = (cid:26) ( x, y, z ) ∈ S (cid:12)(cid:12)(cid:12) x > − f f (cid:27) , while S r = (cid:26) ( x, y, z ) ∈ S (cid:12)(cid:12)(cid:12) < x < − f f (cid:27) . The set S therefore consists of a repelling middle sheet S r and two attracting sheets S a ∓ thatmeet S r along L ± , respectively; see Figure 3. From the above, it is apparent that L − alwayscoincides with the z -axis, whereas variation in f and f translates L + , therefore “stretching” or“compressing” M . (Clearly, variation in α , β and µ has no effect on the geometry of M .) Finally,the elements of the sets Q ∓ defined by Q ∓ = (cid:8) ( x, y, z ) ∈ L ∓ | f ( x, y, z ) = 0 = g ( x, y, z ) (cid:9) are called the folded singularities of M on L ∓ , respectively [21]; for (1), these sets are the singletons Q − = { q − } and Q + = { q + } , with q − : (0 , ,
0) at the origin and q + located at x + q = − f f , y + q = 4 f f , and z + q = 4 βf f − αf f . (13) (a) (b) Figure 3: (a) The critical manifold M as the set of equilibria for the fast flow of (8); the fast fibresare parallel to the x -direction. (b) The supercritical manifold M as the set of equilibria for theintermediate flow of (19); the intermediate fibres are confined to M and evolve on planes with z = const.Finally, we consider the reduced problem on M , as given by (7), with δ sufficiently small;Equation (7) is then singularly perturbed with respect to the small parameter δ , written in theintermediate time t . To classify the folded singularities q ∓ of M , we project the flow of (7) onto M [21]: using the algebraic representation of M in (9), we can apply the chain rule to find y (cid:48) = F (cid:48) ( x ) x (cid:48) = (cid:0) f x + 3 f x (cid:1) x (cid:48) ;from (7), we therefore obtain F (cid:48) ( x ) x (cid:48) = αx + βF ( x ) − z, (14a) z (cid:48) = δ ( µ + φ ( x, F ( x ) , z )) (14b)or x (cid:48) = αx + βF ( x ) − z, (15a) z (cid:48) = δF (cid:48) ( x ) ( µ + φ ( x, F ( x ) , z )) (15b)after a rescaling of time with a factor of δF (cid:48) ( x ), which reverses the direction of the flow on S r .The folded singularities of Equation (1) then correspond to equilibria of (15); specifically, for δ > q ∓ are folded nodes [21, 18]. Their strong and weak stablemanifolds define “funnel regions” on the corresponding sheets S a ∓ , which essentially determine thebasins of attraction to q ∓ on S a ∓ . Here and in the following, we focus on the flow of Equation (1)in the vicinity of the fold line L − ; with regard to the strong stable manifold of the folded node q − ,we hence have the following result: Lemma 1.
Let G ( x , x ; z ; µ ) = (cid:90) x x F (cid:48) ( σ ) ( µ + φ ( σ, F ( σ ) , z )) ασ + βF ( σ ) − z d σ, (16) where F is defined as in (10). Then, for δ sufficiently small, the strong stable manifold of the originin Equation (15) can be written as the graph z = δ G (0 , x ; 0 , µ ) + O ( δ ) for x ∈ I, (17) where I is an appropriately defined, fixed interval about x = 0 . Proof.
Given a trajectory of (15) with initial condition ( x , y , z ) on S a ∓ , i.e., with y = F ( x ),let s denote the displacement in the x -direction of that trajectory under the corresponding flow.Then, in a first approximation, the displacement in the z -direction is given by δ G ( x , x + s ; z ; µ ),where G is defined as in (16); see [14] for details. The result is obtained by setting x = 0 = z in the resulting expression, which corresponds to the unique trajectory of (15) that passes throughthe origin.An analogous representation can be obtained for the strong stable manifold of the folded node q + . From the above, we conclude in particular that the funnels of the folded singularities q ∓ are“stretched” as δ decreases. In the limit of δ = 0, q ∓ are folded saddle-nodes of type II ; see again[21, 18] for details. For future reference, we note that the associated strong manifolds (“strongcanards”) correspond to the unique intermediate fibres on S a ∓ that cross q ∓ , respectively, while thecorresponding weak manifolds (“weak canards”) can be locally approximated by the supercriticalmanifold M which is introduced in the following subsection. M We can view the differential-algebraic systems in (7) and (15) as slow-fast vector-fields on M . Thelayer problem corresponding to (7) therefore reads0 = − y + f x + f x , (18a)˙ y = αx + βy − z, (18b)˙ z = 0 (18c)or 0 = − y + F ( x ) , (19a) F (cid:48) ( x ) x (cid:48) = αx + βF ( x ) − z, (19b)˙ z = 0; (19c)we will refer to the above as the intermediate flow , and to the corresponding trajectories as the intermediate fibres ; see panel (b) of Figure 3.Rewriting Equation (1) in the slow time s = δt , we have εδx (cid:48) = − y + f x + f x , (20a) δy (cid:48) = ( αx + βy − z ) , (20b) z (cid:48) = δ ( µ + φ ( x, y, z )) ; (20c)the reduced system that is obtained from (20) is given by0 = − y + F ( x ) , (21a)0 = αx + βF ( x ) − z, (21b) z (cid:48) = µ + φ ( x, F ( x ) , z ) , (21c)which we will refer to as the slow flow of Equation (1). The supercritical manifold M is the setof equilibria for (19), and is given by M : = (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) f ( x, y, z ) = 0 = g ( x, y, z ) (cid:9) = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) z = G ( x ) (cid:9) , (22)where we define G ( x ) = αx + βF ( x ) . (23)The manifold M can be written as the union M = Z ∪ F M , where Z = (cid:26) ( x, y, z ) ∈ M (cid:12)(cid:12)(cid:12) dgdx ( x, F ( x ) , G ( x )) (cid:54) = 0 (cid:27) (24)is normally hyperbolic and the set F M := (cid:26) ( x, y, z ) ∈ M (cid:12)(cid:12)(cid:12) dgdx ( x, F ( x ) , G ( x )) = 0 (cid:27) = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) α + 2 βf x + 3 βf x = 0 (cid:9) (25)is degenerate. Equation (25) implies F M = { p − , p + } , with p ∓ = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) x = x ∓ p (cid:9) , (26)where x ∓ p = − βf ± (cid:112) β f − αβf βf , y ∓ p = F (cid:0) x ∓ p (cid:1) , and z ∓ p = G (cid:0) x ∓ p (cid:1) . (27)The points p ∓ are called the fold points of M . Equation (27) immediately implies Proposition 1.
The manifold M admits1. exactly two fold points if and only if β f − αβf > ;2. exactly one fold point if and only if β f − αβf = 0 ; and3. no fold points if and only if β f − αβf < . Remark 1.
Under the conditions stated in Proposition 1, the fold points p ∓ of M are “inherited”from the fold lines L ∓ of M , in the sense that G ( x ) in (22) is a cubic polynomial because F ( x ) in(9) is. We note that a necessary (but not sufficient) condition for M to have two fold points is therequirement that β (cid:54) = 0. If M admits two fold points, then the normally hyperbolic portion Z of M consists of three branches: Z = Z − ∪ Z ∪ Z + , where Z − = (cid:8) ( x, y, z ) ∈ Z | x < x − (cid:9) , Z + = (cid:8) ( x, y, z ) ∈ Z | x > x + (cid:9) , and Z = (cid:8) ( x, y, z ) ∈ Z | x − < x < x + (cid:9) . (28) Proposition 2.
Assume that M admits two fold points, i.e., that β f − αβf > , by Proposi-tion 1. If β < ( β > ), then the middle branch Z of M in (28) is repelling (attracting) underthe flow of Equation (19), while the outer branches Z ∓ of M are attracting (repelling).Proof. The statement follows directly from (25) and (28).If β <
0, we may hence write Z = Z a − ∪ Z r ∪ Z a + , whereas for β >
0, we may write Z = Z r − ∪ Z a ∪ Z r + . We emphasise that Proposition 2 refers to the flow on M before desin-gularisation, cf. Equations (14) and (19), whereas the direction of the flow is reversed on S r afterdesingularisation; see (15) and Figure 3. Remark 2.
By the above, the folded singularities q ∓ of M are located at the intersections between M and L ∓ . In the three-timescale limit of ε = 0 = δ , these singularities coincide with the foldedsingularities of M in the two-timescale limit ( ε = 0 , δ > ) in our case, which stems from the factthat the fast and intermediate Equations (1a) and (1b) do not depend on δ . In this subsection, we describe the position of the folded singularities q ∓ of M relative to eachother, as well as of the fold points p ∓ of M – assuming that a pair of such points exists – relativeto the fold lines L ∓ . Proposition 3.
Assume that M admits two fold points, i.e., that β f − αβf > , by Proposi-tion 1.1. If αβ < , then both fold points of M lie on S r ;2. if αβ > , then one fold point of M lies on S a − , while the other fold point lies on S a + ;3. if α = 0 , then one fold point of M lies on L − , while the other fold point lies on L + .Proof. The result follows from a comparison of the values of x ∓ in (27) with the x -coordinates of L ∓ in the three cases where αβ < αβ >
0, and α = 0, respectively.The statements of Proposition 2 and Proposition 3 are summarised in Figure 4. We remarkthat the symmetry described in Proposition 3 breaks down when O ( x )-terms are included in theintermediate Equation (1b); see [7] for an example. If β = 0, then the projection of the criticalmanifold M onto the ( x, z )-plane is a straight line; that case has been studied in [14, 3, 4]. (a) β < α <
0. (b) β < α = 0. (c) β < α > β > α <
0. (e) β > α = 0. (f) β > α > Figure 4: Projection of the supercritical manifold M and of the fold lines L ∓ of the criticalmanifold M onto the ( x, z )-plane: in dependence on the parameters α and β , the pair of foldpoints p ∓ of M lies either on S r (panels (c) and (d)), on S a ∓ (panels (a) and (f)), or on L ∓ (panels (b) and (e)).We now turn our attention to the location of the folded singularities of M relative to eachother and with respect to the fast and intermediate fibres defined previously; recall Figure 3. We1first define planes that contain the folded singularities and that are perpendicular to the fold lines L − and L + , as follows. Definition 1.
Denote by P ∓ the planes P ∓ = (cid:8) ( x, y, z ) ∈ R | z = z ∓ (cid:9) , where z ∓ are the z -coordinates of the folded singularities q ∓ of M on L ∓ , respectively. We will refer to P ∓ as normalplanes in the following. Definition 2.
The folded singularities q ∓ of M are said to be1. aligned if P − ≡ P + ;2. connected if they are not aligned and if P ∓ ∩ Z ± (cid:54) = ∅ ; and3. remote if they are neither aligned nor connected, i.e., if P − (cid:54)≡ P + and P ∓ ∩ Z ± = ∅ . In dependence of the parameters α , β , f , and f in Equation (1), we have the following resulton the position of q − and q + relative to each other: Proposition 4. empty text, merely to move following to new line1. For αβ < , the folded singularities q ∓ of M are aligned if αβ = f f , connected if αβ > f f ,and remote if αβ < f f .2. For αβ ≥ with β (cid:54) = 0 , the folded singularities q ∓ are connected.3. For β = 0 with α (cid:54) = 0 , the folded singularities q ∓ are remote.Proof. The statements follow from Equation (13) and the properties of G ( x ) in (23); see panels (c)and (d) of Figure 5 for cases corresponding to the first statement, and panels (a), (b), (e), and (f)for cases corresponding to the second statement.In what is to come, we will restrict our attention to the case that is illustrated in panel (c) ofFigure 4: Assumption 1.
In the following, we assume that α > and β < in Equation (1). Assumption 1 is made for three reasons. First, it is consistent with the Koper model, Equa-tion (2), after transformation to the prototypical Equation (1). (In particular, it follows that thescenarios illustrated in panels (b) and (e) of Figure 4 cannot be realised in (2).) Second, remotesingularities can only be present when αβ <
0. Third, given Assumption 1, the outer branches of M are attracting, while the middle branch is repelling, which allows for the construction of closedsingular periodic orbits (“cycles”), as will become apparent in the following subsection. We now consider the reduced flow on M . We impose the following assumption on the function φ ( x, y, z ) in the slow Equation (1c): Assumption 2.
The function φ ( x, y, z ) in Equation (1c) is such that φ ( x − q , F ( x − q ) , G ( x − q )) = 0 , φ ( x, F ( x ) , G ( x )) > for x < x − q , φ ( x + q , F ( x + q ) , G ( x + q )) ≤ , and φ ( x, F ( x ) , G ( x )) < for x > x + q . (a) Remote singularities. (b) Aligned singularities. (c) Connected singularities. Figure 5: Relative geometry of the folded singularities q ∓ of M according to Definition 2 (toprow); bifurcation of the resulting singular cycles, as described in Proposition 5 (bottom row).The properties of the reduced flow on the portion Z a ∓ of M therefore depend on µ ; in partic-ular, we have that for µ = 0, a true global equilibrium of the system coincides with q − ; see Section3. Assumption 1 and Assumption 2 together imply the existence of singular cycles in Equation (1),the properties of which depend on the relative position of the folded singularities q ∓ of M , asclassified in Proposition 4. Here, these cycles are defined as the concatenation of singular orbits forthe corresponding limiting systems in (8), (19), and (21), respectively. Proposition 5.
Assume that Assumption 1 and Assumption 2 hold.1. If the folded singularities q ∓ of M are remote, then there exist a singular cycle evolving on P − , a singular cycle evolving on P + , and a family of singular cycles in between; each of thecycles in that family evolves on a plane parallel to P ∓ that lies between P − and P + . Thesecycles are “two-scale”, in the sense that the singular dynamics on them alternates between thefast timescale and the intermediate timescale (on M \M ).2. If q ∓ are aligned, then there there exists exactly one singular cycle that evolves on the plane P := P − ≡ P + . This cycle is “two-scale”, in the sense that the singular dynamics on italternates between the fast timescale and the intermediate timescale (on M \M ).3. If q ∓ are connected, then there exists exactly one singular cycle that evolves on a subset of P − ∪ P + ∪ Z a ∓ . This cycle is “three-scale”, in the sense that the singular dynamics on italternates between the fast timescale, the intermediate timescale (on M \M ), and the slowtimescale (on M ). Definition 2 and Proposition 5 are summarised in Figure 5, where we recall that the fast, inter-mediate, and slow dynamics are given by the limiting systems in (8), (19), and (21), respectively.3In Figure 6, we illustrate how the singular dynamics in the three-timescale limit differs from its twotwo-timescale counterparts. (a) ε = 0, δ = O (1). (b) ε = 0 = δ . (c) ε = O (1), δ = 0. Figure 6: Examples of singular cycles for the cases of (a) one fast and two slow variables; (b) onefast, one intermediate, and one slow variable; and (c) two fast and one slow variables. In (a), thereduced flow on M is not affected by a supercritical manifold M : trajectories are not necessarilyattracted to folded singularities, which are now folded nodes; rather, they can jump as soon as theyreach L ∓ . In (c), trajectories are not affected by fold lines L ∓ of M ; rather, they are attractedto Z a ∓ , jumping as soon as they reach the fold points p ∓ . In the three-timescale limit in (b), theflow exhibits characteristics of the two scenarios described in (a) and (c). In this section, we discuss the correspondence between the families of singular cycles constructedin Proposition 5 and the MMO trajectories which perturb from those cycles for ε and δ positive,but sufficiently small, in Equation (1). In particular, we give a qualitative characterisation of theresulting MMO dynamics in dependence of system parameters. For ε and δ positive and sufficiently small in (1), MMO trajectories can be composed from com-ponents that evolve close to fast, intermediate, and slow segments of the corresponding singularcycles, as discussed in Section 2. In a first approximation, where the fast and intermediate segmentsare approximated as straight lines – the latter in the ( x, z )-plane – trajectories are attracted to thevicinities of both folded singularities q ∓ if these are aligned or connected, and only to one of themif they are remote, as can be seen from Figure 5. (In Section 2, we showed that the funnels ofthe folded nodes q ∓ for Equation (1) stretch with decreasing δ ; in the three-timescale limit as δ approaches zero, these funnels can be viewed as having been stretched “infinitely” in one direction.)From the well-established theory of two-timescale singular perturbations, it is known that SAOsarise in the passage past folded singularities [6, 4, 22] under the perturbed flow; the underlying localmechanisms are well understood. We discuss the three-timescale analogues of these mechanismsin Section 3.2 below. In particular, we conclude that SAOs are observed “above” or “below”, inthe language of Figure 1, depending on which folded singularity of Equation (1) trajectories areattracted to; double epochs of SAOs can occur when trajectories are attracted to both folded singu-larities q ∓ . The mixed-mode dynamics of Equation (1) can hence naturally be classified accordingto whether the folded singularities q ∓ are remote, aligned, or connected; cf. Figure 7 and Figure 8.In a first step we note that, under Assumption 1 and Assumption 2, (1) undergoes Hopf bi-furcations for two values of the parameter µ in (1c), which we denote by µ ∓ SH ; these bifurcations4are referred to as “singular Hopf bifurcations” in the literature [13, 6] and separate the regions ofoscillatory dynamics from steady-state behaviour in Equation (1). We note that, for ε = 0 = δ ,true equilibria of (1) cross the folded singularities q ∓ at µ ∓ SH = − φ ( x ∓ q , y ∓ q , z ∓ q ); in particular, As-sumption 2 implies that 0 = µ − SH < µ + SH in the double singular limit. We consider the asymptoticsof µ ∓ SH for ε positive and sufficiently small in the following Section 3.2.In the case where the folded singularities of M in (1) are remote, it follows that the perturbedflow for µ ∈ ( µ − SH , µ + SH ) exhibits either MMOs with single epochs of SAOs, or “two-timescale”relaxation oscillation where the flow alternates between the fast and the intermediate dynamics.In [14], it was shown that to leading order in δ , the µ -values which separate the correspondingparameter regimes can be determined by requiring that the intermediate flow on S a − is “balanced”by that on S a + . Thus, by Lemma 1, the value µ − r that is associated with µ − SH is found by solving G ( x , x max ; 0; µ ) + G ( x ∗ max ,
0; 0; µ ) = 0 (29)for µ ; here, x max = − f f , x ∗ max = f f , and x = − f f . We emphasise that µ − r is independent of ε and δ to the order considered in (29), by construction. Since Equation (16) implies µ − r = − (cid:82) x ∗ max F (cid:48) ( σ ) φ ( σ,F ( σ ) , ασ + βF ( σ ) d σ + (cid:82) x max x F (cid:48) ( σ ) φ ( σ,F ( σ ) , ασ + βF ( σ ) d σ (cid:82) x ∗ max F (cid:48) ( σ ) ασ + βF ( σ ) d σ + (cid:82) x max x F (cid:48) ( σ ) ασ + βF ( σ ) d σ , (30)the position of µ − r relative to µ − SH depends on the properties of φ ( x, y, z ) in (1); in particular,SAOs below may not occur if µ − r falls into the steady-state regime. To distinguish between the twoscenarios, we first show Lemma 2.
If the folded singularities q ∓ of M are remote, i.e., if αβ < f f , then < (cid:90) x max x F (cid:48) ( σ ) ασ + βF ( σ ) d σ < (cid:90) x ∗ max F (cid:48) ( σ ) ασ + βF ( σ ) d σ (31) Proof.
For σ ∈ ( x ∗ max , F (cid:48) ( σ ) < ασ + βF ( σ ) < (cid:82) x ∗ max F (cid:48) ( σ ) ασ + βF ( σ ) d σ >
0. For σ ∈ ( x max , x ), we again note that F (cid:48) ( σ ) <
0, as well as that (cid:82) x max x F (cid:48) ( σ ) ασ + βF ( σ ) d σ = − (cid:82) x x max F (cid:48) ( σ ) ασ + βF ( σ ) d σ ; moreover, we claim that ασ + βF ( σ ) >
0. To see that, we recall that ασ + βF ( σ ) = ασ + βf σ + βf σ = σ ( α + βf σ + βf σ ). The quadratic function p ( σ ) = α + βf σ + βf σ is concave up due to βf >
0, which implies that p ( σ ) is positive for σ > x max ifand only if p ( x max ) ≥
0. That condition is equivalent to αβ ≤ f f , which is satisfied since we assumethat the folded singularities q ∓ are remote; recall Proposition 4. It follows that σp ( σ ) > σ ∈ ( x max , x ) and, hence, that (cid:82) x max x F (cid:48) ( σ ) ασ + βF ( σ ) d σ >
0. The ordering in (31) is straightforward.Lemma 2 immediately implies the following.
Proposition 6.
If the function φ ( x, y, z ) is such that the numerator in Equation (30) is negative,then we have µ − SH < µ − r . If the numerator in Equation (30) is negative, then there exists µ − r > µ − SH that separates MMOswith single epochs of SAOs from relaxation oscillation; see panel (a) of Figure 7. If, on the otherhand, that numerator is positive, then the value µ − r < µ − SH is in the steady-state regime. For µ > µ − SH , the slow drift (in z ) after one return – which is given by the left-hand side in (29) –5is positive, which means that trajectories do not reach Z a − again; therefore, relaxation oscillationis observed, as shown in panel (b) of Figure 7. Finally, for ε = 0, it follows from Lemma 1,Equation (29), and Lemma 2 that the implicit function theorem applies; hence, µ − r will perturbsmoothly in δ such that a singular (in ε ) cycle passes through q − and q + for δ sufficiently small.A µ -value µ + r that is associated with µ + SH can be obtained in a similar fashion; an illustration isgiven in panels (c) and (d) of Figure 7. As our focus is on non-trivial MMO dynamics in Equation (1)here, and motivated by the Koper model which we study in Section 4 below, we will henceforthalso assume the following: Assumption 3.
The function φ ( x, y, z ) in (1c) is such that, for ε = 0 = δ , µ − SH < µ − r < µ + r < µ + SH (32) holds. (a) (b)(c) (d) Figure 7: Dynamics of Equation (1) in the case of remote singularities, with f , f , α , and β fixed and ε = 0 = δ : there exist µ -values µ ∓ SH that distinguish between oscillatory dynamics andsteady-state behaviour; additionally, there exist two values µ ∓ r which potentially separate MMOswith single epochs of SAOs from relaxation oscillation, in dependence of the properties of φ ( x, y, z )in (1c).When the folded singularities of M are aligned or connected, “double epochs” of slow dynamicsare observed for µ ∈ ( µ − SH , µ + SH ) in (1); see Figure 8 and the bottom row in Figure 9. We note thatthe “double epoch” regime can be further divided into cases where SAOs occur “above” and “below”;SAOs are seen “above” with SAO-less slow dynamics below or vice versa; or “three-timescale”relaxation oscillation is found with the flow alternating between fast, intermediate, and slow SAO-less dynamics. A precise characterisation is dependent on the properties of the function φ ( x, y, z )in (1c) and hence requires a case-by-case study; see also Section 3.3 below.6Figure 8: Dynamics of Equation (1) in the case of aligned/connected singularities, with f , f , α ,and β fixed and ε = 0 = δ : there exist µ -values µ ∓ SH that distinguish between oscillatory dynamicsand steady-state behaviour.In summary, the emergence of MMO dynamics in (1) can thus be understood as follows. Bystandard GSPT [10], the normally hyperbolic portions S a ∓ and Z a ∓ of M and M , respectively,perturb to S a ∓ εδ and Z a ∓ εδ , respectively. Given an initial point ( x, y, z ) ∈ S a − εδ , the correspondingtrajectory will follow the intermediate flow on S a − εδ until it is either attracted to Z a − εδ or until itreaches the vicinity of L − . If the trajectory is attracted to Z a − εδ , then it follows the slow flow thereonand can undergo SAOs; if it reaches the vicinity of L − , then there is no slow dynamics, and thetrajectory jumps to the opposite attracting sheet S a + εδ , resulting in a large excursion. The abovesequence then begins anew; see Figure 9 for a schematic illustration: depending on the relativegeometry of the folded singularities q ∓ of M , single or double epochs of slow dynamics will occur,as indicated in Figure 8. In the remainder of this section, we will outline in some detail theunderlying local mechanisms; then, we will comment on the persistence of the global classificationof the resulting MMO trajectories from Figure 7 and Figure 8 for ε and δ positive, but sufficientlysmall. In this subsection, we discuss the emergence of SAOs in a vicinity of L ∓ in (1) when trajectoriesare attracted to Z a ∓ , respectively; we focus on describing the properties of Z a − close to L − here,as the description of Z a + near L + is analogous.We first consider the partially perturbed fast Equation (6) with ε sufficiently small and δ = 0: x (cid:48) = − y + f x + f x , (33a) y (cid:48) = ε ( αx + βy − z ) , (33b) z (cid:48) = 0 . (33c)By standard GSPT [10, 15], we can define slow manifolds S a,rε for (33) as surfaces that are foliatedby orbits within { z = z } , with z constant. Since the steady states of (33) correspond to portionsof the supercritical manifold M , it follows that Z a,rε ≡ Z a,r , i.e., that the geometry of Z a,rε ≡ Z a,r is, in fact, ε -independent. However, since it will become apparent that the stability properties of Z a,rε do depend on ε , we will not suppress the ε -subscript in our notation.For ε sufficiently small, Equation (33) undergoes a Hopf bifurcation at a point p − DH = (cid:0) x − DH , y − DH ,z − DH (cid:1) ; the periodic orbits that arise in that bifurcation cease to exist at z − CN , where a connecting7 (a) Remote singularities. (b) SAOs “below” ( k = − . λ = − . k = − . λ = 0 . Figure 9: Schematic illustration of the emergence of MMO trajectories in (1).Figure 10: Stability of the supercritical manifold M on the various portions of Z a − , for ε suffi-ciently small and δ = 0: at z − DN ∓ , the real eigenvalues of the linearisation of Equation (33) about M in (38) become complex, with a corresponding change from nodal to focal attraction or repul-sion, and vice versa; at z − DH , a Hopf bifurcation occurs which gives rise to small-amplitude periodicorbits. These orbits cease to exist at z − CN , where a connecting trajectory between S a − ε and S rε isfound. We note that the corresponding z -interval is of width O ( ε ), while the spiralling regime is O ( √ ε ) wide overall.8trajectory between S a − ε and S rε is found. In other words, S a − ε and S rε intersect transversely withinthe hyperplane P − CN : { z = z − CN } which lies O ( ε )-close to p − DH in the z -direction. Moreover, twodegenerate nodes p − DN ∓ are located on M around p − DH at an O ( √ ε )-distance; see Figure 10. Theasymptotics (in ε ) of these objects is summarised below. Lemma 3.
A Hopf bifurcation of Equation (33) occurs at p − DH : (cid:0) x − DH , y − DH , z − DH (cid:1) ∈ Z a − ε , where x − DH = − β f ε + O ( ε ) , y − DH = β f ε + O ( ε ) , and z − DH = − αβ f ε + O ( ε ) . (34) Two degenerate nodes p − DN ∓ are located at x − DN ∓ = ∓ (cid:18) √ αf √ ε + β f ε (cid:19) + O ( ε ) , y − DN ∓ = αf ε + O ( ε ) , and z − DN ∓ = ∓ α f √ ε + αβf (cid:18) ∓ (cid:19) ε + O (cid:0) ε (cid:1) , (35) while a canard trajectory is contained in the hyperplane P − CN : { z = z − CN } , with z − CN = z − DH + αβ f − − αf )4 (1 + f ) f ε + O ( ε ) . (36) Proof.
The hyperplane { z = z − CN } , which contains the transverse intersection between S a − ε and S rε , can be obtained by Melnikov-type calculations; see [15, 18]. The remaining estimates followby considering the Jacobian matrix of the linearisation of (33) along M , J = (cid:18) f x + 3 f x − εα εβ (cid:19) , (37)the eigenvalues of which are ν , = 12 (cid:20) βε + 2 f x + 3 f x ± (cid:113)(cid:0) βε + 2 f x + 3 f x (cid:1) − (cid:0) αε + 2 βεf x + 3 βεf x (cid:1)(cid:21) . (38) Remark 3.
The Hopf bifurcation at p − DH is “inherited” from the fact that M and L − intersect inthe folded singularity q − : for ε = 0 = δ , the trace of the Jacobian J vanishes at that point. It is therefore now apparent how the stability properties of Z a,rε depend on ε in spite of thegeometry being identical to that of Z a,r , respectively. In addition, we recall that Equation (1) un-dergoes Hopf bifurcations for µ ∓ SH = − φ ( x ∓ DH , y ∓ DH , z ∓ DH ), where x ∓ DH , y ∓ DH , and z ∓ DH are estimatedin Lemma 3.Motivated by Lemma 3, we introduce the following notation: for δ = 0, we define the intervals I innod = (cid:16) −∞ , z − DN − (cid:17) , I inspir = (cid:16) z − DN − , z − DH (cid:17) , and I can = (cid:0) min (cid:8) z − DH , z − CN (cid:9) , max (cid:8) z − DH , z − CN (cid:9)(cid:1) . (39)Then, it follows that91. the manifold S a − ε connects to Z a − for z < min (cid:8) z − DH , z − CN (cid:9) , while S rε connects to Z a − for z > max (cid:8) z − DH , z − CN (cid:9) ;2. for f < (1 − αf ) ( f > (1 − αf )), i.e., for z − CN > z − DH ( z − CN < z − DH ), the Hopf bifur-cation at p − DH is supercritical (subcritical), with the resulting periodic orbits the ω -limit sets( α -limit sets) of trajectories on S a − ε ( S rε ).The resulting geometry is illustrated in Figure 10; we emphasise that analogous objects p + DH , P + CN , and p + DN ± , which are located symmetrically to the above, exist on Z a + .As has already been pointed out in [18], the Hopf point p − DH and the canard point p − CN on Z a − collapse to the origin in the limit of ε = 0; correspondingly, the origin is referred to as the “canarddelayed Hopf singularity” in the singular limit of ε = 0 = δ . As a result, the folded singularity at q − displays characteristics of both a Hopf point – in that the trace of the Jacobian in (37) vanishes– and a canard point – in that S a − and S r meet along a fold – in the double singular limit.We briefly describe the associated two mechanisms – bifurcation delay and sector-type dynam-ics – in the following; we remark that the former is common in two-timescale systems with twofast variables, while the latter typically occurs in two-timescale systems with two slow variables.Therefore, the coexistence of these mechanisms in three-timescale systems is due to the fact thatsuch systems can simultaneously be viewed as having two fast and one slow variables, as well asas one fast and two slow variables. (For four-dimensional two-timescale systems with two fast andtwo slow variables, that interplay has been documented in [2].) Bifurcation delay is typically encountered in two-timescale systems with two fast variables andone slow variable. In the context of Equation (1), it is realised when trajectories are attractedto Z εδ (cid:12)(cid:12) z ∈I innod + O ( δ ) or Z εδ (cid:12)(cid:12) z ∈I inspir + O ( δ ) , recall (39) and Figure 10. Following the slow flow on Z a − εδ ,trajectories experience a delay in being repelled away from Z a − εδ when crossing the Hopf bifurcationpoint p − DH , as the accumulated contraction to Z a − εδ needs to be balanced by the total expansion from Z a − εδ [16]. Specifically, given some point p in = ( x in , y in , z in ) ∈ Z a − εδ , one obtains the x -coordinate ofthe corresponding point p out from (cid:90) x out x in (cid:60) { ν , ( x ) } µ + φ ( x, F ( x ) , G ( x )) d x = 0; (40)here, ν , are the eigenvalues of the linearisation of Equation (33) about Z a − εδ , as defined in (38).Trajectories that are attracted to Z εδ (cid:12)(cid:12) I inspir typically exhibit “dense” SAOs with initially decreasingand then increasing amplitude; see panel (e) of Figure 12 below for an illustration in the contextof the Koper model, Equation (2). By contrast, trajectories that are attracted to Z εδ (cid:12)(cid:12) I innod arecharacterised by very few SAOs that are followed by a large excursion; cf. Figure 12(c).The case where trajectories enter the spirally attracting regime I inspir is naturally studied inthe “rescaling chart” κ which is introduced as part of a blow-up analysis in [15, 18], since thatregime is bounded by the degenerate nodes p − DN ∓ and, thus, of width O ( √ ε ). In that case, theeigenvalues ν , in (38) are complex conjugates, which implies that the corresponding trajectory of(1) undergoes damped oscillation towards Z a − εδ .On the other hand, when trajectories enter the nodally attracting regime I innod , the correspond-ing entry point is typically O ( ε c ) away from the folded singularity q − , with c < /
2. One may0therefore refer to the unscaled system, Equation (1), for the study of that case. The eigenval-ues ν , in (38) correspond to strong and weak eigendirections: specifically, for z < z − DN − , theeigenvalue ν represents the weak eigendirection, while the eigenvalue ν corresponds to the strongeigendirection; that correspondence is reversed for z > z + DN − . Due to the hierarchy of timescalesin (1), trajectories are first attracted to S a − εδ and then to Z a − εδ . Therefore, for initial conditions( x, y, z ) ∈ S a − εδ , trajectories approach Z a − εδ along the weak eigendirection, while for ( x, y, z ) ∈ S rεδ ,trajectories are repelled from Z a − εδ along the strong eigendirection. It hence seems reasonable tobalance the accumulated contraction and expansion using solely ν in (40). Since the accumulatedcontraction on the intermediate timescale has to be balanced by expansion on the fast timescale,we make the following Claim 1.
Assume that Assumption 1 and Assumption 2 hold, and consider ( x in , y in , z in ) ∈ Z a − εδ | I inspir ∪I innod .Then, the exit point ( x out , y out , z out ) that is defined by (40) satisfies x out < x − DN + + o (1) , y out < y − DN + + o (1) , and z out < z − DN + + o (1) . Remark 4.
In [18], the weak contraction towards Z a − εδ is balanced by the weak expansion therefromvia (cid:90) x DN x in (cid:60) { ν } d x + (cid:90) x out x DN (cid:60) { ν } d x = 0 . In that context, the fold point p − was in fact identified as the “buffer point” at which trajectorieshave to leave Z a − εδ , which is not in agreement with Claim 1. One possible reason for the discrepancyis that the underlying blow-up analysis in [18] is performed entirely in the rescaling chart; such ananalysis can, however, only be valid o (1) -close to p − DH , which is not the case for nodal entry. Remark 5.
The estimates on the entry point p out in Claim 1 can be refined under the additionalassumption that the slow flow of Equation (1) is constant, i.e., that φ ( x, y, z ) = 0 : as in [4], it thenfollows from (40) that x out = x DH − x in . Sector-type dynamics is typically encountered in two-timescale systems with one fast variable andtwo slow variables; it can be described by exploiting the near-integrable structure of Equation (1)in a vicinity of the canard point p − CN [14, 4]. Sector-type dynamics is realised when trajectoriesare attracted to Z (cid:12)(cid:12) z ∈I can + O ( δ ) , where I can is given by (39). (We emphasise that, for δ sufficientlysmall, S a − εδ and S rεδ intersect in a canard trajectory that provides a connection between the twomanifolds; recall Section 3.2.) For ε and δ sufficiently small and z in ∈ I can + O ( δ ), trajectoriesremain “trapped” and undergo SAOs (“loops”), taking O ( µδ √− ε ln ε ) steps in the z -direction untilthey reach a point p out at which they can escape following the fast flow of Equation (1). The z -coordinate of that point can hence be approximated by z out = z − CN + o (1) . (41)The number of SAOs that is observed in the corresponding trajectory is determined by the passagethereof through sectors of rotation [14], the boundaries of which are so-called “secondary” canards.Trajectories that are attracted to this regime typically exhibit few SAOs of near-constant amplitude;1see panel (b) of Figure 13 below, where sector-type SAOs are seen in between delay-type segments.A detailed study of sector-type dynamics in Equation (1) is part of work in progress; see again [14]for an in-depth discussion in the context of their prototypical model, Equation (3). Remark 6.
Let z q + denote the z -coordinate of the folded singularity q + , as defined in Equation (13).If z q + ≥ z − CN + O ( δ ) , as in Figure 5(a), a “complete” canard explosion occurs at z − CN . On the otherhand, for z − q ≤ z − CN + O ( δ ) , cf. Figure 5(b), the canard explosion is “incomplete” [5]: the layerproblem, Equation (8), has two equilibria, with the equilibrium corresponding to Z r being a saddle.Hence, a homoclinic connection is formed from that saddle to itself in the intersection between S a − ε and S rε . The implications of these two different scenarios, both for the local near-integrabledynamics of Equation (1) and the global oscillatory dynamics, are currently under investigation.We emphasise that, in Equation (4) [18], a homoclinic connection is always present between thecorresponding slow manifolds; our analysis differs from that in [18, Section 4.3] due to the presenceof an O ( x ) -term in (1a), as is also evident from (36). Finally, we combine the non-local and local analyses presented in this section thus far in orderto give a holistic description of the MMO dynamics in Equation (1); in particular, we conjecturethe persistence of the families of singular cycles constructed in Section 2 for ε and δ positive andsufficiently small. (We recall that Figure 7 is merely indicative for the perturbed dynamics of (1),as µ ∓ r are approximated in the double singular limit of ε = 0 = δ .)We first focus on the case where the folded singularities of M are remote: Conjecture 1.
Assume that Assumption 1, Assumption 2 and Assumption 3 hold, and that thesingularities of (1) are remote. Then, there exist ε and δ positive and small such that, for ( ε, δ ) ∈ (0 , ε ) × (0 , δ ) , Equation (1) exhibits either1. MMOs with single epochs of SAOs and associated Farey segments of the type L s or L s , re-spectively, with s > , for µ in appropriate subintervals of ( µ − SH , µ − r + O ( ε + δ )) ∪ ( µ + r + O ( ε + δ ) , µ + SH ) ; or2. two-timescale relaxation oscillation for µ in an appropriate subinterval of ( µ − r , µ + r ) + O ( ε + δ ) . In other words, we conjecture that the corresponding families of singular cycles will persist ina full neighborhood not only of δ = 0, but also of ε = 0, in (1), uniformly in ε and δ . Recall thatthe asymptotics of µ ∓ SH for ε and δ sufficiently small is given in Section 3.2 above. In that case,we are also able to derive quantitative results on the structure of the resulting MMO trajectories;thus, for instance, we have the following result on the number of LAOs that will occur after anSAO segment: Proposition 7.
Let ε and δ be sufficiently small, and assume that the folded singularities of M are remote and that Assumption 1 and Assumption 2 hold. Given a trajectory that is repelled awayfrom Z a − εδ at some point p out with z out > , let L denote the number of large excursions that followbefore the trajectory is again attracted to Z a − εδ .1. If z in < , then L = 1 ;2. if < z in < z out , then L = 1 + (cid:22) z out δ ( G ( x , x max , µ ) + G ( x ∗ max , , µ )) (cid:23) , (42) where (cid:98) (cid:99) denotes the floor function; and
3. if z in ≥ z out , then the trajectory undergoes relaxation oscillation.Proof. For z in <
0, the given trajectory is attracted to Z a − and then undergoes SAOs before beingrepelled away again, which implies L = 1. On the other hand, if 0 < z in < z out , then the trajectoryreaches L − and undergoes large excursions, or “jumps”; the estimate for the total number of thoseis given by the ratio of the total drift in the z -direction, i.e., of z out , over the drift z in − z out aftereach step. Finally, if z in ≥ z out , then the drift in z is positive, which implies that trajectories reach L − after every return. The situation is summarised in Figure 9.Further such results can be derived in a similar fashion.Similarly, we make the following conjecture for the case where the folded singularities in (1) areeither aligned or connected: Conjecture 2.
Assume that Assumption 1 and Assumption 2 hold, and that the singularitiesof (1) are aligned or connected. Then, there exist ε and δ positive and small such that, for ( ε, δ ) ∈ (0 , ε ) × (0 , δ ) and µ in an appropriate subinterval of ( µ − SH , µ + SH ) , Equation (1) exhibitseither1. MMOs with double epochs of SAOs and associated Farey segment of the type s k , with s, k > ;2. MMOs with one epoch of SAOs and one non-SAO epoch of slow dynamics and associatedFarey segment of the type s or s for s > ; or3. three-timescale relaxation oscillation and associated Farey segment of the type . Conjecture 1 and Conjecture 2 are associated with panels (a) and (c) of Figure 7 and Figure 8,respectively, which correspond to the singular limit of ε = 0 = δ and which hence classify toleading order the MMO dynamics of Equation (1). While these conjectures are substantiated bothconceptually and numerically here, a rigorous proof, which is left for future work, would necessitatea detailed description of the corresponding transition maps and their asymptotics on the fast,intermediate, and slow timescales, as was also done in [14].It might not be immediately apparent why aligned singularities yield MMO dynamics withdouble epochs of SAOs; however, we recall that, for ε and δ sufficiently small, trajectories jumpunder the fast flow of (1) after crossing a bifurcation point p ∓ DH , recall (40), which implies thatthey will likewise be attracted to Z a ∓ εδ on the opposite sheet S a ∓ εδ .Moreover, we emphasise that the “double epoch” regime in panel (b) of Figure 8 does notnecessarily imply MMO dynamics with two epochs of SAOs but, rather, with double epochs of perturbed slow dynamics of the corresponding singular cycles. That is, MMO trajectories areattracted to the vicinity of both branches Z a ∓ and hence exhibit slow dynamics; however, whetherSAOs will occur depends on which regime of Z trajectories enter, by (39). In particular, if atrajectory is attracted to the spiralling region on both Z a − and Z a + , then two epochs of SAOsare observed; cf. Figure 12(c). On the other hand, trajectories that are first attracted to thespiralling region on, say, Z a − , and are then attracted to and repelled from the nodal region on Z a + , feature SAOs below and mere slow dynamics above; see Figure 16(a). (The correspondingsegment of the associated Farey sequence would be 1 s , with s > Z a − and Z a + features no SAOs at alland is hence a relaxation oscillation with fast, intermediate, and slow components; the associatedFarey sequence would be 1 . In the transition between remote and connected singularities, exotic3MMO trajectories may occur which exhibit segments of two-timescale relaxation oscillation, SAOsabove, and SAOs below; cf. Figure 16(c). (The associated Farey sequence would be 1 s L k , with L, s, k > φ in (1c); itis hence not feasible to further subdivide that region in Figure 8 exclusively on the basis of systemparameters in Equation (1). Rather, a case-by-case study is required.It is essential to note that the above conjectures refer to a fixed geometry of Equation (1)and, in particular, to the function φ satisfying Assumption 2; furthermore, we emphasise that theresulting MMO trajectrories cannot, strictly speaking, be viewed as perturbations of particularsingular cycles but, rather, that the latter can yield qualitative information on the former. Finally,we remark on the role of the ratio between the scale separation parameters ε and δ for the dynamicsof Equation (1). Locally, in order for the system to exhibit three timescales and for the iterativereduction from the fast via the intermediate to the slow dynamics to be accurate, ε and δ need tobe sufficiently small, which is akin to asking “When is ε small enough?” in a two-timescale system.We recall that, by Lemma 3, the width of the various regimes on Z is either O ( ε ) or O ( √ ε ). By[14] and Lemma 1, the “step” in the z -direction taken by trajectories after a large excursion andre-injection is O ( δ ); it therefore follows that if δ = O ( ε c ) (0 < c < Z a − | I can , since the width of the latter is O ( ε ). Hence, delay-type SAOs areexpected to dominate in that case; see Figure 13. As the first realisation of our extended prototypical model, Equation (1), we consider the Kopermodel from chemical kinetics, Equation (2): after an affine transformation, the latter can be writtenin the form of (1), with ε = (cid:15) | k | , f = 3 | k | , f = − | k | , (43a) α = 1 , β = − , (43b) µ = k + λ + 2 k , and φ ( x, y, z ) = − y − z ; (43c)henceforth, we will refer to (1) with the above choice of parameters as the Koper model; here, k < λ ∈ R will be our bifurcation parameters.From Section 2, it is apparent that the effect of the parameter k on the dynamics is moresubstantial than that of λ , since variation in k simultaneously affects the timescale separation(through ε ) and the singular geometry (through f and f ), as well as the slow flow and the globalreturn (through µ ). Given k < λ only affects the slow flowand the global return (through µ ). It is therefore the parameter k that determines whether thefolded singularities in the Koper model are remote or aligned/connected, and whether the modelcan exhibit single or double epochs of SAOs. For given k <
0, the parameter λ can differentiatebetween steady-state and oscillatory dynamics, as well as between MMOs and relaxation in thecase of remote singularities. Remark 7.
Alternatively, the Koper model can be written in the symmetric form (cid:15) ˙ x = y − x + 3 x, ˙ y = kx − y + λ ) + z, ˙ z = δ ( λ + y − z ) , which is invariant under the transformation ( x, y, z, λ, k, t ) → ( − x, − y, − z, − λ, k, t ) [6]. In the following, we will restrict to the case where λ > L − only: by Remark 7, the behaviour near L + for λ < t -axis for k fixed and λ → − λ . The critical and supercritical manifolds M and M , respectively, for the Koper model are givenby M = (cid:26) ( x, y, z ) ∈ R (cid:12)(cid:12)(cid:12) y = x − x | k | (cid:27) and M = (cid:26) ( x, y, z ) ∈ M (cid:12)(cid:12)(cid:12) z = x − x − x | k | (cid:27) ;see Section 2. The critical manifold M is normally hyperbolic at S = S a ∓ ∪ S r , (44)where S a − = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) x < (cid:9) , S a + = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) x > (cid:9) , and S r = (cid:8) ( x, y, z ) ∈ M (cid:12)(cid:12) < x < (cid:9) . The fold lines of M are located at L − = (cid:8) ( x, y, z ) ∈ R (cid:12)(cid:12) x = 0 , y = 0 (cid:9) and L + = (cid:26) ( x, y, z ) ∈ R (cid:12)(cid:12)(cid:12) x = 2 , y = 4 | k | (cid:27) ; (45)the corresponding folded singularities q ∓ are found at q − = (0 , ,
0) and q + = (cid:18) , | k | , − | k | (cid:19) , (46)respectively. For the relative position of the folded singularities q ∓ , we have the following Proposition 8.
Let ε = 0 = δ . Then, the folded singularities of the Koper model are connected for − < k < , aligned for k = − , and remote when k < − .Proof. The statement follows from Proposition 4 and (43), or by comparison of the z -coordinatesof q − and q + .The supercritical manifold M is normally hyperbolic everywhere except at the fold points p ∓ ,where x ∓ p = 1 ± (cid:113) − | k | , y ∓ p = (cid:16) ± (cid:113) − | k | (cid:17)(cid:16) ∓ (cid:113) − | k | (cid:17) | k | , and z ∓ p = 1 ∓ (cid:113) − | k | − (cid:16) ± (cid:113) − | k | (cid:17)(cid:16) ∓ (cid:113) − | k | (cid:17) | k | . (47)Based on the above, we have the following5 Proposition 9. If − < k < , then M admits two fold points which are located between thepoints of intersection of M with L ∓ , i.e., on the repelling sheet of M . If k < − , then M admits no fold points. We reiterate that, due to k <
0, the fold points p ∓ in the Koper model cannot cross L ∓ , andthat the corresponding singular geometry is therefore as depicted in Figure 4(c). Remark 8.
In [1, Example 4.3], the manifold M is characterised as normally hyperbolic every-where, in spite of its graph being S -shaped. Proposition 9 above shows that M can, in fact, admittwo fold points at which normal hyperbolicity is lost. Here, we classify the dynamics of the Koper model in the three-timescale context for various choicesof the parameters k and λ . In particular, we hence construct the two-parameter bifurcation diagramshown in Figure 11. (A two-timescale analogue of Figure 11, for one fast and two slow variablesin Equation (2), is presented in [6].) Given the definition of µ in (43), we consider λ as a functionof k here when retracing the analysis from Section 3, in particular in relation to the classificationin Figure 7 and Figure 8; the requisite calculations are simplified due to the symmetry of (2), byRemark 7.Figure 11: Two-parameter bifurcation diagram of the Koper model ε and δ sufficiently small:oscillatory dynamics is restricted to the triangular region of the ( k, λ )-plane that is bounded by λ ∓ SH ( k ); MMO dynamics is separated from relaxation oscillation by the curves λ ∓ r = λ ∓ r ( k ); toleading order, the mixed-mode regime is subdivided into regions of either single or double epochsof SAOs at k = − p ∓ DH coincide with true, global equilibria in the governing equations: Proposition 10.
Let ε and δ be sufficiently small, and fix k < . If k < − , then the Koper model undergoes (singular) Hopf bifurcations for λ = λ ∓ SH ( k ) , where λ − SH ( k ) = − (2 + k ) + k ε + O ( δ + ε ) = − (2 + k ) + | k | (cid:15) + O ( δ + (cid:15) ) and (48a) λ + SH ( k ) = − λ − SH ( k ) . (48b) Proof.
The equilibria of Equation (1) are obtained by solving0 = − y + f x + f x = αx + βy − z = µ + φ ( x, y, z ) . (49)Substituting in for x − DH , y − DH , and z − DH from (34), using (43), and solving for λ , one finds (48a).It hence follows that oscillatory dynamics is restricted to the triangular area illustrated in Fig-ure 11. A further subdivision of that area is obtained by noting that MMO dynamics is separatedfrom relaxation oscillation by two curves λ − r ( k ) and λ + r ( k ) = − λ − r ( k ); these are found by substi-tuting (43) into (29) and solving for λ . While analytical expressions for λ ∓ r ( k ) can be obtained bydirect integration, they are quite involved algebraically, and are hence not included here. Theseexpressions imply that, for ε = 0 = δ and k < − λ + SH ( k ) < λ + r ( k ) < λ − r ( k ) < λ − SH ( k ), as wellas that λ ∓ r are asymptotically parallel to λ ∓ SH , respectively, for | k | sufficiently large. Numericalexperiments show that, for ε = O (10 − ) and δ = O (10 − ), the transition between MMO dynamicsand relaxation occurs at λ ∓ r ( k ) + O ( δ ), as is to be expected from (29).Finally, the resulting, chevron-shaped region in which MMO dynamics is observed is furtherdivided into subregions in which either single or double epochs of SAOs are found; to leading orderin ε and δ , that division occurs at k = −
4. Geometrically, the division is due to the fact that thefolded singularities q ∓ in the Koper model are remote for k < −
4, while they are connected when − < k <
0; recall Figure 9. We emphasise that, in the two-timescale context of ε sufficiently smalland δ = O (1) in the Koper model, MMOs with double epochs of SAOs occur in a very narrowregion of the ( k, λ )-plane; that region corresponds to the regime where both folded singularities areof folded-node type and trajectories are attracted to both of them through the associated funnels, by[6]. Here, we have shown that, in the three-timescale context, that parameter regime is “stretched”;hence, trajectories can reach both folded singularities as long as they are attracted to M on both S a ∓ , i.e., as long as the folded singularities are aligned or connected. Remark 9.
Comparing panels (a) and (c) of Figure 7 and Figure 8 with Figure 11, we note thatthe former two figures are combined in the latter, as one-parameter diagrams (in µ ) are merged intoone two-parameter diagram in ( k, λ ) ; correspondingly, parallel lines with µ constant in Figure 7 andFigure 8 are “bent”, and hence intersect, in Figure 11. (Here, we reiterate that k determines thesingular geometry of the Koper model, while λ impacts the resulting flow.) In this subsection, we verify our classification of the three-timescale dynamics of the Koper modelfor various representative choices of the parameters k and λ , as indicated in Figure 11. We initiallyfix ε = 0 .
01 = δ and λ = 0 .
5, and we vary k . We recall that the Koper model is symmetric in λ ,and that it hence suffices to consider positive λ -values; see Remark 7 and Figure 1.For k = − . k = − . q ∓ are connected in that regime. We note that the dynamicson Z a − differs from that on Z a + due to the properties of φ ( x, y, z ) given in (43) in spite of thesingular geometry being symmetric; see Figure 12(c). In addition, trajectories “jump” close to the7degenerate nodes z ∓ DN + , which is in agreement with Claim 1. For k = − . k = − . δ = O ( ε ), their prototypical model, Equation (3), can exhibit MMOswhich contain SAO segments that are the product of bifurcation delay alternating with sector-typedynamics. In Figure 13, we present an example that indicates sector-delayed-Hopf-type dynamicsin the Koper model. We remark that the canard point p − CN coalesces with p − DH when k = − I can hence vanishes, which follows by substitution of (43) into (36). Since,in addition, z CD = O ( √ ε ), a crude requirement for the existence of such mixed dynamics is that δ = O ( | z CN | c ) for c ≥
1. Finally, we note that we typically observe more LAO segments betweenepochs of SAOs for smaller values of δ than for larger ones, by Equation (42) of Proposition 7; seeagain Figure 13.We emphasise that the MMO trajectories described here cannot be viewed, strictly speaking,as perturbations of individual singular cycles, as described in Section 2. Rather, we have shownthat if the folded singularities of (1) are remote, then there exist ε and δ positive and sufficientlysmall such that the Koper model exhibits MMOs with single epochs of SAOs; correspondingly,we observe double epochs of SAOs if those singularities are aligned or connected, as postulated inConjecture 1 and Conjecture 2. The above statement is corroborated by numerical continuation,as illustrated in Figure 14, where multiple periodic orbits seem to coexist for k , λ , ε , and δ fixed.(A similar observation was made in the context of the two-timescale Koper model, i.e., for δ = 1in Equation (2c) [6, Figure 19].) A more detailed study of the properties of these periodic orbits inrelation to the MMO dynamics of Equation (1) is part of work in progress. In our second example, we outline how the results obtained for Equation (1) can be applied toa three-dimensional reduction of the Hodgkin-Huxley equations that was derived by Rubin andWechselberger in [20]: ε ˙ v = ¯ I − (cid:0) v − ¯ E Na (cid:1) m ∞ ( v ) h − ¯ g k (cid:0) v − ¯ E k (cid:1) n − ¯ g l (cid:0) v − ¯ E L (cid:1) , (50a)˙ n = 1 τ n t n ( v ) ( n ∞ ( v ) − n ) , (50b)˙ h = 1 τ h t h ( v ) ( h ∞ ( v ) − h ) , (50c)where t x ( v ) = 1 α x ( v ) + β x ( v ) and x ∞ ( v ) = α x ( v ) α x ( v ) + β x ( v ) for x = m, h, n. Moreover, the functions α x ( v ) and β x ( v ) ( x = m, h, n ) are defined as α m ( v ) = ( v + 40) / − e − ( v +40) / , α h ( v ) = 7100 e − ( v +65) / , α n ( v ) = ( v + 55) / − e − ( v +55) / ,β m ( v ) = 4 e − ( v +65) / , β h ( v ) = 11 + e − ( v +35) / , and β n ( v ) = 14 e − ( v +65) / ;finally, the corresponding parameters in (50) are defined as¯ I = Ik , ¯ g k = 0 . , ¯ g l = 0 . , ¯ E Na = 0 . , ¯ E k = − . , ¯ E L = − . , and ε = 0 . , I is the applied current in µ A / cm in the original Hodgkin-Huxley equations [9], while k = (cid:0)
120 mS / cm (cid:1) k v , with k v = 100 mV . Following [9], we set τ n = 1 in (50) and assume that τ h (cid:29) δ = τ − h ,we obtain the three-timescale system ε ˙ v = ¯ I − (cid:0) v − ¯ E Na (cid:1) m ∞ ( v ) h − ¯ g k (cid:0) v − ¯ E k (cid:1) n − ¯ g l (cid:0) v − ¯ E L (cid:1) , (51a)˙ n = 1 t n ( v ) ( n ∞ ( v ) − n ) (51b)˙ h = δ t h ( v ) ( h ∞ ( v ) − h ) , (51c)where v , n , and h are the fast, intermediate, and slow variables, respectively. (A similar geometryis obtained for τ h = 1 and τ n (cid:29) h is the intermediate variable and n is theslow variable, while v is still fast; see again [9] for details.)In the notation of Section 2, the critical and supercritical manifolds M and M of (51) aregiven by M = (cid:110) ( v, n, h ) ∈ ( ¯ E k , ¯ E Na ) × (0 , (cid:12)(cid:12) V ( v, n, h ) = 0 (cid:111) and M = (cid:8) ( v, n, h ) ∈ M (cid:12)(cid:12) V ( v, n ∞ ( v ) , h ) = 0 (cid:9) , respectively, where V ( v, n, h ) = ¯ I − (cid:0) v − ¯ E Na (cid:1) m ∞ ( v ) h − ¯ g k (cid:0) v − ¯ E k (cid:1) n − ¯ g l (cid:0) v − ¯ E L (cid:1) . The critical manifold M is three-dimensional, S -shaped, and has a fold set F M which is the unionof two curves L − and L + : F M = L − ∪ L + = (cid:26) ( v, n, h ) ∈ M (cid:12)(cid:12)(cid:12) ∂V∂v ( v, n, h ) = 0 (cid:27) ;see Figure 15. (In particular, it is shown in [20] that L − and L + are disjoint as long as h isbounded away from zero.) The manifold M hence consists of a repelling middle sheet S r andtwo attracting outer sheets S a ∓ . The supercritical manifold M is two-dimensional, S -shaped,and admits a pair of fold points. Figure 16 illustrates the projection of the fold lines L − and L + and of the supercritical manifold M onto the ( v, h )-plane for three different values of the appliedcurrent I , which is identified as the natural bifurcation parameter in [20]. (We note that they statetheir results in terms of the original current I , rather than of its rescaled counterpart ¯ I , for ease ofcomparison with [9].)In Figure 16, we choose ε = 0 . τ h = 45, which implies δ h = 0 . I is varied in the reduced three-timescaleHodgkin-Huxley model, Equation (51). We remark that the time series illustrated in panels (a),(c), (e), and (g) of Figure 16 for various values of I have been documented in [9]; however, theunderlying geometry was not emphasised there. We also remark that “exotic”, and potentially evenchaotic, mixed-mode dynamics is possible in the transition between the different geometries, seeFigure 16(c); the corresponding trajectories consist of an alternation of SAOs “above” and “below”with relaxation oscillation and bursting-type segments, and are observed at the transition fromMMOs with double epochs of SAOs to those with single epochs. Finally, one observes that M L − than it does near L + , which is due to the difference in slope of M there. An in-depth geometric analysis of themultiple-timescale Hodgkin-Huxley equations, in their formulation due to [9], is part of work inprogress; cf. the upcoming article [12]. In the present work, we have introduced an extended prototypical example of a three-dimensional,three-timescale system, Equation (1). We have focused on the various types of oscillatory dynamicsthat can arise in dependence of the geometry of our system, and we have shown how transitionsbetween those occur. In particular, in Section 3, we identified the geometric mechanism that isresponsible for the transition from MMOs with single epochs of SAOs to those with double epochs,and we argued that those are robust in the three-timescale context. Specifically, we claim that, ifthe folded singularities of (1) are remote, then there exist ε and δ sufficiently small such that thesystem exhibits MMOs with single epochs of SAOs, whereas double epochs of SAOs can be observedif the singularities are aligned or connected; cf. Proposition 4. In Section 4, we demonstrated ourresults for the Koper model from chemical kinetics [13], which represents one particular realisationof our prototypical system; in particular, we constructed the two-parameter bifurcation diagram inFigure 11 on the basis of results obtained in Section 3, thus classifying in detail the mixed-modedynamics of the three-timescale Koper model. Then, in Section 5, we indicated how some of ourfindings extend to a three-dimensional reduction of the Hodgkin-Huxley equations that were derivedby Rubin and Wechselberger [20].A posteriori, it is evident that the local dynamics of our extended prototypical example, Equa-tion (1), is similar to that of the canonical system, Equation (4), proposed in [18]; however, dueto the presence of cubic x -terms in (1a), MMO trajectories are found in (1), but not in (4). Theprototypical system in Equation (3), on the other hand, can only exhibit MMOs with single epochsof SAOs, as opposed to our extended Equation (1), which is due to the absence of y -terms in (3b);recall Proposition 4. We additionally remark that the singular geometry considered here is rela-tively specific due to its symmetry properties. In particular, we have not considered explicitly thescenario where the fold points p ∓ of M cross the fold lines L ∓ of M ; recall, in particular, panels(b) and (e) of Figure 4. While that scenario is not realised in the Koper model, Equation (2), ithas been shown to give rise to interesting local dynamics through the interaction of p ∓ with L ∓ ; arecent, relevant example can be found in [7].Our analysis in Section 3 shows that, in parameter regimes where both M and M are normallyhyperbolic, standard GSPT [10] implies that an iterative reduction of timescales can be applied. Inthe fully perturbed Equation (1) with ε and δ sufficiently small, it follows that Z a,rεδ lie O ( δ )-close totheir unperturbed counterparts Z a,r , since Z a,rε are ε -independent. Since, moreover, the manifolds S a,rε lie O ( ε )-close to S a,r [10], any fibers of Z a,rε,δ that lie on S a,rεδ are O ( ε + δ )-close to S a,r . (Thatestimate is in disagreement with [1]; however, we note that, away from Z a,rεδ , S a,rεδ are O ( ε )-closeto S a,r .) Under Assumption 2, trajectories that are attracted to Z a,rεδ follow the slow flow of (15)and potentially undergo SAOs. In the context of (1), the mechanisms that generate these SAOsare “bifurcation delay” [16, 4, 18] and “sector-type” dynamics [14, 4].With regard to regions where normal hyperbolicity of M is lost, we reiterate that the dynamicsof Equation (1) combines features of two-timescale slow-fast systems with either two slow variablesand a fast one, or one fast variable and two slow ones. As shown in Section 3, the correspondingmechanisms hence coexist and interact, giving rise to complex local dynamics in the vicinity ofthe fold lines L ∓ in (1). We briefly sketched the implications of that interaction; in particular, we0related the emergence of canard-type SAOs to the perturbation of an integrable system [14]. Amore rigorous description of the resulting near-integrable system in the context of the Koper model,Equation (2), is part of work in progress. Of particular interest here is the investigation of Shilnikov-type homoclinic phenomena, as well as the further classification of MMOs with single epochs ofSAOs; specifically, we conjecture that the bifurcation diagram in Figure 11 may be refined, in thatone can identify regions of chaotic mixed-mode dynamics in dependence of the various parametersin the model, as well as on the ratio of ε and δ .We emphasise that, strictly speaking, the MMO trajectories described in Section 3 do notcorrespond to perturbations, for ε and δ positive, of the individual singular cycles constructedin Section 2. Rather, the latter determine the qualitative properties of the former, for ε and δ sufficiently small, as is evident from Figure 14 in the context of the three-scale Koper model, whereseveral periodic orbits seem to coexist for a given choice k , λ , ε , and δ . Similarly, in the reducedHodgkin-Huxley model, Equation (51), exotic Farey sequences are observed, in addition to the onesillustrated in Section 4, which again attests to the fact that one cannot consider MMO trajectoriesas perturbations of individual cycles. An example is given in Figure 16(c); see again the upcomingarticle [12] for an in-depth discussion. Acknowledgements
The authors thank Martin Krupa for his critical reading of previous versions of the manuscript andfor constructive feedback, as well as for insightful discussions and relevant references. PK wouldalso like to thank Hinke Osinga for insightful discussions.PK is supported by the Principal’s Career Development Scholarship for PhD Studies of theUniversity of Edinburgh.
References [1]
P. T. Cardin and M. A. Teixeira , Fenichel theory for multiple time scale singular pertur-bation problems , SIAM Journal on Applied Dynamical Systems, 16 (2017), pp. 1425–1452.[2]
R. Curtu and J. Rubin , Interaction of canard and singular Hopf mechanisms in a neuralmodel , SIAM Journal on Applied Dynamical Systems, 10 (2011), pp. 1443–1479.[3]
P. De Maesschalck, E. Kutafina, and N. Popovi´c , Three time-scales in an extendedBonhoeffer–van der Pol oscillator , Journal of Dynamics and Differential Equations, 26 (2014),pp. 955–987.[4]
P. De Maesschalck, E. Kutafina, and N. Popovi´c , Sector-delayed-Hopf-type mixed-mode oscillations in a prototypical three-time-scale model , Applied Mathematics and Compu-tation, 273 (2016), pp. 337–352.[5]
P. De Maesschalck and M. Wechselberger , Neural excitability and singular bifurca-tions , The Journal of Mathematical Neuroscience (JMN), 5 (2015), p. 16.[6]
M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H. M. Osinga, andM. Wechselberger , Mixed-mode oscillations with multiple time scales , SIAM Review, 54(2012), pp. 211–288.1[7]
M. Desroches and V. Kirk , Spike-adding in a canonical three-time-scale model: superslowexplosion and folded-saddle canards , SIAM Journal on Applied Dynamical Systems, 17 (2018),pp. 1989–2017.[8]
E. J. Doedel, A. R. Champneys, F. Dercole, T. F. Fairgrieve, Y. A.Kuznetsov, B. Oldeman, R. Paffenroth, B. Sandstede, X. Wang, and C. Zhang , AUTO-07P: Continuation and bifurcation software for ordinary differential equations . , 2007. Accessed on 09/09/2020.[9] S. Doi, S. Nabetani, and S. Kumagai , Complex nonlinear dynamics of the Hodgkin–Huxleyequations induced by time scale changes , Biological Cybernetics, 85 (2001), pp. 51–64.[10]
N. Fenichel , Geometric singular perturbation theory for ordinary differential equations , Jour-nal of Differential Equations, 31 (1979), pp. 53–98.[11]
J. Guckenheimer and I. Lizarraga , Shilnikov homoclinic bifurcation of mixed-mode oscil-lations , SIAM Journal on Applied Dynamical Systems, 14 (2015), pp. 764–786.[12]
P. Kaklamanos, N. Popovi´c, and K. U. Kristiansen , Geometric singular perturbationanalysis of the multiple-timescale Hodgkin-Huxley equations , in preparation, (2020).[13]
M. T. Koper , Bifurcations of mixed-mode oscillations in a three-variable autonomous Vander Pol-Duffing model with a cross-shaped phase diagram , Physica D: Nonlinear Phenomena,80 (1995), pp. 72–94.[14]
M. Krupa, N. Popovi´c, and N. Kopell , Mixed-mode oscillations in three time-scale sys-tems: a prototypical example , SIAM Journal on Applied Dynamical Systems, 7 (2008), pp. 361–420.[15]
M. Krupa and P. Szmolyan , Extending geometric singular perturbation theory to nonhy-perbolic points—fold and canard points in two dimensions , SIAM Journal on MathematicalAnalysis, 33 (2001), pp. 286–314.[16]
M. Krupa and M. Wechselberger , Local analysis near a folded saddle-node singularity ,Journal of Differential Equations, 248 (2010), pp. 2841–2888.[17]
C. Kuehn , On decomposing mixed-mode oscillations and their return maps , Chaos: An Inter-disciplinary Journal of Nonlinear Science, 21 (2011), p. 033107.[18]
B. Letson, J. E. Rubin, and T. Vo , Analysis of interacting local oscillation mechanismsin three-timescale systems , SIAM Journal on Applied Mathematics, 77 (2017), pp. 1020–1046.[19]
P. Nan , Dynamical Systems Analysis of Biophysical Models with Multiple Timescales , PhDthesis, ResearchSpace@ Auckland, 2014.[20]
J. Rubin and M. Wechselberger , Giant squid-hidden canard: the 3D geometry of theHodgkin–Huxley model , Biological Cybernetics, 97 (2007), pp. 5–32.[21]
P. Szmolyan and M. Wechselberger , Canards in R , Journal of Differential Equations,177 (2001), pp. 419–453.[22] M. Wechselberger , Existence and bifurcation of canards in R in the case of a folded node ,SIAM Journal on Applied Dynamical Systems, 4 (2005), pp. 101–139.2 (a) k = − . λ = 1 .
5. (b) k = − . λ = 1 . k = − . λ = 1 .
5. (d) k = − . λ = 1 . k = − . λ = 1 .
5. (f) k = − . λ = 1 . k = − . λ = 1 .
5. (h) k = − . λ = 1 . Figure 12: Verification of the bifurcation diagram in Figure 11 for representative choices of k , with λ = 1 . k decreases, one observes a transition from (a) steady-state behaviour via (c)MMO trajectories with double epochs of SAOs and (e) single epochs of SAOs to (g) relaxationoscillation. The corresponding singular geometry in phase space is shown in panels (b), (d), (f),and (h), respectively.3 (a) δ = 0 . O ( √ ε ). (b) δ = 0 .
001 = O ( ε ). (c) δ = 0 . O ( ε ). Figure 13: Mixed-mode time series in the Koper model for ε = 0 .
01 fixed and varying δ : as δ decreases, one typically observes more LAOs between SAO segments; additionally, for theseparticular parameter values, the model seems to exhibit sector-delayed-Hopf-type dynamics [4], asis particularly apparent in panel (b).Figure 14: Numerical continuation of periodic orbits in the Koper model with auto-07p [8] for λ = 1 . ε = 0 . δ : one observes coexistence of multiple periodic orbits, as evidenced by theoverlap between the corresponding k -intervals.4Figure 15: The critical manifold M of the three-dimensional, three-timescale Hodgkin-Huxleymodel, Equation (51).5 (a) I = 23. (b) I = 23.(c) I = 25 .
6. (d) I = 25 . I = 26 .
25. (f) I = 26 . I = 27. (h) I = 27. Figure 16: Bifurcations of MMOs in the reduced Hodgkin-Huxley model, Equation (51): as I isincreased, one observes a transition from (a) MMOs with double epochs of SAOs via (c) exoticMMOs and (e) those with single epochs of SAOs to (f) relaxation oscillation. The correspond-ing trajectories, projected onto the ( v, hv, h