Mixed-mode oscillations in a three-timescale coupled Morris–Lecar system

Mixed-mode oscillations (MMOs) are complex oscillatory behaviors of multiple-timescale dynamical systems in which there is an alternation of large-amplitude and small-amplitude oscillations. It is well known that MMOs in two-timescale systems can arise either from a canard mechanism associated with folded node singularities or a delayed Andronov–Hopf bifurcation (DHB) of the fast subsystem. While MMOs in two-timescale systems have been extensively studied, less is known regarding MMOs emerging in three-timescale systems. In this work, we examine the mechanisms of MMOs in coupled Morris–Lecar neurons with three distinct timescales. We investigate two kinds of MMOs occurring in the presence of a singularity known as canard-delayed-Hopf (CDH) and in cases where CDH is absent. In both cases, we examine how features and mechanisms of MMOs vary with respect to variations in timescales. Our analysis reveals that MMOs supported by CDH demonstrate significantly stronger robustness than those in its absence. Moreover, we show that the mere presence of CDH does not guarantee the occurrence of MMOs. This work yields important insights into conditions under which the two separate mechanisms in two-timescale context, canard and DHB, can interact in a three-timescale setting and produce more robust MMOs, particularly against timescale variations.


I. INTRODUCTION
Mixed mode oscillations (MMOs) are frequently perceived in the dynamical systems involving multiple timescales 11 ; these are complex oscillatory dynamics characterized by the concatenation of small-amplitude a) ngocanh-phan@uiowa.edub) yangyangwang@brandeis.eduoscillations (SAOs) and large-amplitude excursions in each periodic cycle.
Theoretical analysis of MMOs in systems with two distinct timescales has been well developed with the implementation of the geometric singular perturbation theory (GSPT) 13 ; see 11 for review.Two common mechanisms leading to the occurrence of MMOs in multiple timescale problems are canard dynamics associated with the twisting of slow manifolds due to folded singularities 40,49 and a slow passage through the delayed Andronov-Hopf bifurcation (DHB) of the fast subsystem 2,16,32,33 .While in two-timescale settings, these two mechanisms remain separated, they can coexist and interact in three-timescale regime 9,29,41,43 .
Compared with the extremely well-studied MMOs in two-timescale problems, the theory of MMOs in the three-timescale settings has been less well-developed.Traditionally, three-timescale problems are simplified to two-timescale problems which is the natural setting for geometric singular perturbation theory 3 .However, many real-world systems have more than two timescales 6,9,18,35,43,[46][47][48] .It has also been established that a two-timescale decomposition fails to capture certain aspects of the system's dynamics 31 .Therefore, classifying three timescales into two groups is not a sufficient approach for modelling and analysis.
In this paper, we contribute to the investigation of MMOs in three-timescale settings by considering a model of 4-dimensional coupled Morris-Lecar neurons 30,36 that was introduced by 31 .The model equations are given by with Table I lists the parameter values for the model chosen to ensure that (2) exhibits three distinct timescales where V 1 is fast, w 1 , V 2 are slow and w 2 is superslow.In a more biologically realistic model for calcium and voltage interactions, V 1 might represent membrane potential, while V 2 might represent intracellular calcium concentration with appropriate adjustments to parameter units and functional terms (see, e.g., 46,47 ).For the physiological description of functions in (2) and (3), we refer readers to 31 for details.
In the absence of coupling g syn = 0, (V 1 , w 1 ) is excitable with an attracting critical point at relatively low V 1 value, whereas (V 2 , w 2 ) is oscillatory with an attracting limit cycle independent of the value of g syn and the dynamics of (V 1 , w 1 ).To analyze the three-timescale coupled Morris-Lecar neurons, Ref. 31 extended two approaches previously developed in the context of GSPT for the analysis of two-timescale systems to the threetimescale setting and showed these two approaches complemented each other nicely.By varying g syn in system (2), the authors identified various solution features that truly require three timescales, thus demonstrating the functional relevance of three timescales in the model.While system (2) exhibits both the fast subsystem Hopf and folded nodes that can support MMOs, MMOs were not observed within the parameter regime examined by Ref. 31 (e.g., g syn = 4.1 and 5.1 in Figure 1).The goal of this work is the analysis of MMOs and their robustness in three-timescale systems by focusing on a coupled Morris-Lecar system (2).Based on our simulations, we have selected two MMO solutions on which to focus our analysis.Specifically, we consider g syn = 4.3 and g syn = 4.4 (with the unit of mS/cm 2 ), as highlighted in Figure 1.To provide further insight into our choice of g syn values, we perform a bifurcation analysis to explore the effect of g syn on a singularity called canard-delayed-Hopf (CDH) that was firstly introduced by 29 (see Figure 2, blue).As noted above, this singularity plays a crucial role in organizing MMOs within the three-timescale setting.We show in Sections II that the full system (2) may exhibit two CDH points -one at larger V i values, i ∈ {1, 2} (denoted as upper CDH) and the other at smaller V i (denoted as lower CDH, see Figure 4).Similarly, (2) may exhibit an upper DHB and a lower DHB.However, we demonstrate in Sections III and IV that only the upper CDH or DHB can support MMOs.
In Figure 2, we plot the bifurcation curves of the upper CDH and the upper DHB in the (g syn , V 2 ) plane, both of which approach vertical asymptotes as g syn decreases.
Our first choice of g syn = 4.3 represents g syn values between the two asymptotes, at which the upper DHB exists but there is no upper CDH.In contrast, g syn = 4.4 lies to the right of the CDH asymptote and serves as a representative scenario where both the upper CDH and DHB exist.When g syn < 4.2628 (e.g., g syn = 1 and 4.1 as considered in 31 ), the system (2) does not produce MMOs.
Although both CDH and DHB are present for g syn = 5.1, the absence of MMOs as shown in Figure 1 is due to the real nature of eigenvalues in a relevant subsystem during the transition from V 1 spikes to a V 1 plateau 31 .Our analysis suggests that the two MMO solutions at g syn = 4.3 and g syn = 4.4 arise from distinct mechanisms, resulting in remarkably different sensitivities to variations in timescales (i.e., varying C 1 and φ 2 ), as illustrated in Figure 3. Utilizing the extended GSPT 13,31 , we show that there is no interaction of different MMO mechanisms due to the lack of a nearby CDH singularity when g syn = 4.3 (see Figure 2).Instead, the MMOs at g syn = 4.3 solely depend on the delayed Hopf mechanism and are sensitive to variations in timescales (Figure 3A).In contrast, there exists a CDH in the middle of the SAOs when g syn = 4.4 as discussed above.We demonstrate that this CDH allows the fast subsystem Hopf and a canard point to coexist and interact to comodulate properties of the local oscillatory behavior, resulting in MMOs with significantly stronger robustness than g syn = 4.3 (Figure 3B).In summary, our findings reveal that MMOs near a CDH exhibit stronger robustness compared to those governed by a single mechanism, and that the CDH singularity is a determining factor in whether the two distinct MMO mechanisms can interact or not.However, it is essential to note that not all CDH singularities can support local MMOs.Specifically, we demonstrate that the lower CDH does not support MMOs.
Our work is novel in two main aspects.First, to the best of our knowledge, our study is the first to investigate the geometric conditions that lead to robust occurrences of MMOs in three-timescale systems.It is worth noting that while Ref. 20 also considered the robustness of MMOs in a three-timescale system, their focus was specifically on MMOs with double epochs of SAOs.Second, we discovered that the CDH singularities do not always enable the two MMO mechanisms to interact and produce MMOs.This is different from past studies 29,43 where the CDH always leads to occurrence of MMOs.From analyzing system (2), we found that CDH singularities that lie close to the actual fold of the superslow manifold (defined later by ( 14)) do not support MMOs regardless of perturbation sizes ε and δ.As the first step of our timescale decomposition approach, we perform a dimensional analysis of (2) to reveal the important timescales.This transforms (2) to the following three-timescale problem where ε = 0.1, δ = 0.053, t s is the slow dimensionless time variable, f 1 , f 2 , g 1 and g 2 are O(1) functions specified in (A1) in the Appendix A which include details of the nondimensionalization procedure.For simplicity, we did not rescale V 1 and V 2 in (4) as the scalings of voltage have no influence on the timescales.
We call system (4) that is described over the slow timescale the slow system in which V 1 evolves on a timescale of O(ε −1 ), (w 1 , V 2 ) on O(1) and w 2 on O(δ).Introducing a superslow time t ss = δt s yields an equivalent description of dynamics: which evolves on the superslow timescale and is called the superslow system.Alternatively, defining a fast time t f = t s /ε, we obtain the following fast system: which evolves on the fast timescale.
The paper is organized as follows.In Section II, we perform a geometric singular perturbation analysis of the 3-timescale problem (2) by treating ε as the principal perturbation parameter while keeping δ fixed, treating δ as the principal perturbation parameter while keeping ε fixed, and by treating ε and δ as two independent perturbation parameters.We review both mechanisms for MMOs and discuss their interaction at the double singular limit (ε, δ) → (0, 0).Notation, subsystems, construction of singular orbits at various singular limits and other preliminaries relating to the method of GSPT are all presented in Section II.In Section III, we investigate the mechanism and sensitivity of MMOs when g syn = 4.3 to variations in perturbation sizes ε or δ (i.e.varying C 1 and φ 2 in (2)).We explain the transitions between MMO and non-MMO dynamics as we vary one perturbation parameter while keeping the other fixed at its default value, as illustrated by the two lines in Figure 3A.Our analysis indicates that MMOs persist for δ ≤ O(ε).Increasing δ via increasing φ 2 or decreasing ε via decreasing C 1 to a degree where δ > O(ε) leads to multiple MMOs/non-MMOs transitions.MMOs are completely lost when δ ≥ O(ε 1 3 ).While Figure 3A also shows transitions occurring when both ε and δ are relatively large, the analysis of these transitions is beyond the standard GSPT and falls outside the scope of this paper.In Section IV, we uncover the dynamic mechanism underlying MMOs from (2) when g syn = 4.4.In this case, the existence of a CDH in the middle of the SAOs enables two different mechanisms to co-modulate properties of MMOs.We explain why MMOs organized by a CDH singularity as seen in the case of g syn = 4.4 exhibit remarkable robustness against variations in timescales (see Figure 3B).Our analysis demonstrate that when δ = O(ε), MMOs exhibit both canard and DHB features.Upon tuning δ ≥ O( √ ε), DHB-like features disappear and the canard mechanism dominates.Finally, we conclude in Section V with a discussion.

II. GEOMETRIC SINGULAR PERTURBATION ANALYSIS
In this section, we apply the extended geometric singularity perturbation analysis 13,31 to the three-timesale coupled Morris-Lecar system (4) by treating ε as the only singular perturbation parameter 40,49 , treating δ as the only singular perturbation parameter 2,16,32,33 , and finally treating ε and δ as two independent singular perturbation parameters 29,31,43 .
Although the detailed GSPT analysis and derivation of subsystems have been previously presented in 31 , we provide a brief overview in this paper for the sake of completeness.However, the focus of our current work is on the investigation of MMOs, which is distinct from the emphasis of 31 .Specifically, we concentrate on reviewing and discussing the canard mechanism in subsection II B, delayed Hopf bifurcation in subsection II C, and their interactions in subsection II D.
Fixing δ > 0 and taking ε → 0 in the fast system (6) yields the one-dimensional (1D) fast layer problem, a system that describes the dynamics of the fast variable, V 1 , for fixed values of the other variables, The set of equilibrium points of the fast layer problem is called the critical manifold and is denoted as M S : Although M S is a three-dimensional (3D) manifold in R 4 space, it does not depend on w 2 .We can solve and can therefore represent M S as for a function F 1 .It is well known that, for sufficiently small ε > 0, normally hyperbolic parts of M S each perturb to a locally invariant manifold called a slow manifold, on which w 1 is given by an O(ε)-perturbation of F 1 13 ; we simply use M S as a convenient numerical approximation of these slow manifolds.
M S is a 3D folded manifold with two-dimensional (2D) fold surface, L s , given by or equivalently where F 1V1 and f 1V1 denote the partial derivatives of F 1 and f 1 with respect to V 1 .The fold surface divides the critical manifold M S into attracting upper (M U S ) and lower (M L S ) branches where F 1V1 < 0 and repelling middle branch (M M S ) where F 1V1 > 0. Taking the same limit, i.e., ε → 0 with δ > 0, in the slow system (4) yields the 3D slow reduced problem, a system that describes the dynamics of w 1 , V 2 , w 2 along M S , where f 1 = 0.
Alternatively, fixing ε > 0 and taking δ → 0 in the slow system (4) yields the 3D slow layer problem in the form where the superslow variable w 2 is a parameter.The set of equilibrium points of the slow layer problem ( 13) is defined to be the superslow manifold and is denoted as Similarly to M S , the normally hyperbolic parts of M SS perturb to nearly locally invariant manifolds for δ sufficiently small.Later in subsection II C, we will discuss the bifurcations of the slow layer problem (13), i.e., nonhyperbolic regions on M SS where Fenichel's theory (GSPT) breaks down.
Taking the same singular limit in the superslow system (5) leads to the 1D superslow reduced problem where f 1 = g 1 = f 2 = 0.The superslow motions of trajectories of (15) are slaved to M SS until nonhyperbolic points are reached.
Both the slow reduced problem (12) and the slow layer problem (13) still include two distinct timescales.Further taking the limit δ → 0 in (12) or taking the limit ε → 0 in (13) yields the same slow reduced layer problem, (16)   which describes the slow motion along M S and the superslow variable w 2 is fixed as a constant.
It follows that the double singular limits lead to three subsystems: the fast layer problem (7), the slow reduced layer problem (16) and the superslow reduced problem (15).In addition to the naturally expected fast/slow transitions and slow/superslow transitions, transitions directly from superslow to fast dynamics and from superslow to fast-slow relaxation oscillations have also been observed in 31 .

B. Slow Reduced Problem and Canard Dynamics
To investigate canard dynamics, we project the slow reduced problem (12) onto (V 1 , V 2 , w 2 ) to obtain a complete description of the dynamics along M S .To this end, we differentiate the graph representation of M S given by where F 1V2 := ∂F1 ∂V2 .Note that the reduced system ( 17) is singular at the fold surfaces L s (10).Nonetheless, this singular term can be removed by a time rescaling t s = −F 1V1 t d and we obtain the following desingularized system (18) We observe that the desingularized system (18) is equivalent to (17) on the attracting branch, i.e, for F 1V1 < 0, but has the opposite orientation on the repelling branch, i.e, for F 1V1 > 0.
The desingularized system (18) has two kinds of singularities: ordinary and folded singularities.The ordinary singularities are the true equilibria of the full system (2), which is defined by (19) For the chosen parameter set in Table I, E always lies on the repelling branch of M S and hence is unstable.In contrast to the ordinary singularities, the folded singularities are not equilibria of the full system.They lie on one-dimensional curves along the fold surface L s defined by (20) Folded singularities are special points that allow trajectories of (17) to cross the fold L s with nonzero speed.Such solutions are called singular canards 40 .Note that when projecting to (V 1 , w 1 , V 2 )-space, the condition F = 0 is redundant and the fold surfaces L s become curves that overlap with the folded singularity curves M.
Since there is a curve of folded singularities, the Jacobian of ( 18) evaluated along M (denoted as J D , see Appendix B) always has a zero eigenvalue and the eigenvector corresponding to this zero eigenvalue is tangent to the curve of folded singularities.Generically, the other two eigenvalues (λ w , λ s ) where |λ w | < |λ s | have nonzero real part and are used to classify the folded singularities.Folded singularities with two real eigenvalues with the same sign (resp., with opposite signs) are called folded nodes (resp., folded saddles).Those with complex eigenvalues are called folded foci, which does not produce canard dynamics.
In the stable folded node case, we have strong and weak eigenvalues 0 > λ w > λ s .The singular strong canard is the unique solution corresponding to the strong stable manifold tangent to the strong eigendirection.For each folded node, the corresponding strong canard and the fold surface L s form a two dimensional trapping region (the funnel) on the attracting branch of M S such that all solutions in the funnel converge to that folded node.The funnel family of all folded nodes of M and the fold surface L s then form a three dimensional funnel volume.Trajectories that land inside the funnel volume will be drawn into one of the folded nodes, passing through the fold surface from an attracting M S to a repelling M S due to a cancellation of a simple zero in (17), and such solutions are so-called singular canards.

Folded saddle node (FSN)
In (18), a degenerate singularity arises when a second eigenvalue, λ w = 0, becomes zero.This folded singularity is referred to as a "folded saddle node" (FSN), and is characterized by the condition ).The yellow curve consists of mostly folded foci points and small segments of other singularities (e.g., folded node, folded saddle) that are barely visible and hence are not displayed here.In the top two panels (A, B) when δ = 0, an FSN 1 point (blue star) is O(δ) close to a CDH singularity (blue diamond), whereas an FSN 2 (cyan star) is far away from any CDH.In the lower two panels (C, D) at the singular limit δ = 0, the FSN 1 point becomes a CDH singularity (blue star overlapping with blue diamond).The center subspace of an FSN 1 (resp., an FSN 2 ) is denoted by a pink plane (resp., a yellow plane).It follows that the center manifolds of both FSN 1 and FSN 2 are transverse to Ls.
where Q and P are defined in (B2) in Appendix B, which contains a detailed derivation of the FSN condition.Similar to 43 , we demonstrate in the appendix that our system can exhibit an FSN in two different ways: either where K 1 , K 2 and K 3 are defined in Appendix B. We discuss below in Remark II.1 that both FSN 1 and FSN 2 are novel types of FSN 29 .
Remark II.1 In our parameter regime, the ordinary singularity point lies in the middle branch of critical manifold and is not involved in any bifurcations of the folded singularities.Hence, the FSN singularities (21) are neither of type II nor type III 27,29,37,45 .While both FSN 1 and FSN 2 are saddle-node bifurcations of folded singularities, the corresponding center manifolds of these FSN singularities are transverse to the fold surface at the singularties (see Figure 4C and D, yellow and pink planes).
We follow 29 to denote both FSN 1 and FSN 2 as novel types of FSN.
Remark II.2The condition (22) suggests that an FSN 1 is O(δ) close to the intersection point of the superslow manifold M ss and the fold surface L s , which was defined as the canard-delayed-Hopf (CDH) singularity in 29 .In contrast, FSN 2 is always far away from a CDH. Figure 4 shows the positions of FSN 1 (blue star), FSN 2 (cyan star) and CDH (blue diamond) in (V 1 , V 2 , w 2 )-space, for δ = 0 (top panels) and the singular limit δ = 0 (bottom panels).Recall the bifurcation of the upper CDH with respect to g syn is shown in Figure 2, which explains why there is no upper CDH for g syn = 4.3.It is worth noting that a CDH point of (2) is always a folded singularity because the critical manifold M S does not depend on the superslow variable w 2 .

C. Slow Layer Problem and Delayed Hopf Bifurcations
In this subsection, we turn to the slow layer problem (13) resulting from the δ → 0 singular limit, which exhibits delayed Hopf bifurcations that allow for interesting dynamics.
Let J SL denote the Jacobian matrix of ( 13) evaluated along the superslow manifold M SS , which is given by where the nonzero entries denote partial derivatives.The eigenvalues of J SL are given by f 2V2 and the eigenvalues of .
Thus, the Hopf bifurcation points on M SS are given by tr(J) = 1 ε f 1V1 + g 1w1 = 0 and det J > 0. The former defining condition can be rewritten as Remark II.3It follows from ( 25) that an M H SS is O(ε) close to the intersection of M SS and L s , i.e., a CDH singularity.The subsystem HB bifurcation M H SS is also known as delayed Hopf bifurcation (DHB).
The isolated fold bifurcation points on M SS are located by letting det J SL = 0.That is, The fold points L ss that satisfy the former condition (denoted as L 1 ss ) are the folds of the V 2 -nullcline (see Figure 5, green circle and green triangle), which correspond to the transition between superslow dynamics along M SS and slow jumps.L ss given by the latter condition (denoted as L 2 ss ) corresponds to the actual fold of M SS when projected to (V 1 , w 1 , V 2 )-space.Since g 1w1 < 0, f 1w1 < 0 and g 1V1 > 0, it follows that L 2 ss lies on the middle branch of M S (f 1V1 = f 1w1 g 1V1 /g 1w1 > 0) and hence will not play a role in dynamics.At the double singular limit (ε, δ) → (0, 0), the L 2 ss fold point of M SS will occur at the fold curve of the critical manifold and become a CDH.This can be shown by analyzing the slow reduced layer problem ( 16) obtained from the double singular limits.

D. Interaction between canard and delayed hopf
To investigate the interaction between the canard and the delayed Hopf mechanisms in the double limit case (ε → 0, δ → 0), we need to examine the slow reduced layer problem (16).The corresponding desingularized system is given by where F and G are defined in (18).Note that ( 27) is the δ → 0 limit of the desingularized system (18) from subsection II B. The folded singularities of ( 27) are exactly the same as M given by ( 20), whereas the ordinary singularities of ( 27) are relaxed to be M SS .The FSN condition at the double singular limit can be obtained from letting δ → 0 in the FSN 1 condition (22) or the FSN 2 condition ( 23).This implies that, at the double singular limit, an FSN 1 becomes a CDH singularity (see Figure 4C and D) whereas an FSN 2 singularity is characterized by According to Remarks II.2 and II.3, an FSN 1 singularity from the ε viewpoint converges to a CDH as δ → 0 and a DHB M H SS from the δ viewpoint converges to a CDH as ε → 0. It is natural to expect that a CDH singularity point should serve as the interplay between the canard dynamics and the delayed Hopf bifurcation to produce MMOs, as seen in 29 .However, this is not always the case.Specifically, we find that while the CDH on the upper fold surface L s (see Figure 4B and D) supports MMOs with a high level of robustness due to the coexistence and interaction of two distinct MMO mechanisms, no MMO dynamcis were observed near the lower CDH.It is worth highlighting that both the upper and lower CDH points in our system are of the same type as the CDH investigated in 29 .
The CDH points in system (2) are FSN 1 singularities at the double singular limit.We prove in Appendix C that the CDH singularity in (2) is a novel type of saddlenode bifurcation of folded singularities as described in 29 , with the center manifold of the CDH transverse to the fold L s of the critical manifold.This is further confirmed in Figure 4C and D, as discussed in Remark II.1.
In the case of g syn = 4.3 (see the left panels of Figure 4), there is no upper CDH.As a result, there is no coexistence and interaction of canard and delayed Hopf mechanisms, leading to MMOs that are sensitive to variation of timescales (see Section III).For g syn = 4.4, an upper CDH exists.We show in Section IV that this CDH serves as an organizing center for the local smallamplitude oscillatory dynamics, which results in robust MMOs through the interplay of the DHB and canard mechanisms.In both cases, there exist CDH points on lower L s .However, MMOs are not observed in the neighborhood of any lower CDH (see Section III for more discussions).

E. Singular orbit construction
To understand the dynamics of (2) using GSPT 13 , we need to construct singular periodic orbits by concatenating solution segments of singular limit systems.According to GSPT, such a singular oscillation will generically perturb to a periodic solution of the full problem as we move away from the singular limit.
We now use our analysis from previous subsections to construct singular approximations of (2) in order to understand the full system trajectory Γ (ε,δ) .We let Γ F (0,δ) and Γ S (0,δ) denote trajectories of the fast layer problem ( 7) and the slow reduced problem ( 12) obtained from the ε → 0 singular limit.Let Γ S (ε,0) and Γ SS (ε,0) denote trajectories of the slow layer problem ( 13) and the superslow reduced problem (15) obtained from the δ → 0 singular limit.The process of constructing singular orbits at the double singular limits is more complicated since there are more than two singular limit systems.We let Γ x (0,0) , x ∈ {F, S, SS} denote the fast, slow and superslow flows at the double singular limits.FIG. 5. Projection of the singular orbit (green) and the solution trajectory Γ (ε,δ) (black) of the full system (2) onto the phase plane of (V2, w2) system with parameters given in Table I.The red curve is the V2-nullcline, which is the projection of the superslow manifold MSS.The dark and light green symbols mark the key transitions between the slow and superslow sections of the singular orbit and the perturbed oscillation, respectively: the star and square indicate the transitions from the slow to the superslow motions, and the circle and triangle mark the transitions from superslow to slow sections at the fold of the V2-nullcline.The circled numbers indicate four phases of the oscillations: superslow excursions along MSS during 1 ○ and 3 ○ and slow jumps at the fold of MSS during 2 ○ and 4 ○.
We start by showing the singular orbit construction for the (V 2 , w 2 )-subsystem (Figure 5), which is a relaxation oscillation that is independent of g syn , (V 1 , w 1 ), and ε.The singular orbit (Figure 5, green trajectory) consists of a superslow excursion along the left branch of V 2 -nullcline (Γ SS (ε,0) , phase 1 ○), a slow jump at the lower fold of V 2nullcline up to its right branch (Γ S (ε,0) , phase 2 ○), a superslow excursion through the active phase (Γ SS (ε,0) , phase 3 ○) and a slow jump back to its left branch (Γ S (ε,0) , phase 4 ○).Green symbols mark points at the key transition between the four different sections of the oscillation and will be used for later analysis.When projected onto (V 2 , w 2 )space, the V 2 -nullcline (red) corresponds to the superslow manifold M SS and the singular orbit Γ S (ε,0) ∪ Γ SS (ε,0) overlaps with the singular orbit from the double singular limits due to its independence on ε.For sufficiently small perturbation δ, Γ S (ε,0) ∪ Γ SS (ε,0) perturbs to the full system trajectory Γ (ε,δ) shown by the black curve.
Figure 6A illustrates the singular orbit Γ S (ε,0) ∪ Γ SS (ε,0) -space, together with the superslow manifold M SS (red solid curve M a SS : attracting; red dashed curve M r SS : repelling).In other panels, the critical manifold M S (blue surface) and its folds L s (blue curves) are also plotted.M S is separated into three sheets by the folds L s , in which the upper (M U S ) and lower (M L S ) branches are stable and the middle branch (M M S ) is unstable.
Starting from the green star in panel (A), the singular orbit is in phase 1 ○ and evolves along the lower branch of M SS under the superslow reduced problem (15).After hitting the DHB where the stability of M SS changes, the slow layer problem (13) takes over but with V 2 still evolving on a superslow timescale (see Figure 5, phase 1 ○).Thus, the singular orbit during the rest of this phase is a continuum of (V 1 , w 1 ) relaxation oscillations.As the evolution speed of V 2 changes from superslow to slow at the green circle (beginning of phase 2 ○), a few more spikes occur before the slow flow Γ S (ε,0) travels to the green square on M SS , at which phase 3 ○ begins.After that, the superslow reduced problem takes over until the singular orbit Γ SS (ε,0) reaches the upper DHB at which M SS becomes unstable.As such, a family of (V 1 , w 1 ) oscillations emerges as we observed during phase 1 ○.As phase 4 ○ begins at the green triangle, several additional V 1 spikes occur before the slow flow Γ S (ε,0) travels back to the green star, thus completing a full cycle.
Figure 6B shows the singular orbit at the double singular limits.It closely resembles the orbit in panel (A), with the notable exception that there is no longer a superslow segment along the upper M SS .Instead, we observe a continuum of V 1 spikes throughout phase 3 ○.This is because the upper M SS becomes unstable for all V 2 values as ε → 0, which will be discussed further in subsection III B 1. In (B), the singular orbit consists of Γ F (0,0) (triple arrow) that are fast V 1 jumps from L s , Γ S (0,0) (double arrow) which travels along stable branches of M S or the intersection of M S and the V 2 -nullcline, and Γ SS (0,0) (single arrow) that follows the stable branch of M SS .
While the V 2 value at the slow/superslow transition (green circle or green triangle) is uniquely determined, the corresponding (V 1 , w 1 ) values can assume arbitrary positions on the upper or lower sheet of M S for that fixed V 2 .Thus, infinitely many singular orbit segments can be constructed during phases 2 ○ and 4 ○, although we only plot one in panel (B) for clarity.For the same reason, the singular orbits in (A) are also not unique since there exist infinitely many ways to choose a starting position for Γ S (ε,0) during phases 2 ○ and 4 ○.Comparing (A) and (B) indicates that to obtain MMO dynamics in the perturbed system, ε cannot be too small.To explain the MMO dynamics for g syn = 4.3 in section III, we mainly make use of the singular orbit in panel (A) with 0 < ε ≪ 1.This orbit, aside from the upper superslow segment where the stability of M SS differs between the two panels, can be viewed as an O(ε) perturbation of the singular orbit in panel (B).Hence, when discussing the fast-slow oscillations in the full system, we still refer to different segments of them as being governed by (7)  which describes fast V 1 jumps and ( 16) which describes the slow motion along M S .
Lastly, the singular orbit for g syn = 4.4 at the double singular limits is shown in Figure 6D.There are two major differences between g syn = 4.3 and g syn = 4.4: Firstly, in contrast to the g syn = 4.3 case where the upper DHB vanishes, it now converges to the upper CDH at the dou-ble singular limits in (D).Secondly, Γ SS (0,0) in panel (D) follows the upper stable branch of M SS until reaching a saddle-node bifurcation at the green triangle where M SS becomes unstable.This results in a unique construction of singular orbit segment during phase 4 ○ due to its unique starting position.Nonetheless, the construction of singular orbits during phase 2 ○ remains non-unique as discussed in the g syn = 4.3 case.Later in section IV, we show how this singular orbit perturbs to MMOs at g syn = 4.4 for various perturbations.

III. ANALYSIS OF MMOS WHEN gsyn = 4.3 mS/cm 2
In this section, we study MMOs that arise when there is no CDH singularity in the middle of the smallamplitude oscillations (SAOs).We show that the only existing mechanism for MMOs at g syn = 4.3 is the delayedhopf bifurcation (DHB) (see subsection III A) and explain why the absence of an upper CDH leads to the sensitivity of MMOs to timescale variations (see subsections III B and III D).In particular, we explain the complex transitions between MMOs and non-MMOs due to changes of ε or δ in subsection III D (see Figure 3A).We also discuss why there is no MMOs near the lower CDH in subsection III C.

A. Relation of the trajectory to MS and MSS
The MMO solution of (2) for g syn = 4.3 from Figure 1 is projected onto the (V 1 , w 1 , V 2 )−space in Figure 7.The full system equilibrium (black circle) lies on M M S and is unstable.The stability of the superslow manifold M SS changes at the two DHBs that are subcritical (red circles).In particular, as V 2 decreases, the upper branch of M SS changes from stable-focus (with one negative real eigenvalue and a pair of complex-conjugate eigenvalues whose real parts are negative) to saddle-focus (one negative real eigenvalue and a pair of complex-conjugate eigenvalues whose real parts are positive).Note that the upper fold L s always lies above the upper branch of the superslow manifold M SS so there is no upper CDH, whereas the lower fold intersects the lower branch of M SS at a CDH singularity (blue diamond in the lower inset).This CDH is a folded focus and will become an FSN 1 as δ → 0 as discussed in §II D.Moreover, the nearby HB bifurcation (red circle in the lower inset) is O(ε) close to this CDH.
Figure 7 shows that the singular orbit from Figure 6A (green curve) is a suitable predictor of the full trajectory (black curve).For the sake of clarity, we choose to not display the entire singular orbit but instead focus on regions where small amplitude oscillations (SAOs) emerge.During phase 3 ○, the upper superslow segment Γ SS (ε,0) perturbs to Γ (ε,δ) which displays two SAOs around the stable branch of M SS .These SAOs soon transition to large-amplitude oscillations before crossing the DHB at the red circle to reach the unstable branch of M SS .As phase 4 ○ begins at the light green triangle, the slow jump down of V 2 brings the trajectory to a region of M L S where there is a nearby stable M SS that attracts the trajectory.From the green star, the trajectory follows M SS on the superslow timescale to the lower fold of M S (see Figure 7, lower inset), where it jumps up to M U S on the fast timescale under (7), which corresponds to the onset of spikes in V 1 .The spikes persist until reaching the green circle, at which phase 2 ○ begins and the slow jump up of V 2 brings the trajectory to the green square, completing a full cycle.
We claim that the only mechanism underlying the MMOs at g syn = 4.3 is the DHB mechanism.This is not surprising given that the upper Γ SS (ε,0) switches to a continuum of big spikes when the upper DHB vanishes at the singular limit ε = 0, as shown in Figure 6B.To understand why canard dynamics are not involved, we view the trajectory in (V 1 , V 2 , w 2 )-space (see Figure 8).Unlike Figure 7 where L s and curves of folded singularity (M) overlap, the (V 1 , V 2 , w 2 )-projection captures the structure of the fold surfaces L s and lets us examine the position of the solution trajectory relative to M. To illustrate how SAOs arise from the DHB mechanism, we look at the projection of the trajectory onto (V 2 , V 1 )-plane, which includes the periodic orbits born at the upper Hopf bifurcation (see Figure 9).The solution from Figure 7 projected onto (V1, V2, w2)-space.Also shown are the projections of the fold surface Ls (blue surface), superslow manifold MSS (red curve) and folded singularities M (green and yellow curves).The cyan star denotes FSN 2 defined in (23).Color codings of the folded singularity curves are the same as in Figure 4. Other symbols have the same meanings as in Figure 7.

Canard mechanism does not contribute to MMOs
In Figure 8, the solution of (2) for g syn = 4.3 is projected onto the space (V 1 , V 2 , w 2 ) with two separate fold surfaces L s (blue) and two branches of folded singularities.In the lower folded singularity curve, there is an FSN 2 (cyan star) separating the folded singularities that are mostly folded foci (yellow curve) and folded saddle (dashed green).The inset around the lower CDH (blue diamond) shows that the trajectory crosses the lower fold at a regular jump point and hence immediately jumps up to M U S .The upper folded singularity curve also has an FSN 2 (cyan star) marking the boundary between folded saddles (green dashed curve) and stable folded nodes (green solid curve).However, as shown in the upper inset, the trajectory crosses the upper L s at the normal jump points that are distant from folded nodes.Hence, the emergence of SAOs in the MMOs is not due to the canard mechanism.

MMOs arise from the delayed Hopf mechanism
To examine how the DHB mechanism engages in organizing the MMOs, we turn to the subsystems obtained by treating δ as the main singular perturbation parameter (see §II C). Figure 9A shows the projection of the solution, M SS and the bifurcation diagram of the slow layer problem (13) onto the (V 2 , V 1 )-plane.As discussed above, the upper Γ SS (ε,0) (green curve) perturbs to SAOs followed with large spikes (black curve).Starting from the green square, the SAOs gradually decrease in magnitude as the trajectory moves towards the upper DHB.After two such oscillations, the trajectory crosses the unstable periodic orbit branch, whose amplitude also decreases as it approaches the upper DHB.Upon crossing, the trajectory undergoes a sudden jump to the outside large-amplitude periodic orbit branch, giving rise to large spikes.
Recall at default parameters, δ ≈ 0.053.Figure 10 shows that decreasing the perturbation δ improves the singular limit predictions and increases the amount of time the full system trajectory Γ (ε,δ) spends near the attracting branch of upper M SS (denoted as M a SS ).In particular, panel (B) indicates the trajectory with a slightly smaller perturbation δ ≈ 0.042 is able to pass over the upper DHB (red circle) to the repelling side of M SS (denoted as M r SS ) and there is a delay in which the trajectory traces M r SS before it jumps away.However, we recognize that Γ (ε,δ) under the default perturbation δ ≈ 0.053 does not exhibit a similar delay (see Figure 9A), which is contrary to what one would expect in DHBinduced MMOs.While this might initially suggest that δ ≈ 0.053 is distant from the singular limit, we show below that this is not the case.Moreover, this specific value of δ aligns with the perturbation used in 31 , where GSPT analysis has been successfully employed to elucidate the dynamics across various g syn values.
It is worth emphasizing two interesting points: Firstly, the impact of decreasing δ on SAOs is non-monotonic.In certain cases, smaller perturbations can lead to fewer SAOs with even larger amplitudes.This effect is related to how the slow flow during phase 2 ○ approaches M a SS , a detailed analysis of which is provided in subsection III D. Therefore, the absence of a delay in Figure 9 does not imply a significant deviation from the singular limit.Rather, it is mainly due to the manner in which the trajectory approaches M a SS during phase 2 ○, resulting in small oscillations that cross the unstable inner periodic orbit branch before passing through the DHB.As discussed above (also see Figure 14), a slight reduction or increase in the perturbation size can induce a delay phenomenon.Secondly, the plateauing behavior of the trajectory after passing the Hopf bifurcation is somewhat different from what one would expect to see in a typical DHB fashion, which typically involves oscillations with diminishing and then increasing amplitude.This is because the variable V 2 switches from superslow to slow timescale at the green triangle shortly after passing the HB, and there is insufficient time for the trajectory to oscillate.Hence, the associated pattern after HB is plateauing, and the amount of time the trajectory spends near M r SS is significantly shorter than that near M a SS .
In the following three subsections, we demonstrate how the absence of the interaction between canard and DHB mechanisms, specifically due to the lack of an upper CDH, can result in the sensitivity of the MMOs to timescale variations.To achieve this, we first explore how changes in the singular perturbation parameters ε and δ can induce transitions between MMOs and non-MMOs by analyzing their impacts on the two MMO mechanisms -canard and DHB.Additionally, we provide an explanation for why the lower CDH singularity does not guarantee the occurrence of SAOs.When g syn = 4.3, the CDH singularity only exists on the lower fold surface of M S .We demonstrate that this leads to different effects of ε on the upper and lower DHB points, respectively (see Figure 11).We also examine the effect of δ on the lower FSN points (see Figure 12).

Effect of ε on DHB
The effects of ε on the DHBs M H SS are summarized by the 2-parameter bifurcation diagrams of (13) projected onto (V 2 , C 1 )-space (see Figure 11).As C 1 decreases (or equivalently, as ε decreases), the upper Hopf moves to larger V 2 values and eventually vanishes for ε small enough (see Figure 11A).This explains why the upper singular orbit Γ SS (ε,0) in Figure 6A for ε = 0 switches to a continuum of spikes in Figure 6B as ε → 0. On the other hand, since there is a CDH on the lower L s , the lower Hopf will converge to that CDH as ε → 0 (see Fig-

Effect of δ on FSN points
There is no CDH or FSN 1 on the upper fold, hence, we only examine the effect of δ on FSN singularities on the lower L s .Figure 12A shows as δ → 0 (or equivalently, φ 2 → 0), the FSN 1 singularity converges to the lower CDH as demonstrated in our analysis (see (28)).On the other hand, the FSN 2 singularity is significantly distant from the CDH (see Figure 12B    Before examining the effects of varying ε and δ on MMO dynamics based on their impact on DHB and FSN points, we discuss briefly in this subsection why there are no SAOs near the lower CDH. As discussed above (also see Remarks II.2 and II.3), the lower CDH is O(δ) close to an FSN 1 and O(ε) close to a DHB.One would naturally expect to observe SAOs arising from Canard and/or DHB mechanisms near this lower CDH.Nonetheless, there is no SAOs near the lower fold regardless of ε and δ values.Below we explain why none of the two mechanisms produces SAOs.
As discussed earlier, the reason why there is no canardinduced SAOs when ε and δ are at their default values is because the trajectory crosses the fold surface at a regular jump point near which the folded singularities including the CDH are folded focus.This remains to be the case as ε varies or as δ increases.While as δ decreases the trajectory will follow M SS more closely and hence cross the fold surface somewhere near or at a folded node, that folded node is very close to an FSN where canard theory breaks down.As a result, the existence of canard solutions for smaller δ is not guaranteed.
On the other hand, with default ε and δ values, the trajectory jumps away from the lower fold before reaching the Hopf bifurcation and hence there is no DHB-induced small oscillations.Increasing δ will make the DHB less relevant, whereas increasing ε will move the DHB further away from the lower L s and eventually vanish upon crossing the actual fold of M SS (see Figure 11B).Thus, we do not expect to detect MMOs with increased ε or δ.As ε or δ decreases, the trajectory should pass closer to the lower DHB point.This is because decreasing ε moves the Hopf bifurcation closer to the CDH singularity (see Figure 11B) and reducing δ pushes the trajectory to travel along M SS more closely.However, the reason that no SAOs are induced by the passage through the lower HB is that this HB is relatively close to the actual fold of M SS , i.e., close to a double zero eigenvalue at a BT bifurcation of the (V 1 , w 1 , V 2 ) subsystem (see Figure 9A and Figure 11B).As a result, the branch of unstable small-amplitude periodic orbits born at the lower HB is almost invisible (see the inset of Figure 9A) and there is only a small region of M SS along which the Jacobian matrix J SL (24) of the slow layer problem (13) has complex eigenvalues.Figure 9B shows the real part of the first two eigenvalues λ of ( 24) evaluated along M SS , excluding the third eigenvalue given by f 2V2 .The eigenvalues are real when there are two branches of curves for Re λ and become complex when the curves coalesce and become a single branch.In panel (B), the eigenvalues on the stable lower branch of M SS are initially real and negative.That is, the trajectory approaches the attracting M SS along stable nodes of the slow layer problem.As the superslow flow brings the trajectory towards the Hopf bifurcation, the eigenvalues become complex.However, this region of complex eigenvalues is short and λ becomes real again shortly after.As a result, the trajectory has insufficient time to oscillate and we do not observe any small-amplitude oscillations before the trajectory jumps up to the outer periodic orbit branch.

D. Effects of varying ε and δ on MMOs
Recall when g syn = 4.3, the only mechanism available for MMOs is the DHB mechanism.While there exist folded node singularities, they do not play any significant role in generating MMOs.In this subsection, we explore the effects of ε and δ on the dynamics of the full system (4) by mainly examining their effects on the upper DHB around which SAOs are observed.As before, we will vary ε by changing C 1 and vary δ by changing φ 2 .Recall that increasing (resp., decreasing) C 1 or ε slows down (resp., speeds up) the fast variable V 1 , whereas increasing (resp., decreasing) φ 2 or δ speeds up (resp., slows down) the superslow variable w 2 .Other (slow) variables are not affected.
Figure 3A summarizes the effects of (C 1 , φ 2 ) on MMOs when g syn = 4.3.Our findings suggest that MMOs with only the DHB mechanism are robust to changes that slow down either the fast variable or the superslow variable, but they are vulnerable to perturbations that speed up either timescale to a degree where δ is approximately greater than O(ε).Specifically, we observe the following: • For fixed δ = 0.053 at φ 2 = 0.001, slowing down the fast variable V 1 by increasing C 1 from its default value C 1 = 8 (Figure 3A, vertical line above the red star) preserves the MMOs (also see Figure 13A, from top to bottom row).Moreover, we observe more characteristics of DHBs in the SAOs as ε increases.
• For fixed δ at φ 2 = 0.001, speeding up the fast variable V 1 via decreasing C 1 leads to a total of 3 transitions between MMOs and non-MMOs (see Figure 13B).These transitions correspond to the 3 crossings between the vertical black line and the yellow/blue boundary in Figure 3.
• For fixed ε = 0.1 at C 1 = 8, slowing down the superslow variable w 2 by decreasing φ 2 from the default value φ 2 = 0.001 (Figure 3A, red star) preserves MMOs.However, the number of small oscillations in the MMOs does not exhibit a simple monotonic increase or decrease, but rather shows an alternating pattern of increase and decrease as φ 2 decreases.Additionally, the amplitudes of the small oscillations also display a similar nonmonotonic behavior (see Figure 14A).
• For fixed ε at C 1 = 8, speeding up the superslow variable w 2 by increasing φ 2 leads to a total of five transitions between MMOs and non-MMOs (see Figure 14B).These transitions correspond to the five crossings between the horizontal black line and the yellow/blue boundary in Figure 3A.
Next, we discuss the above four scenarios separately.

MMOs are robust to increasing C1
Increasing C 1 slows down the fast variable V 1 and hence moves the three timescale (1F, 2S, 1SS) problem closer to (3S, 1SS) separation.As a result, the critical manifold M S and the folded singularities become less relevant with increased C 1 .Nonetheless, this does not affect the existence of MMOs, as they arise from the upper Hopf bifurcation.As C 1 is increased, the upper HB moves to smaller V 2 values (see Figure 11A).This change allows the trajectory to travel a longer distance along the stable part of M SS during phase 3 ○ and generate more SAOs, as shown in Figure 15.To simplify the presentation, we omit singular orbits and focus only on M SS and bifurcation diagrams that are essential for organizing SAOs in the full model, as demonstrated earlier.
Recall that the small oscillations occurring during phase 3 ○ switch to large spikes upon crossing the inner unstable periodic orbit branch that is born at the upper HB (see Figure 9A).Increasing C 1 causes the upper HB in (V 2 , V 1 )-space (Figure 15, red circle) to move to the left and become further away from the maximum of V 2 (green square) which remains unchanged.As a result, the trajectory with larger C 1 begins small oscillations with decaying amplitude at a greater distance from the Hopf bifurcation point.This greater distance results in the trajectory crossing the inner unstable periodic orbit at a smaller value of V 2 (compare Figures 9 and 15).As C 1 is increased to 11, the trajectory passes over the HB to M r SS and experiences a delay before jumping away (see Figure 15, right panel).As explained earlier, the absence of small oscillations with growing amplitudes after the upper HB is due to the slow jump down of V 2 at the green triangle.The initial decrease of C 1 (e.g., from C 1 = 8 to C 1 = 7, see Figure 16A) results in the loss of SAOs and thus a transition to non-MMOs due to the opposite effect of the mechanism discussed in the case of increasing C 1 .Interestingly, a further decrease of C 1 to a range of C 1 ∈ (1.46, 4.2), which causes the upper HB to cross the green square and eventually vanish, results in the recovery of MMOs characterized by small oscillations with increasing amplitude (see Figure 16B and C).This is because, for C 1 ∈ (1.46, 4.2) and V 2 near the green square, (13) exhibits either unstable periodic orbits with negligible amplitudes or saddle foci equilibrium characterized by a negative real eigenvalue and complex eigenvalues with positive real parts.As a result, the SAOs during phase 3 ○ grow in amplitude as the trajectory spirals away from M SS .When C 1 is reduced to be below 1.46, the voltage V 1 exhibits rapid spikes that occur immediately after the green square, consistent with the singular orbits shown in Figure 6B and C. As a result, there is no more MMOs (see Figure 16D, C 1 = 1.4).
To better understand why the SAOs for C 1 ∈ (1.46, 4.2) grow in amplitude (see Figure 16)B and C), we project the trajectory when 17).The blue triangle denotes a saddle focus of the (V 1 , w 1 , V 2 ) subsystem near the green square.During phase 2 ○ from the green circle to the green square, the trajectory travels towards the saddlefocus (blue triangle) along its stable manifold (magenta curve) on the slow timescale.After phase 3 ○ begins at the green square, SAOs grow in amplitude as the trajectory moves upwards and spirals away from the equilibrium curve along its unstable manifold (not shown).Similar dynamical behaviors have also been observed near a subcritical Hopf-homoclinic bifurcation 14 and a singular Hopf bifurcation in two-timescale settings 11 .
FIG. 17.The solution (black curve) from Figure 16C projected onto (V1, w1, V2) space, together with MSS (red curve).The blue triangle denotes a saddle-focus equilibrium on MSS with maximum V2, which has one negative real eigenvalue and a pair of complex eigenvalues with positive real parts (0.0027 ± 0.31i).A stable manifold branch associated with the saddle-focus at the blue triangle is denoted by the magenta curve.

Decreasing φ2 preserves MMOs but causes non-monotonic effects on SAOs
Decreasing the perturbation parameter φ 2 (i.e., reducing δ) moves the system closer to the 3-slow/1-superslow splitting and manifests the DHB mechanism.This preserves the DHB-induced MMOs as expected.However, a non-monotonic effect on the small-amplitude oscillations is also observed, as shown in Figure 14A, where the amplitude and the number of SAOs exhibit an alternation between increase and decrease.
Smaller perturbations should cause the solution trajectory to follow M SS more closely and at a slower rate.Intuitively, one may expect that this leads to an increase in the number of SAOs and a decrease in their amplitudes.Indeed, we observe such changes as φ 2 decreases from 0.001 to 0.0008, as shown in Figure 14A (top three rows) and also in Figure 10 as we discussed before.However, to our surprise, we find that for φ 2 = 0.0007, the MMOs exhibit less SAOs with larger amplitudes than the SAOs at φ 2 = 0.0008 (see Figure 14A, the third and the fourth row).The number of SAOs increases and their amplitudes decrease again as φ 2 decreases from 0.0007 to 0.0006.This alternating pattern of changes in the number and amplitude of SAOs repeats as φ 2 continues to decrease (see Figure 14A).
Our analysis reveals that as φ 2 decreases, there will be more SAOs with smaller amplitudes if no additional big (full) spike is generated.However, if an additional full spike is gained during the process of decreasing φ 2 , the changes to the SAOs will be reversed; that is, there will be fewer SAOs and they will exhibit larger amplitudes.This is because the additional spike before SAOs can push the trajectory away from M SS at the beginning of phase 3 ○, leading to fewer SAOs with larger amplitudes.Hence, the amplitude and number of SAOs in the full system are not only determined by the size of the perturbation but also by how the flow approaches M a SS during phase 2 ○.As previously discussed, infinitely many singular orbit segments can be constructed during this phase.This leads to different ways for the full trajectory to reach M a SS under varying perturbations.In a sense, this alternating pattern of changes in SAOs occurs due to a spike-adding like mechanism.We refer the readers to Appendix D for a more detailed discussion on why decreasing φ 2 causes non-monotonic effects on SAOs.
When φ 2 increases, the superslow variable w 2 speeds up, moving the system closer to the 1-fast/3-slow splitting and making the DHB mechanism less relevant.Since there is no canard mechanism, MMOs should be eliminated for φ 2 large enough.Indeed, we observe a total of five transitions between MMOs and non-MMOs as φ 2 increases, and eventually, MMOs are lost for φ 2 > 0.007 (i.e., δ ≥ O(ε 1 3 )).The mechanism driving these MMOs/non-MMOs transitions over the increase of φ 2 is similar to the mechanism underlying the non-monotonic effects on SAOs when φ 2 is decreased.Specifically, if no additional big (full) spike before phase 3 ○ is lost with the increase of φ 2 , there will be fewer SAOs with larger amplitudes or no MMOs as one would naturally anticipate (e.g., when φ 2 increases from 0.0012 to 0.0016, see Figure 14B, top two rows).In contrast, if one full spike is lost during the process of increasing φ 2 , changes to the SAOs will be reversed such that there will be more SAOs with smaller amplitude (e.g., when φ 2 increases from 0.001 to 0.0012, see the top row in Figure 14A and B).Eventually, MMOs will be completely lost when φ 2 ≥ 0.008 for which the HB is no longer relevant.
For a more detailed discussion, we refer the readers to Appendix E.
IV. ANALYSIS OF MMOS WHEN gsyn = 4.4 mS/cm 2 This section explores MMOs that occur when an upper CDH singularity is present.In this scenario, we show in subsection IV A that both canard and DHB mechanisms coexist and interact to produce MMOs that exhibit significant robustness to timescale variations as shown in Figure 3B.We explain the robust occurrence of MMOs in subsection IV C and show that the two MMO mechanisms can be modulated by adjusting ε and δ.Specifically, increasing ε manifests the DHB-like characteristics, while an increase in δ leads to dominance of the canard mechanism.

A. Relation of the trajectory to MS and MSS
The solution of (2) for g syn = 4.4 in Figure 1 is projected onto the (V 1 , w 1 , V 2 )−space, together with the singular orbit from Figure 6D (green curve), critical manifold M S (blue surface), fold L s (blue curve) and the superslow manifold M SS (red curves) (see Figure 18).As the coupling strength g syn increases from 4.3 to 4.4, two new features regarding the upper M SS emerge.First, the stability of the upper M SS changes at a fold point L 1 ss (the green triangle) rather than a DHB.Specifically, as V 2 increases, the equilibrium along the upper M SS switches from unstable saddle-focus (characterized by one positive real eigenvalue and a pair of complex eigenvalues whose real parts are negative) to stable focus (with one negative real eigenvalue and a pair of complex eigenvalues whose real parts are negative).Second and more importantly, the upper fold L s now intersects the upper branch of M SS at a CDH singularity (the blue diamond), which is a folded node at the default parameter values given in Table I.As proved in subsection II D, this CDH point is located O(δ) close to a folded saddle-node singularity FSN 1 (upper blue star) and O(ε) close to a HB (upper red circle).This is further confirmed in Figure 19, which shows that the upper FSN 1 and HB point converge to the same CDH point on the upper L s in the double singular limit (ε, δ) → (0, 0).
Throughout the remainder of this section, we concentrate on elucidating the emergence and robustness of small amplitude oscillations (SAOs) in the vicinity of the upper CDH (see Figure 18).During phase 3 ○, the singular orbit (green) traces M a SS and jumps down at the fold point.The full trajectory Γ (ε,δ) does not immediately jump at the fold.Instead, there is a delay before Γ (ε,δ) jumps to the lower attracting manifold M L S .Small amplitude oscillations are observed as the trajectory passes through the neighborhood of the canard and DHB points.For smaller ε or δ perturbations (e.g., top row in Figure 23 and Figure 24), the delay is more substantial, but the small oscillations are very small and difficult to observe due to the stronger attraction to M a SS during phase 3 ○.We omit explaining the solution dynamics of (2) for g syn = 4.4 during other phases as they are similar to those observed when g syn = 4.3.In particular, the absence of SAOs near the lower CDH point is due to the same mechanism as discussed in subsection III C.
At default parameters, δ = O(ε).The emergence of SAOs for g syn = 4.4 is governed by both the canard dynamics due to folded node singularities and the slow passage effects associated with the DHB on the upper fold L s , as seen in examples considered in 29,43 .To understand why folded node singularities play an important role in the occurrence of SAOs, we view the trajectory and folded singularities projected onto (V 1 , V 2 , w 2 )-space, as illustrated in Figure 20, as well as draw the funnel volume associated with the folded node singularity curve on the upper fold (see Figure 21).On the other hand, to grasp the importance of the DHB towards the generation of SAOs, we examine the trajectory on the (V 1 , V 2 ) projection, which includes the periodic orbit branches born at the upper Hopf bifurcation (see Figure 22).For simplicity, we omit depicting singular orbits in the following figures and focus only on M S , M SS and their relations to the full trajectory.

B. MMOs are organized by both canard and DHB mechanisms 1. Canard Dynamics
Figure 20 shows the solution trajectory of (2) projected onto (V 1 , V 2 , w 2 )-space when g syn = 4.4.The upper folded singularity curve M comprises folded singularities of two types: folded nodes (solid green) and folded saddles (dashed green).The blue star denotes the saddlenode bifurcation FSN 1 , which occurs O(δ) away from the upper CDH at the blue diamond.Different from the previous case (g syn = 4.3) where folded singularities did not contribute to the SAO dynamics, the current solution trajectory crosses the upper M at a folded node, which can play a critical role in organizing the MMOs.
To confirm that the MMOs exhibit the characteristics of canard dynamics due to the folded node, we plot the funnel volume corresponding to the folded node singularity curve on the upper fold L s (see Figure 21).As discussed in §II B, the funnel for a folded node point represents a two-dimensional trapping region on M S .The 2D funnels for all folded nodes on the upper fold together form a three-dimensional funnel volume.In the (V 1 , V 2 , w 2 ) projection (see the top panel of Figure 21), the funnel volume is bounded by the singular strong ca- nard surface (shown in magenta) and the upper fold surface (in blue).The bottom panel shows the projection onto the (V 1 , V 2 )-plane, where the funnel is indicated by the shaded region.Trajectories initiating inside the funnel (e.g., the cyan and black curves) are filtered through the CDH region and there exist SAOs, whereas trajectories starting outside the funnel (e.g., the yellow curve) cross the fold L s at a regular jump point and there are no SAOs.These observations suggest that the MMOs for g syn = 4.4 exhibit canard-like features and are organized by the canard mechanism.Next, we elucidate that the DHB mechanism also contributes to the occurrence of SAOs in the MMO solution.

Delayed Hopf bifurcation mechanism
Figure 22 shows the projection of the solution trajectory and the bifurcation diagram of the slow layer problem (13) onto the (V 2 , V 1 )-plane.Starting at the green square, the trajectory exhibits SAOs as it follows the upper branch of M SS towards the left.As the trajectory passes through the attracting region of M SS , the oscillation amplitude decays in a typical DHB fashion.Moreover, the orbit experiences a delay along the repelling M SS for an amount of time as it passes through the HB.It is worth noting that after the trajectory enters the unstable part of M SS , there are no symmetric oscilla-  18 (black) and the singular 3D funnel volume corresponding to the curve of folded nodes onto (w2, V2, V1)-space (top) and (V2, V1)plane (bottom).The 3D funnel volume is bounded between the singular canard surface (magenta surface) and the fold surface Ls.The black, cyan and yellow curves denote solutions of (2) with different initial conditions.Other color coding and symbols have the same meaning as in Figure 20.
tions with respect to the DHB, i.e., there are no SAOs with increasing amplitudes.This is due to the fact that V 2 switches from a superslow to a slow timescale at the green triangle, as discussed earlier when g syn = 4.3.
Thus, for δ = O(ε), the MMOs for g syn = 4.4 exhibit characteristics of both the canard and DHB mechanisms.To further confirm this, we performed two perturbations on the system.Firstly, we increased ε by raising the value of C 1 to make the canard mechanism less relevant.This is because increasing C 1 slows down the evolution speed of V 1 , which in turn drives the three-timescale system (2) closer to (3S, 1SS) splitting.We observed that MMOs persisted for C 1 as large as 80 (ε = O( 1)), at which the folded singularities were no longer relevant (see Figure 23).Secondly, we increased δ by raising φ 2 to drive the system closer to (1F, 3S) splitting, and observed that SAOs persisted for φ 2 as large as 0.01 (δ = O( 1)), at which the DHB mechanism was no longer relevant.The persistence of SAOs even when one of the two mechanisms vanishes further highlights the coexistence and interplay of the canard and DHB mechanisms for supporting SAOs.

C. Effects of varying ε and δ on MMOs
Unlike the sensitivity of MMOs at g syn = 4.3 to timescale variations, the interaction of canard and DHB mechanisms due to the existence of the upper CDH when g syn = 4.4 makes MMOs much more robust, as discussed above and illustrated in Figure 3B.Specifically, MMOs persist over biologically relevant ranges of C 1 ∈ (0.1, 80) and φ 2 ∈ (1e − 4, 0.01) (also see Figures 23 and 24).
The robustness of MMOs or SAOs to decreasing ε or δ is expected, as it moves either the DHB point or the folded saddle-node singularity closer to the CDH, causing them to move into the midst of the small oscillations (see Figure 23A and Figure 24A).In this subsection, we only explain the effects of increasing ε or δ on the features of small oscillations within MMOs, as decreasing them yields analogous effects but reversed.We find that MMOs with δ = O(ε) exhibit both canard-and DHBlike features, whereas by tuning δ ≥ O( √ ε), DHB-like features diminish and the canard mechanism dominates.
Below we summarize the effects of increasing the singular perturbations: • For fixed φ 2 = 0.001, increasing ε (i.e., increasing C 1 ) enhances the DHB-like features of the SAOs (see Figure 23).

Increasing C1 makes DHB dominate
Increasing C 1 drives the three-timescale system (2) closer to (3S, 1SS) splitting.As a result, the critical manifold M s and folded singularities including the CDH point become less meaningful and eventually irrelevant for sufficiently large C 1 (e.g., C 1 = 80).Despite this, MMOs persist due to the existence of the DHB mechanism.Moreover, the DHB mechanism becomes more dominant in controlling the features of the small oscillations as C 1 increases, while the influence of the canard mechanism becomes less significant.
Figure 23 shows the effect of increasing C 1 on voltage traces of the full system and the bifurcation diagrams of the fast subsystem (13).As C 1 increases from panel (A2) to panel (C2), the upper DHB (red circle) moves away from the upper CDH (blue diamond) to smaller V 2 values, whereas the CDH points and slow/superslow timescale transitions denoted by yellow and green symbols all remain unaffected by C 1 .As a result, the trajectory with larger C 1 begins small oscillations at a larger distance from the Hopf bifurcation point, similar to what we observed in the case of g syn = 4.3.Furthermore, we have noticed that trajectories for larger C 1 exhibit more pronounced DHB-like characteristics, including SAOs with decreasing amplitudes and a more extended travel distance along the unstable branch of M SS .As discussed before, due to a switch of the V 2 timescale at the green triangle, there are no oscillations with growing amplitudes as one would expect in a typical DHB fashion.

Increasing φ2 makes the canard mechanism dominate
Increasing φ 2 speeds up the superslow variable w 2 and hence drives the three-timescale system (2) closer to (1F, 3S) splitting.As a result, the superslow manifold M SS and the DHB points become less relevant and eventually no longer meaningful for sufficiently large φ 2 (e.g., φ 2 = 0.01).Nonetheless, MMOs continue to persist due to the existence of the canard mechanism.Figure 24 reveals the presence of canard-like features for φ 2 values across a wide range (from 0.0005 to 0.01).That is, trajectories within the funnel volume (cyan and black curves) exhibit SAOs near the CDH, whereas those outside the funnel (yellow curves) display no oscillations.
Moreover, as φ 2 increases, small oscillations tend to pull away from M SS and lose their DHB-like features (i.e., oscillations with decaying amplitude and a delay after passing the HB).Specifically, MMOs with δ = O(ε) (e.g., Figures 20 and 24A) exhibit both canard and DHB characteristics.When δ ≈ 1 2 √ ε (Figure 24B), there are still some DHB-like features.Further increasing δ to δ = 0.006 ≈ √ ε or δ = 0.01 ≈ ε 1 4 (Figure 24C and D), the trajectories no longer closely follow M SS and the amplitudes of small oscillations become almost constant, which reflects the absence of DHB-like features.
In summary, increasing φ 2 makes the canard mechanism dominate and the SAOs exhibit fewer DHB-like features.Conversely, decreasing φ 2 brings the solution and M SS closer together and amplifies the DHB characteristics of the sustained MMOs (see Figure 24).

V. DISCUSSION
Mixed mode oscillations (MMOs) are commonly exhibited in dynamical systems that involve multiple timescales.These complex oscillatory dynamics have been observed in various areas of applications and are well studied in two-timescale settings 2,5,11,16,32,33,38,40,49 .In contrast, progress on MMOs in three-timescale problems has been made only in the recent past (see e.g., [18][19][20][21]24,26,29,43 ). In this wor, we have contributed to the investigation of MMOs in three-timescale settings by considering a four-dimensional model system of coupled Morris-Lecar neurons that exhibit three distinct timescales.We have investigated two types of MMO solutions obtained with different synaptic strengths (g syn = 4.3 mS/cm 2 and g syn = 4.4 mS/cm 2 ).Applying geometric singular perturbation theory and bifurcation analysis 13,31 , we have revealed that the two MMOs exhibit different mechanisms, leading to remarkably different sensitivities to variations in timescales (see Figure 3).In particular, the MMO solution at g syn = 4.3 only involves one mechanism (the delayed Hopf (DHB)) and appears to be vulnerable to certain timescale perturbations, whereas for g syn = 4.4, two separate mechanisms (canard and DHB) coexist and interact to produce more robust MMOs.
The existence of three distinct timescales leads to two important subjects: the critical (or slow) manifold M S and the superslow manifold M SS .The point where a fold of the critical manifold L s intersects the superslow manifold M SS is referred to as the canard-delayed-Hopf (CDH) singularity, which naturally arises in problems that involve three different timescales.Ref. 29 considered a common scenario in three-timescale systems where the CDH singularity exists and proved the existence of canard solutions near the CDH singularity for sufficiently small ε and δ.Moreover, small-amplitude oscillations (SAOs) constantly occur in the vicinity of a CDH singularity 29,43 .
In this work, we have reported several key findings that have not been previously observed in three-timescale systems.First, although we have identified the same type of ).Increasing φ2 preserves MMOs with more canard-like feature.Codings and symbols on the right panels are the same as in Figure 21.
as the further reduction of ε resulted in the cross of the DHB with the green square or a complete vanish of the DHB, we observed a recovery of MMOs originating from the saddle foci equilibria along the superslow manifold M SS (see Figures 16B,C and 17).Eventually, SAOs disappeared entirely when V became so rapid that it failed to remain in proximity of M SS to generate small oscillations.
As one would expect, MMOs when g syn = 4.3 is also sensitive to increasing δ, which speeds up the superslow variable and thus makes the DHB mechanism less relevant.Interestingly, however, increasing δ does not just simply eliminate MMOs.Instead, it led to a total of five transitions between MMOs and non-MMO states (see Figures 3 and 14B).Our analysis suggests that these complex transitions occur due to a spike-adding like mechanism.When no big spike is lost with the increase of δ, a transition from MMOs to non-MMOs will take place.
However, if an entire big spike is lost during this process, changes to the SAOs will be reversed and MMOs will recover again.Ultimately, MMOs will be completely lost as δ is increased to a point where the DHB mechanism is no longer relevant.
On the other hand, MMOs at g syn = 4.3 show strong robustness to increasing ε or decreasing δ.This is not surprising as both of these changes manifest the DHB mechanism by moving the three timescale problem closer to (3S, 1SS) separation.As a result, we observed more DHB characteristic in the MMOs as demonstrated in the case of increasing ε (see Figures 13A and 15).Nonetheless, decreasing δ led to an interesting non-monotonic effect on SAOs, where the amplitude and the number of SAOs exhibit an alternation between increase and decrease as δ is reduced (see Figure 14A).Our analysis showed that the mechanism underlying such non-monotonic effects on SAOs over the decrease of δ is similar to the mechanism that drives multiple MMOs/non-MMOs transitions as δ is increased.
Unlike g syn = 4.3, an upper CDH occurs near the SAOs at g syn = 4.4.In this case, we showed that this CDH allowed for the coexistence and interaction of canard and DHB mechanisms, resulting in MMOs with strong robustness against timescale variations.Since both MMO solutions for g syn = 4.3 and 4.4 exhibit DHB mechanisms, they show similar responses and robustness to increasing ε and decreasing δ, both of which lead to more DHB characteristic in the MMOs.In contrast to g syn = 4.3, MMOs at g syn = 4.4 are also robust against changes that speed up fast or superslow variables.This is due to the existence of the canard mechanism in addition to the DHB.Instead of eliminating the upper DHB in the case of g syn = 4.3 (Figure 11A), decreasing ε at g syn = 4.4 brings the DHB point closer to the CDH (Figure 19, right panel) and, consequently, to the midst of small oscillations (Figure 23).Moreover, this timescale change manifests the canard mechanism by moving the system closer to the ε → 0 singular limit.Hence, MMOs at g syn = 4.4 persist as ε decreases and are organized by both mechanisms.On the other hand, increasing δ diminishes the relevance of the DHB mechanism and moves the trajectory further away from the superslow manifold, resulting in MMOs with more canard-like features.
To summarize, the coexistence of canard and DHB mechanisms for MMOs at g syn = 4.4, due to the presence of a nearby upper CDH, leads to significantly enhanced robustness against timescale variations compared with g syn = 4.3, where only one mechanism is present.While we have not examined the case when only a canard mechanism is present, based on our analysis, we expect that MMOs with only canard mechanism would show more sensitivities to timescale variations that reduces the relevance of the canard mechanism such as increasing ε or decreasing δ.Similarly, such MMOs should show stronger robustness to decreasing ε or increasing δ.It would be of interest to explore such a scenario for future investigation.Furthermore, we did not notice any non-monotonic effects on the features of SAOs as φ 2 is decreased when g syn = 4.4, unlike what we observed in g syn = 4.3.This difference is likely attributed to the presence of an additional canard mechanism which may have hindered the occurrence of complex non-monotonic behaviors.A complete analysis of this phenomenon could be investigated in future work.where From system (A1), we can see that the voltage V 1 evolves on a fast timescale ( C1 gmax = 1 ms), (w 1 , V 2 ) evolve on a slow timescale ( min(τw(V1)) φ1 ≈ C2 gmax = O(10) ms), whereas w 2 evolves on a superslow timescale ( min(τw(V2)) φ2 = O(100) ms).
Appendix B: Folded saddle-node (FSN) singularities of (18)   To derive the conditions for FSN singularities, we note that the eigenvalues λ of the folded singularities satisfy the following algebraic equation ) where tr(J D ) and det(J D ) denote the trace and determinant of the Jacobian matrix J D of the desingularized system (18), respectively.As discussed before, det(J D ) ≡ 0 along the folded singularity curve M.This also directly follows from the following calculations det where the entries in J D denote partial derivatives.The last equality in the above equations holds because ∂V1∂V2 .Hence the remaining two nontrivial eigenvalues λ w and λ s satisfy It follows that an FSN bifurcation is given by Plugging in functions F , G and H from (18), we can rewrite the above condition as the equation ( 21) in subsection II B, namely, where Next we prove there are two ways an FSN can occur.In case 1, suppose Q = 0.An FSN can occur when f 2 = δK 1 (V 1 , V 2 , w 2 ), where K 1 (V 1 , V 2 , w 2 ) := P/Q.Note that FSN is a folded singularity point (20), that is, F = g 1 − F 1V2 f 2 = 0.This implies g 1 = F 1V2 f 2 = δK 2 (V 1 , V 2 , w 2 ), where K 2 (V 1 , V 2 , w 2 ) := F 1V2 K 1 (V 1 , V 2 , w 2 ).Thus, the first way that an FSN can occur is described as the following.
In case 2, suppose f 2 = 0.An FSN can occur when Q(V 1 , V 2 , w 2 ) = δK 3 (V 1 , V 2 , w 2 ) where K 3 (V 1 , V 2 , w 2 ) = P/f 2 .As a result, the second way that an FSN can occur is defined as which is the equation (23) in subsection II B. This implies that an FSN 2 is O(δ) close to the intersection of M and Q(V 1 , V 2 , w 2 ) = 0.
Appendix C: CDH at double singular limit exhibits two linearly independent critical eigenvectors Recall the Jacobian matrix J D at an FSN 1 , which becomes a CDH at the double singular limit, has two zero eigenvalues.We prove that the center subspace associated with the two zero eigenvalues is two dimensional and is transverse to the fold surface of the critical manifold (see subsection II D and Figure 4).
The Jacobian matrix J D of the desingularized system at the double singular limit (ε → 0, δ → 0) is given by in which G V1 = −F 1V1V1 f 2 and G V2 = −F 1V1V2 f 2 .At a CDH, f 2 = g 1 = 0.It follows that G V1 = G V2 = 0. Thus, the Jacobian matrix J D at a CDH singularity becomes which has a nullity of 2, implying the existence of two linearly independent critical eigenvectors associated with zero eigenvalues.Moreover, the center subspace at the CDH, given by the plane F V1 V 1 + F V2 V 2 + F w2 w 2 = 0 (see pink planes in Figure 4C and D), is transverse to the fold surface L s of the critical manifold.This is because the normal vector of L s at the CDH, which is given by n f = ns |ns| where n s = (F 1V1V1 , F 1V1V2 , 0) is not parallel to the normal vector of the center subspace given by (F V1 , F V2 , F w2 ).Thus, the CDH singularity (an FSN 1 at the double singular limit) considered in system (2) is a saddle-node bifurcation of folded singularities, with the center manifold of the FSN 1 transverse to the fold of the critical manifold.
Finally, we verify a previous claim that we made in subsection II D that the eigenvector associated with the first trivial zero eigenvalue of J D (denoted by v 0 ) is always tangent to L s .Through direct computations, we obtain v 0 = u0 |u0| with The dot product of v 0 and the unit normal vector of L s at the CDH is given by It follows directly that v 0 is tangent to L s as expected.
Appendix D: Decreasing φ2 when gsyn = 4.3 causes non-monotonic effects on SAOs In this subsection, we explain the mechanism underlying the non-monotonic effects of decreasing φ 2 on SAOs (see Figure 14A in subsection III D).We claim that (1) if an additional (full) big spike is generated during the process of decreasing φ 2 , the number of SAOs will decrease and their amplitudes will increase (e.g., when φ 2 decreases from 0.0008 to 0.0007); (2) if no additional big spike is generated during the process, the changes to the SAOs will be reversed (e.g., there are more SAOs with smaller amplitudes when φ 2 decreases from 0.0007 to 0.0006).Below, we examine the two opposite effects using the projections of the corresponding solutions and bifurcation diagrams onto (V 2 , V 1 ) space (see Figure 25).
Figure 25A shows that the amplitude of SAOs with φ 2 = 0.0008 (black) is smaller than SAOs with φ 2 = 0.0007 (magenta).There are also more black SAOs when φ 2 = 0.0008 (see Figure 14A), which however is not obvious in the (V 2 , w 2 , V 1 ) projection as the black SAOs are hardly visible.Compared with the black solution, the magenta one has a slower evolving rate of V 2 , which is slaved to w 2 , during phase 1 ○. Thus, one extra full spike occurs for the magenta trajectory compared to the black trajectory during the jump at phase 2 ○.As a result, the black trajectory, after its last big spike, approaches the black square at its maximum V 2 along M SS .In contrast, the magenta trajectory, due to the extra full spike that occurs during the jump, approaches the maximum V 2 from the outside of the periodic orbit branch.Consequently, the magenta trajectory that is being pushed further away from M SS exhibits fewer small oscillations with larger amplitudes compared with the black trajectory, which remains close to M SS .
When φ 2 is reduced from 0.0007 to 0.0006, the solution trajectory changes to the green curve in Figure 25B.Slowing down of w 2 even further keeps the number of full spikes between the magenta and green trajectories the same, but now all five full green spikes occur within phase 1 ○, similar to the case when φ 2 = 0.0008.As a result, the green and black SAOs exhibit similar characteristics as shown in Figure 25B, resulting in more SAOs with smaller amplitudes than the magenta solution in Figure 25A.In this subsection, we explain the mechanism underlying the MMOs/non-MMOs transitions induced by an increase of φ 2 in the case of g syn = 4.3 (see Figure 14 in subsection III D), which is similar to the mechanism that causes the non-monotonic effects on SAOs when φ 2 is decreased (see Appendix D).
During the increase of φ 2 which speeds up w 2 , the big spikes produced during phase 1 ○ will gradually decrease.If one big (full) spike before phase 3 ○ is lost during this process, there will be more SAOs with smaller amplitudes (e.g., when φ 2 increases from 0.001 to 0.0012 in Figure 26A) or the earlier lost MMOs will reappear (e.g., when φ 2 increases from 0.0016 to 0.0018 in Figure 26C).As φ 2 increases from 0.001 to 0.0012 in Figure 26A, the number of full spikes during phases 1 ○ and 2 ○ decreases by one.Consequently, all three green full spikes occur within phase 1 ○, while there is a large black spike during phase 2 ○.As explained in the previous subsection, this causes the green trajectory to approach the maximum V 2 along M SS and remain close to it after reaching the green square.As a result, more SAOs with smaller amplitudes are observed in the green trajectory compared to the black trajectory.Similarly, when the increase of φ 2 from 0.0016 to 0.0018 leads to a decrease in the number of big spikes that occur before phase 3 ○, more SAOs are observed and the transition from non-MMOs to MMOs occurs (see Figure 26C).
In contrast, increasing φ 2 from 0.0012 to 0.0016 does not change the number of large spikes before phase 3 ○ (compare the green and magenta trajectories in Figure 26A and B).However, speeding up of φ 2 pushes the third big magenta spike to occur during the jump, similar to the black trajectory.This causes the trajectory to be pushed away from M SS , resulting in fewer SAOs with larger amplitudes.In fact, in this case, SAOs in the magenta trajectory are eliminated and therefore a transition from MMOs to non-MMOs results.

FIG. 2 .
FIG. 2. Bifurcation curves of DHB (red) and CDH singularities (blue) for (2) with respect to gsyn.Specifically, these curves represent the upper DHB and upper CDH, corresponding to larger V1 and V2 values.The lower CDH and DHB, associated with smaller Vi values, are not presented here.The two vertical asymptotes are given by gsyn ≈ 4.2628 and gsyn ≈ 4.3213.

FIG. 3 .
FIG.3.Regions of MMOs (yellow) and non-MMOs (blue) of the full system (2) in (φ2, C1)-space for (A) gsyn = 4.3 and (B) gsyn = 4.4.Increasing C1 slows down the fast variable V1, whereas increasing φ2 speeds up the superslow variable w2.The timescales of w1 and V2 remain unaffected.The red star marks the default parameter values of C1 and φ2 as given in TableI.(A) gsyn = 4.3.While MMOs are robust to increasing C1 and decreasing φ2, decreasing C1 or increasing φ2 leads to multiple transitions between MMOs and non-MMOs (crossings between the dashed lines with the yellow/blue boundaries).(B) gsyn = 4.4.MMOs are robust to changes of both C1 and φ2 over the ranges of 0.1 ≤ C1 ≤ 20 and 0.1 ≤ φ2 ≤ 8e − 3. Note that the MMOs at gsyn = 4.4 will eventually vanish for C1 and φ2 large enough at which there is no more timescale separation (data not shown).

FIG. 4 .
FIG.4.Projections to (V1, V2, w2)-space of the critical manifold fold surfaces Ls (blue surface) for (A, C) gsyn = 4.3 and (B, D) gsyn = 4.4.Also shown are the curves of folded singularities M including folded node (solid green), folded saddle (dashed green), and two types of folded saddle-nodes FSN (blue star: FSN 1 ; cyan star: FSN 2 ).The yellow curve consists of mostly folded foci points and small segments of other singularities (e.g., folded node, folded saddle) that are barely visible and hence are not displayed here.In the top two panels (A, B) when δ = 0, an FSN 1 point (blue star) is O(δ) close to a CDH singularity (blue diamond), whereas an FSN 2 (cyan star) is far away from any CDH.In the lower two panels (C, D) at the singular limit δ = 0, the FSN 1 point becomes a CDH singularity (blue star overlapping with blue diamond).The center subspace of an FSN 1 (resp., an FSN 2 ) is denoted by a pink plane (resp., a yellow plane).It follows that the center manifolds of both FSN 1 and FSN 2 are transverse to Ls.

FIG. 7 .
FIG. 7. Projection of an attracting MMO solution trajectory (black curve) of system (2) for gsyn = 4.3 from Figure 1 to (V1, w1, V2)-space.Also shown are parts of the singular orbit from Figure 6A (green curves), the critical manifold MS (blue surface), folds of MS (blue curves) and the superslow manifold MSS (red curves).The upper inset shows the upper branch of MSS is always below the upper fold of MS, indicating the absence of an upper CDH.The black circle near the upper fold is the true equilibrium of the full system (2), which is unstable.The lower inset shows a magnified view around the lower CDH (blue diamond) at which the lower branch of MSS intersects the lower fold of MS.The lower HB bifurcation (red circle) is O(ε) close to the CDH singularity.Other color coding and symbols have the same meaning as in Figures 5 and 6.

FIG. 8 .
FIG. 8.The solution from Figure7projected onto (V1, V2, w2)-space.Also shown are the projections of the fold surface Ls (blue surface), superslow manifold MSS (red curve) and folded singularities M (green and yellow curves).The cyan star denotes FSN 2 defined in(23).Color codings of the folded singularity curves are the same as in Figure4.Other symbols have the same meanings as in Figure7.

FIG. 9 .FIG. 10 .
FIG. 9.The solutions from Figure 7, and the bifurcation diagram of the slow layer problem (13) projected onto (V2, V1)-plane.(A) Besides the trajectory (black), the singular orbit (green) and MSS (red curve), also shown is the periodic orbit (PO) branches (solid cyan: stable; dashed cyan: unstable) born at the upper HB bifurcation.The upper two yellow circles denoting the saddle-node bifurcations of MSS exhibit the same V2 values as the folds of the V2-nullcline (green circle and green triangle).The lower yellow circle (see the inset) represents the actual fold of MSS (denoted as L 2 ss ) and is not a fold of the V2-nullcline.Other symbols have the same meanings as in Figure 7. (B) Real part of the eigenvalues of the upper triangular block of JSL (24), the Jacobian matrix of the slow layer problem (13) along the superslow manifold MSS from panel (A).The eigenvalues along MSS are real when there are two branches of Re(λ) (blue curves) and complex when there is a single branch of Re(λ) (black curve).
ure 11B and recall Remark II.3).When ε increases, the lower Hopf and the lower fold of M SS meet and coalesce at a Bogdanov-Takens (BT) bifurcation.After the BT bifurcation, the Hopf bifurcation disappears.Unlike the upper DHB, the lower DHB is close to the actual fold of M SS (also see Figure11A, lower yellow circle).

FIG. 18 .FIG. 19 .
FIG. 18. Projection of a solution trajectory (black curve) of system (2) for gsyn = 4.4 from Figure 1B to (V1, w1, V2)-space.Also shown are parts of the singular orbit Γ F (0,0) ∪ Γ S (0,0) ∪ Γ SS (0,0) from Figure 6D (green curve), the critical manifold MS (blue surface), folds Ls of MS (blue curves) and the superslow manifold MSS (red curves).The solid (resp.dashed) red curves represent the attracting (resp.repelling) branches of MSS. Green symbols mark the transitions between slow and superslow pieces of the (V2, w2) oscillations as in Figure 5.The blue diamonds denote CDH singularities -the intersection of Ls and MSS.The upper CDH is a folded node and the lower CDH (almost overlapping with the lower DHB denoted by a red circle) is a folded focus.The upper CDH point is O(δ) close to the folded saddle-node singularity FSN 1 (blue star) and O(ε) close to the upper DHB (red circle).The black circle denotes the isolated ordinary singularity of the full system, whose type is a saddle-focus.

FIG. 20 . 30 FIG. 21 .
FIG. 20.Projection of the solution from Figure 18 onto (w2, V2, V1)-space.Also shown are the projections of the fold surface Ls (blue surfaces), superslow manifold MSS and folded singularities M (curves of green and yellow).Other color coding and symbols have the same meaning as in Figures 4 and 18.

FIG. 22 .
FIG. 22. Projection of the solution trajectory and superslow manifold MSS (red curve) from Figure 18 onto (V2, V1)-plane.The solid (resp.dashed) cyan curves denote stable (resp.unstable) periodic orbit branches.Yellow circles represent saddle-node bifurcations of MSS, in which the upper two have the same V2 values as the folds of the V2-nullcline at green circle and triangle.Other colors and symbols are the same as in Figure 18.

TABLE I .
The values of the parameters in the model given by (2) and (3).

TABLE II .
Effects of varying ε and δ on MMOs.