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

Authors: Ngoc Anh Phan, Yangyang Wang

Abstract

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.

One of the most common types of complex oscillatory dynamics observed in systems with multiple timescales is mixed-mode oscillations (MMOs). MMOs are characterized by patterns that involve the interspersion of small-amplitude and large-amplitude oscillations. Over the years, the theory of MMOs in fast–slow systems has been well developed. Recently, there has been more progress on the analysis of MMOs in three-timescale systems. Nonetheless, MMOs in the latter case are still much less understood. In this work, we contribute to the investigation of MMOs in the three-timescale settings by considering coupled Morris–Lecar neurons. We uncover the properties and geometric mechanisms underlying two different MMO patterns in our three-timescale system. One of them involves the interaction of the two distinct MMO mechanisms, showing a high degree of robustness to timescale perturbations, whereas the other lacks such mechanism and is thus vulnerable to timescale variations. Based on our analysis, we establish conditions that lead to more robust generation of MMOs in three-timescale problems, particularly against perturbations in timescales.

INTRODUCTION

Mixed-mode oscillations (MMOs) are frequently perceived in the dynamical systems involving multiple timescales;^1^ these are complex oscillatory dynamics characterized by the concatenation of small-amplitude oscillations (SAOs) and large-amplitude excursions in each periodic cycle. Such phenomena have been recognized in many branches of sciences including physics, chemistry,^2,3^ and particularly life sciences.^4–17^

Theoretical analysis of MMOs in systems with two distinct timescales has been well developed with the implementation of the geometric singular perturbation theory (GSPT);^18^ see Ref. 1 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^19,20^ and a slow passage through the delayed Andronov–Hopf bifurcation (DHB) of the fast subsystem.^21–24^ While in two-timescale settings, these two mechanisms remain separated, they can coexist and interact in three-timescale regime.^14,25–27^

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.^28^ However, many real-world systems have more than two timescales.^25,27,29–34^ It has also been established that a two-timescale decomposition fails to capture certain aspects of the system’s dynamics.^35^ Therefore, classifying three timescales into two groups is not a sufficient approach for modeling and analysis.

MMOs in three-timescale systems have been studied before (see, e.g., Refs. 10, 30, 11, 27, 25, 26, and 36–42 and references therein). Initial approaches were to consider three-dimensional systems,

εdxdt=f(x,y,z),dydt=g(x,y,z),dzdt=δh(x,y,z), (1)

with special cases δ=ε or δ=ε.^10,25,30,36,41^ MMOs were shown to emerge through an effect analogous to a slow passage through a canard explosion.^25,30,41^ More recently, there has been a growing interest in MMOs with independent singular perturbation parameters ε and δ, as explored in various three-dimensional models.^26,38–40^ In particular, Ref. 26 centered on a novel singularity type denoted as canard-delayed-Hopf (CDH) singularity, which naturally arises in three-timescale settings when the two mechanisms for MMOs (the fast subsystem Hopf and a folded node) coexist and interact. The authors investigated the existence and properties of MMOs near the CDH singularity.

In this paper, we contribute to the investigation of MMOs in three-timescale settings by considering a model of four-dimensional coupled Morris–Lecar neurons^43,44^ that was introduced by Ref. 35. The model equations are given by

dV1dt=(I1−gCam∞(V1)(V1−VCa)−gKw1(V1−VK)−gL(V1−VL)−gsynS(V2)(V1−Vsyn))/C1,dw1dt=ϕ1(w∞(V1)−w1)/τw(V1),dV2dt=(I2−gCam∞(V2)(V2−VCa)−gKw2(V2−VK)−gL(V2−VL))/C2,dw2dt=ϕ2(w∞(V2)−w2)/τw(V2), (2)

with

S(Vi)=α(Vi)/(α(Vi)+β),α(Vi)=1/(1+exp(−(Vi−θs)/σs)),m∞(Vi)=0.5(1+tanh((Vi−K1)/K2)),w∞(Vi)=0.5(1+tanh((Vi−K3)/K4)),τw(Vi)=1/cosh((Vi−K3)/2K4). (3)

Table I lists the parameter values for the model chosen to ensure that (2) exhibits three distinct timescales, where V1 is fast, w1,V2 are slow, and w2 is superslow. In a more biologically realistic model for calcium and voltage interactions, V1 might represent membrane potential, while V2 might represent intracellular calcium concentration with appropriate adjustments to parameter units and functional terms (see, e.g., Refs. 32 and 33). For the physiological description of functions in (2) and (3), we refer readers to Ref. 35 for details.

In the absence of coupling gsyn=0, (V1,w1) is excitable with an attracting critical point at relatively low V1 value, whereas (V2,w2) is oscillatory with an attracting limit cycle independent of the value of gsyn and the dynamics of (V1,w1). To analyze the three-timescale coupled Morris–Lecar neurons, Ref. 35 extended two approaches previously developed in the context of GSPT for the analysis of two-timescale systems to the three-timescale setting and showed these two approaches complemented each other nicely. By varying gsyn 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. 35 (e.g., gsyn=4.1 and 5.1 in Fig. 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 gsyn=4.3 and gsyn=4.4 (with the unit of mS/cm2), as highlighted in Fig. 1. To provide further insight into our choice of gsyn values, we perform a bifurcation analysis to explore the effect of gsyn on a singularity called canard-delayed-Hopf (CDH) that was first introduced by Ref. 26 (see Fig. 2, blue). As noted above, this singularity plays a crucial role in organizing MMOs within the three-timescale setting. We show in Sec. II that the full system (2) may exhibit two CDH points—one at larger Vi values, i∈{1,2} (denoted as upper CDH) and the other at smaller Vi (denoted as lower CDH, see Fig. 4). Similarly, (2) may exhibit an upper DHB and a lower DHB. However, we demonstrate in Secs. III and IV that only the upper CDH or DHB can support MMOs.

In Fig. 2, we plot the bifurcation curves of the upper CDH and the upper DHB in the (gsyn,V2) plane, both of which approach vertical asymptotes as gsyn decreases. Our first choice of gsyn=4.3 represents gsyn values between the two asymptotes, at which the upper DHB exists but there is no upper CDH. In contrast, gsyn=4.4 lies to the right of the CDH asymptote and serves as a representative scenario where both the upper CDH and DHB exist and facilitate the existence of MMOs. When gsyn<4.2628 (e.g., gsyn=1 and 4.1 as considered in Ref. 35), the system (2) does not produce MMOs. When gsyn>4.3213, MMOs may arise through mechanisms similar to those we will thoroughly examine for gsyn=4.4 in Sec. IV. As we shall see later, the emergence of these oscillations is contingent upon the system’s trajectory closely approaching the DHB and CDH points. In situations where the trajectory remains distant from these critical points, such as gsyn=5.1, MMOs are not observed (see Fig. 1). Additionally, the absence of damped oscillations is also influenced by the real eigenvalues within the subsystem that governs the transition from V1 spikes to a V1 plateau, as detailed in Ref. 35. Therefore, as gsyn increases, the system exhibits multiple transitions between MMOs and the absence of MMOs, affected by the varying proximity of the trajectory to the CDH and DHB near the spiking/plateauing transition point. Eventually, MMOs vanish entirely at sufficiently large gsyn values (e.g., gsyn=80), where both CDH and DHB points fall to low V2 values that lie beyond the system’s trajectory.

In this paper, we focus only on MMOs at gsyn=4.3 and gsyn=4.4. To fully understand these MMO dynamics, we employ the extended GSPT,^18,35^ a geometric approach to multiple-timescale systems that enables the prediction of the full system dynamics based on lower-dimensional subproblems. As the first step of the GSPT approach, we perform a dimensional analysis of (2) to reveal the important timescales. This transforms (2) to the following three-timescale

εdV1dts=f1(V1,w1,V2),dw1dts=g1(V1,w1),dV2dts=f2(V2,w2),dw2dts=δg2(V2,w2), (4)

where ε=0.1, δ=0.053, ts is the slow dimensionless time variable, f1, f2, g1, and g2 are O(1) functions specified in (A1) in Appendix A, which include details of the nondimensionalization procedure. For simplicity, we did not rescale V1 and V2 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 V1 evolves on a timescale of O(ε−1), (w1,V2) on O(1) and w2 on O(δ). Introducing a superslow time tss=δts yields an equivalent description of dynamics,

εδdV1dtss=f1(V1,V2,w1),δdw1dtss=g1(V1,w1),δdV2dtss=f2(V2,w2),dw2dtss=g2(V2,w2), (5)

which evolves on the superslow timescale and is called the superslow system. Alternatively, defining a fast time tf=ts/ε, we obtain the following fast system:

dV1dtf=f1(V1,V2,w1),dw1dtf=εg1(V1,w1),dV2dtf=εf2(V2,w2),dw2dtf=εδg2(V2,w2), (6)

which evolves on the fast timescale.

The presence of two independent singular perturbation parameters, ε and δ, indicates there are various ways to implement GSPT, thereby leading to multiple singular limit predictions as we describe in Sec. II. Our analysis suggests that the two MMO solutions at gsyn=4.3 and gsyn=4.4 arise from distinct mechanisms, resulting in remarkably different sensitivities to variations in timescales, as illustrated in Fig. 3.

When gsyn=4.3, we show that there is no interaction of different MMO mechanisms due to the lack of a nearby CDH singularity (see Fig. 2). We demonstrate that only the δ→0 singular limit provides a faithful prediction for the observed MMOs, whereas the ε→0 limit does not. This observation pinpoints the DHB from the δ→0 limit as the only mechanism for the MMOs at gsyn=4.3. As a result, these MMO dynamics are sensitive to variations in timescales [see Fig. 3(a)]. Specifically, we show that MMOs persist for δ≤O(ε). Increasing δ via increasing ϕ2 or decreasing ε via decreasing C1 to a degree where δ>O(ε) leads to multiple MMOs/non-MMOs transitions. MMOs are completely lost when δ≥O(ε13) for which the DHB is no longer present.

In contrast, there exists a CDH in the middle of the SAOs when gsyn=4.4 as discussed above. We demonstrate that this CDH allows the fast subsystem Hopf and a canard point to coexist and interact to co-modulate properties of the local oscillatory behavior. In this case, both the ε→0 and δ→0 singular limits contribute faithful predictions for the observed dynamics, resulting in MMOs with significantly stronger robustness than gsyn=4.3 [Fig. 3(b)]. We demonstrate that when δ=O(ε), MMOs exhibit both canard and DHB features. Upon tuning δ≥O(ε), DHB-like features disappear and the canard mechanism dominates. 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 produce any local oscillations.

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. 39 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^26,27^ where the CDH always leads to the 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 local oscillations regardless of perturbation sizes ε and δ.

The paper is organized as follows. In Sec. 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 Sec. II. In Sec. III, we investigate the mechanism and sensitivity of MMOs when gsyn=4.3 to variations in perturbation sizes ε or δ [i.e., varying C1 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 Fig. 3(a). While Fig. 3(a) 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 Sec. IV, we uncover the dynamic mechanism underlying MMOs from (2) when gsyn=4.4. In this case, the existence of a CDH in the middle of the SAOs enables two different mechanisms to co-modulate the properties of MMOs. We explain why MMOs organized by a CDH singularity as seen in the case of gsyn=4.4 exhibit remarkable robustness against variations in timescales [see Fig. 3(b)]. Finally, we conclude in Sec. V with a discussion.

GEOMETRIC SINGULAR PERTURBATION ANALYSIS

In this section, we apply the extended geometric singularity perturbation analysis^18,35^ to the three-timesale coupled Morris–Lecar system (4) by treating ε as the only singular perturbation parameter,^19,20^ treating δ as the only singular perturbation parameter,^21–24^ and finally treating ε and δ as two independent singular perturbation parameters.^26,27,35^

Although the detailed GSPT analysis and derivation of subsystems have been previously presented in Ref. 35, 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 Ref. 35. Specifically, we concentrate on reviewing and discussing the canard mechanism in Sec. II B, delayed Hopf bifurcation in Sec. II C, and their interactions in Sec. II D.

Singular limits

ε→0 Singular limit

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, V1, for fixed values of the other variables,

dV1dtf=f1(V1,w1,V2). (7)

The set of equilibrium points of the fast layer problem is called the critical manifold and is denoted as MS,

MS:={(V1,w1,V2,w2):f1(V1,w1,V2)=0}. (8)

Although MS is a three-dimensional (3D) manifold in R4 space, it does not depend on w2. We can solve f1(V1,w1,V2)=0 for w1 as a function of V1 and V2 and can, therefore, represent MS as

w1=F1(V1,V2) (9)

for a function F1. It is well known that, for sufficiently small ε>0, normally hyperbolic parts of MS each perturb to a locally invariant manifold called a slow manifold, on which w1 is given by an O(ε)-perturbation of F1 (Ref. 18); we simply use MS as a convenient numerical approximation of these slow manifolds.

MS is a 3D folded manifold with two-dimensional (2D) fold surface, Ls, given by

Ls:={(V1,w1,V2,w2)∈MS:F1V1=0} (10)

or equivalently

Ls:={(V1,w1,V2,w2)∈MS:f1V1=0}, (11)

where F1V1 and f1V1 denote the partial derivatives of F1 and f1 with respect to V1. The fold surface divides the critical manifold MS into attracting upper (MSU) and lower (MSL) branches where F1V1<0 and repelling middle branch (MSM) where F1V1>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 w1,V2,w2 along MS,

dw1dts=g1(V1,w1),dV2dts=f2(V2,w2),dw2dts=δg2(V2,w2), (12)

where f1=0.

δ → 0 Singular limit

Alternatively, fixing ε>0 and taking δ→0 in the slow system (4) yields the 3D slow layer problem in the form

εdV1dts=f1(V1,w1,V2),dw1dts=g1(V1,w1),dV2dts=f2(V2,w2), (13)

where the superslow variable w2 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 MSS,

MSS:={(V1,w1,V2,w2):f1(V1,w1,V2)=g1(V1,w1)=f2(V2,w2)=0}. (14)

MSS is a 1D subset of MS. Similar to MS, the normally hyperbolic parts of MSS perturb to nearly locally invariant manifolds for δ sufficiently small. Later in Sec. II C, we will discuss the bifurcations of the slow layer problem (13), i.e., nonhyperbolic regions on MSS where Fenichel’s theory (GSPT) breaks down.

Taking the same singular limit in the superslow system (5) leads to the 1D superslow reduced problem

dw2dtss=g2(V2,w2), (15)

where f1=g1=f2=0. The superslow motions of trajectories of (15) are slaved to MSS until nonhyperbolic points are reached.

ε→0,δ→0 Double singular limits

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

0=f1(V1,V2,w1),dw1dts=g1(V1,w1),dV2dts=f2(V2,w2), (16)

which describes the slow motion along MS and the superslow variable w2 is fixed as a constant.

It follows that the double singular limits lead to three 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 Ref. 35.

Slow reduced problem and canard dynamics

To investigate canard dynamics, we project the slow reduced problem (12) onto (V1,V2,w2) to obtain a complete description of the dynamics along MS. To this end, we differentiate the graph representation of MS given by w1=F1(V1,V2) to obtain

F1V1dV1dts=g1(V1,w1)−F1V2f2(V2,w2),dV2dts=f2(V2,w2),dw2dts=δg2(V2,w2), (17)

where F1V2:=∂F1∂V2. Note that the reduced system (17) is singular at the fold surfaces Ls (10). Nonetheless, this singular term can be removed by a time rescaling ts=−F1V1td, and we obtain the following desingularized

dV1dtd=F1V2f2(V2,w2)−g1(V1,w1):=F(V1,w1,V2,w2),dV2dtd=−F1V1f2(V2,w2):=G(V1,V2,w2),dw2dtd=−δF1V1g2(V2,w2):=H(V1,V2,w2,δ). (18)

We observe that the desingularized system (18) is equivalent to (17) on the attracting branch, i.e., for F1V1<0, but has the opposite orientation on the repelling branch, i.e., for F1V1>0.

The desingularized system (18) has two kinds of ordinary and folded singularities. The ordinary singularities are the true equilibria of the full system (2), which is defined by

E:={(V1,w1,V2,w2)∈MS:g1(V1,w1)=f2(V2,w2)=g2(V2,w2)=0}. (19)

For the chosen parameter set in Table I, E always lies on the repelling branch of MS 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 Ls defined by

M:={(V1,w1,V2,w2)∈Ls:F(V1,w1,V2,w2)=0}. (20)

Folded singularities are special points that allow trajectories of (17) to cross the fold Ls with nonzero speed. Such solutions are called singular canards.^19^ Note that when projecting to (V1,w1,V2)-space, the condition F=0 is redundant and the fold surfaces Ls 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 JD, 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 (respectively, with opposite signs) are called folded nodes (respectively, 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 Ls form a two-dimensional trapping region (the funnel) on the attracting branch of MS 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 Ls 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 MS to a repelling MS 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

FSN:={(V1,w1,V2,w2)∈M:f2Q(V1,V2,w2)=δP(V1,V2,w2)}, (21)

where Q and P are defined in (B2) in Appendix B, which contains a detailed derivation of the FSN condition. Similar to Ref. 27, we demonstrate in Appendix B that our system can exhibit an FSN in two different

FSN1:={(V1,w1,V2,w2)∈M:f2=δK1(V1,V2,w2)}={(V1,w1,V2,w2)∈Ls:f2=δK1(V1,V2,w2),g1=δK2(V1,V2,w2)}, (22)

or

FSN2:={(V1,w1,V2,w2)∈M:Q(V1,V2,w2)=δK3(V1,V2,w2)}, (23)

where K1,K2, and K3 are defined in Appendix B. We discuss below in Remark II.1 that both FSN1 and FSN2 are novel types of FSN.^26^

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 JSL denote the Jacobian matrix of (13) evaluated along the superslow manifold MSS, which is given by

JSL=(1εf1V11εf1w11εf1V2g1V1g1w1000f2V2), (24)

where the nonzero entries denote partial derivatives.

The eigenvalues of JSL are given by f2V2 and the eigenvalues of

J=(1εf1V11εf1w1g1V1g1w1).

Thus, the Hopf bifurcation points on MSS are given by tr(J)=1εf1V1+g1w1=0 and detJ>0. The former defining condition can be rewritten as

MSSH:={(V1,w1,V2,w2)∈MSS:f1V1=−εg1w1}. (25)

The isolated fold bifurcation points on MSS are located by letting detJSL=0. That is,

Lss:={(V1,w1,V2,w2)∈MSS:f2V2=0 or f1V1g1w1−f1w1g1V1=0}. (26)

The fold points Lss that satisfy the former condition (denoted as Lss1) are the folds of the V2-nullcline (see Fig. 5, green circle and green triangle), which correspond to the transition between superslow dynamics along MSS and slow jumps. Lss given by the latter condition (denoted as Lss2) corresponds to the actual fold of MSS when projected to (V1,w1,V2)-space. Since g1w1<0, f1w1<0 and g1V1>0, it follows that Lss2 lies on the middle branch of MS ( f1V1=f1w1g1V1/g1w1>0) and hence will not play a role in dynamics. At the double singular limit (ε,δ)→(0,0), the Lss2 fold point of MSS 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.

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

dV1dtd=F(V1,w1,V2,w2),dV2dtd=G(V1,V2,w2),dw2dtd=0, (27)

where F and G are defined in (18). Note that (27) is the δ→0 limit of the desingularized system (18) from Sec. 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 MSS. The FSN condition at the double singular limit can be obtained from letting δ→0 in the FSN1 condition (22) or the FSN2 condition (23). This implies that, at the double singular limit, an FSN^1^ becomes a CDH singularity [see Figs. 4(c) and 4(d)],

FSN(ε,δ)→(0,0)1:={(V1,w1,V2,w2)∈M:f2=0}={(V1,w1,V2,w2)∈Ls:f2=g1=0}=MSS∩Ls, (28)

whereas an FSN2 singularity is characterized by

FSN(ε,δ)→(0,0)2:={(V1,w1,V2,w2)∈M:Q(V1,V2,w2)=0}. (29)

According to Remarks II.2 and II.3, an FSN1 singularity from the ε viewpoint converges to a CDH as δ→0 and a DHB MSSH 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 Ref. 26. However, this is not always the case. Specifically, we find that while the CDH on the upper fold surface Ls [see Figs. 4(b) and 4(d)] supports MMOs with a high level of robustness due to the coexistence and interaction of two distinct MMO mechanisms, no MMO dynamics 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 Ref. 26.

The CDH points in system (2) are FSN1 singularities at the double singular limits. We prove in Appendix C that the CDH singularity in (2) is a novel type of folded saddle-node singularity as described in Ref. 26, with the center manifold of the CDH transverse to the fold Ls of the critical manifold. This is further confirmed in Figs. 4(c) and 4(d), as discussed in Remark II.1.

In the case of gsyn=4.3 (see the left panels of Fig. 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 Sec. III). For gsyn=4.4, an upper CDH exists. We show in Sec. IV that this CDH serves as an organizing center for the local small-amplitude 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 Ls. However, MMOs are not observed in the neighborhood of any lower CDH (see Sec. III for more discussions).

Singular orbit construction

To understand the dynamics of (2) using GSPT,^18^ 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 Subsections II A–II D to construct singular approximations of (2) in order to understand the full system trajectory Γ(ε,δ). We let Γ(0,δ)F and Γ(0,δ)S denote trajectories of the fast layer problem (7) and the slow reduced problem (12) obtained from the ε→0 singular limit. Let Γ(ε,0)S and Γ(ε,0)SS 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 Γ(0,0)x, x∈{F,S,SS} denote the fast, slow, and superslow flows at the double singular limits.

We start by showing the singular orbit construction for the (V2,w2)-subsystem (Fig. 5), which is a relaxation oscillation that is independent of gsyn, (V1,w1), and ε. The singular orbit (Fig. 5, green trajectory) consists of a superslow excursion along the left branch of V2-nullcline ( Γ(ε,0)SS, phase ◯1), a slow jump at the lower fold of V2-nullcline up to its right branch ( Γ(ε,0)S, phase ◯2), a superslow excursion through the active phase ( Γ(ε,0)SS, phase ◯3), and a slow jump back to its left branch ( Γ(ε,0)S, phase ◯4). Yellow 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 (V2,w2)-space, the V2-nullcline (red) corresponds to the superslow manifold MSS and the singular orbit Γ(ε,0)S∪Γ(ε,0)SS overlaps with the singular orbit from the double singular limits due to its independence on ε. For sufficiently small perturbation δ, Γ(ε,0)S∪Γ(ε,0)SS perturbs to the full system trajectory Γ(ε,δ) shown by the black curve.

Figure 6(a) illustrates the singular orbit Γ(ε,0)S∪Γ(ε,0)SS for gsyn=4.3 in (V1,V2,w1)-space, together with the superslow manifold MSS (red solid curve MSSa: attracting; red dashed curve MSSr: repelling). In other panels, the critical manifold MS (blue surface) and its folds Ls (blue curves) are also plotted. MS is separated into three sheets by the folds Ls, in which the upper ( MSU) and lower ( MSL) branches are stable and the middle branch ( MSM) is unstable.

Starting from the yellow star in panel (a), the singular orbit is in phase ◯1 and evolves along the lower branch of MSS under the superslow reduced problem (15). After hitting the DHB where the stability of MSS changes, the slow layer problem (13) takes over but with V2 still evolving on a superslow timescale (see Fig. 5, phase ◯1). Thus, the singular orbit during the rest of this phase is a continuum of (V1,w1) relaxation oscillations. As the evolution speed of V2 changes from superslow to slow at the yellow circle (beginning of phase ◯2), a few more spikes occur before the slow flow Γ(ε,0)S travels to the yellow square on MSS, at which phase ◯3 begins. After that, the superslow reduced problem takes over until the singular orbit Γ(ε,0)SS reaches the upper DHB at which MSS becomes unstable. As such, a family of (V1,w1) oscillations emerges as we observed during phase ◯1. As phase ◯4 begins at the yellow triangle, several additional V1 spikes occur before the slow flow Γ(ε,0)S travels back to the yellow star, thus completing a full cycle.

Figure 6(b) 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 MSS. Instead, we observe a continuum of V1 spikes throughout phase ◯3. This is because the upper MSS becomes unstable for all V2 values as ε→0, which will be discussed further in Sec. III B 1. In (b), the singular orbit consists of Γ(0,0)F (triple arrow) that are fast V1 jumps from Ls, Γ(0,0)S (double arrow) which travels along stable branches of MS or the intersection of MS and the V2-nullcline, and Γ(0,0)SS (single arrow) that follows the stable branch of MSS.

While the V2 value at the slow/superslow transition (yellow circle or yellow triangle) is uniquely determined, the corresponding (V1,w1) values can assume arbitrary positions on the upper or lower sheet of MS for that fixed V2. 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 Γ(ε,0)S 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 gsyn=4.3 in Sec. 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 MSS 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 V1 jumps and (16), which describes the slow motion along MS.

For completeness, we also plot the singular orbit for ε=0 but δ≠0 in Fig. 6(c). Instead of a continuum of V1 spikes, we observe a finite number of V1 spikes during phases ◯1 and ◯3 since δ≠0. In this limit, the fast segment ( Γ(0,δ)F, triple arrow) is the same as Γ(0,0)F, whereas the slow segments Γ(0,δ)S are O(δ) perturbations of Γ(0,0)S or Γ(0,0)SS from (b). The latter has also been illustrated in Fig. 5.

Finally, the singular orbit for gsyn=4.4 at the double singular limits is shown in Fig. 6(d). There are two major differences between gsyn=4.3 and gsyn=4.4: First, in contrast to the gsyn=4.3 case where the upper DHB vanishes, it now converges to the upper CDH at the double singular limits in (d). Second, Γ(0,0)SS in panel (d) follows the upper stable branch of MSS until reaching a saddle-node bifurcation at the yellow triangle where MSS 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 gsyn=4.3 case. Later in Sec. IV, we show how this singular orbit perturbs to MMOs at gsyn=4.4 for various perturbations.

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 small-amplitude oscillations (SAOs). We show that the only existing mechanism for MMOs at gsyn=4.3 is the delayed-Hopf bifurcation (DHB) (see Sec. III A) and explain why the absence of an upper CDH leads to the sensitivity of MMOs to timescale variations (see Secs. III B and III D). In particular, we explain the complex transitions between MMOs and non-MMOs due to changes of ε or δ in Sec. III D [see Fig. 3(a)]. We also discuss why there is no MMOs near the lower CDH in Sec. III C.

Relation of the trajectory to MS and MSS

The MMO solution of (2) for gsyn=4.3 from Fig. 1 is projected onto the (V1,w1,V2)-space in Fig. 7. The full system equilibrium (black circle) lies on MSM and is unstable. The stability of the superslow manifold MSS changes at the two DHBs (red circles). In particular, as V2 decreases, the upper branch of MSS 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 Ls always lies above the upper branch of the superslow manifold MSS so there is no upper CDH, whereas the lower fold intersects the lower branch of MSS 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 Sec. 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 Fig. 6(a) (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 Γ(ε,0)SS perturbs to Γ(ε,δ), which displays two SAOs around the stable branch of MSS. These SAOs soon transition to large-amplitude oscillations before crossing the DHB at the red circle to reach the unstable branch of MSS. As phase ◯4 begins at the green triangle, the slow jump down of V2 brings the trajectory to a region of MSL where there is a nearby stable MSS that attracts the trajectory. From the green star, the trajectory follows MSS on the superslow timescale to the lower fold of MS (see Fig. 7, lower inset), where it jumps up to MSU on the fast timescale under (7), which corresponds to the onset of spikes in V1. The spikes persist until reaching the green circle, at which phase ◯2 begins and the slow jump up of V2 brings the trajectory to the green square, completing a full cycle.

We claim that the only mechanism underlying the MMOs at gsyn=4.3 is the DHB mechanism. This is not surprising given that the upper Γ(ε,0)SS switches to a continuum of big spikes when the upper DHB vanishes at the singular limit ε=0, as shown in Fig. 6(b). To understand why canard dynamics are not involved, we view the trajectory in (V1,V2,w2)-space (see Fig. 8). Unlike Fig. 7 where Ls and curves of folded singularity ( M) overlap, the (V1,V2,w2)-projection captures the structure of the fold surfaces Ls 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 (V2,V1)-plane, which includes the periodic orbits born at the upper Hopf bifurcation (see Fig. 9).

Canard mechanism does not contribute to MMOs

In Fig. 8, the solution of (2) for gsyn=4.3 is projected onto the space (V1,V2,w2) with two separate fold surfaces Ls (blue) and two branches of folded singularities. In the lower folded singularity curve, there is an FSN2 (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 MSU. The upper folded singularity curve also has an FSN2 (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 Ls 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 as discussed inSec. II C. The bifurcation diagrams of the slow layer problem (13) are projected onto the (V2,V1)-plane, along with singular and perturbed orbits [see Figs. 9(a)–9(c)]. Both the upper and lower DHBs (red circle) are subcritical. As V2 increases, a branch of unstable periodic orbits emerges in a homoclinic bifurcation involving the middle branch of equilibria and terminates in the lower DHB (see the inset). A second branch of stable periodic orbits of large amplitudes is created in another homoclinic bifurcation, which then terminates at a saddle-node bifurcation of periodic orbits for larger V2 as it coalesces with a third family of unstable periodic solutions born in the upper DHB.

Figures 9(a)–9(c) illustrate how the δ→0 singular orbits from Fig. 6(a) perturb to MMOs for various perturbations δ. With small perturbations [Fig. 9(a), δ≈0.042], the trajectory closely follows the attracting side of MSS (denoted as MSSa), passes over the upper DHB (red circle) to the repelling side of MSS (denoted as MSSr), and undergoes a delay in which the trajectory traces MSSr before it jumps away. As the perturbation size increases, the delay is less substantial and the observed small oscillations around MSSa become larger and fewer [see Fig. 9(b)]. Panel (c) illustrates the projection of the full system trajectory Γ(ε,δ) under the default perturbation δ≈0.053. Starting from the green square, the SAOs gradually decrease in magnitude as the trajectory moves toward the upper subcritical DHB. After only 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. Note that although the default perturbation is only slightly larger than δ≈0.042 in panel (a), the delay is no longer present and the trajectory jumps away from MSSa before reaching the upper DHB. While this might initially suggest that the default perturbation δ≈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 Ref. 35, where GSPT analysis has been successfully employed to elucidate the dynamics across various gsyn values.

It is worth emphasizing two interesting First, 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 MSSa, a detailed analysis of which is provided in Sec. III D. Therefore, the absence of a delay in Fig. 9(c) does not imply a significant deviation from the singular limit. Rather, it is mainly due to the manner in which the trajectory approaches MSSa 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 Fig. 13), a slight reduction or increase in the perturbation size can induce a delay phenomenon. Second, 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 V2 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 MSSr is significantly shorter than that near MSSa.

In Subsections III B–III D, 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.

Effects of varying

ε and δ on DHB and FSN points

When gsyn=4.3, the CDH singularity only exists on the lower fold surface of MS. We demonstrate that this leads to different effects of ε on the upper and lower DHB points, respectively (see Fig. 10). We also examine the effect of δ on the lower FSN points (see Fig. 11).

Effect of

ε on DHB

The effects of ε on the DHBs MSSH are summarized by the two-parameter bifurcation diagrams of (13) projected onto (V2,C1)-space (see Fig. 10). As C1 decreases (or equivalently, as ε decreases), the upper Hopf moves to larger V2 values and eventually vanishes for ε small enough [see Fig. 10(a)]. This explains why the upper singular orbit Γ(ε,0)SS in Fig. 6(a) for ε≠0 switches to a continuum of spikes in Fig. 6(b) as ε→0. On the other hand, since there is a CDH on the lower Ls, the lower Hopf will converge to that CDH as ε→0 [see Fig. 10(b) and recall Remark II.3]. When ε increases, the lower Hopf and the lower fold of MSS 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 MSS [also see Figs. 9(a) and 9(c), lower magenta circle].

Effect of δ on FSN points

There is no CDH or FSN1 on the upper fold, hence, we only examine the effect of δ on FSN singularities on the lower Ls. Figure 11(a) shows as δ→0 (or equivalently, ϕ2→0), the FSN1 singularity converges to the lower CDH as demonstrated in our analysis [see (28)]. On the other hand, the FSN2 singularity is significantly distant from the CDH [see Fig. 11(b) and recall the condition (29)].

Why there are no SAOs near the lower CDH

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 FSN1 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 canard-induced 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 MSS 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 Ls and eventually vanish upon crossing the actual fold of MSS [see Fig. 10(b)]. 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 Fig. 10(b)] and reducing δ pushes the trajectory to travel along MSS 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 MSS, i.e., close to a double zero eigenvalue at a BT bifurcation of the (V1,w1,V2) subsystem [see Figs. 9 and 10(b)]. As a result, the branch of unstable small-amplitude periodic orbits born at the lower HB is almost invisible [see the inset of Figs. 9(a) and 9(c)] and there is only a small region of MSS along which the Jacobian matrix JSL (24) of the slow layer problem (13) has complex eigenvalues. Figure 9(d) shows the real part of the first two eigenvalues λ of (24) evaluated along MSS, excluding the third eigenvalue given by f2V2. 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 (d), the eigenvalues on the stable lower branch of MSS are initially real and negative. That is, the trajectory approaches the attracting MSS along stable nodes of the slow layer problem. As the superslow flow brings the trajectory toward 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.

Effects of varying

ε and δ on MMOs

Recall when gsyn=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 C1 and vary δ by changing ϕ2. Recall that increasing (respectively, decreasing) C1 or ε slows down (respectively, speeds up) the fast variable V1, whereas increasing (respectively, decreasing) ϕ2 or δ speeds up (respectively, slows down) the superslow variable w2. Other (slow) variables are not affected.

Figure 3(a) summarizes the effects of (C1,ϕ2) on MMOs when gsyn=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

∙For fixed δ=0.053 at ϕ2=0.001, slowing down the fast variable V1 by increasing C1 from its default value C1=8 [Fig. 3(a), vertical line above the red star] preserves the MMOs [also see Fig. 12(a), 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 V1 via decreasing C1 leads to a total of three transitions between MMOs and non-MMOs [see Fig. 12(b)]. These transitions correspond to the three crossings between the vertical black line and the yellow/blue boundary in Fig. 3(a). ∙For fixed ε=0.1 at C1=8, slowing down the superslow variable w2 by decreasing ϕ2 from the default value ϕ2=0.001 [Fig. 3(a), 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 non-monotonic behavior [see Fig. 13(a)]. ∙For fixed ε at C1=8, speeding up the superslow variable w2 by increasing ϕ2 leads to a total of five transitions between MMOs and non-MMOs [see Fig. 13(b)]. These transitions correspond to the five crossings between the horizontal black line and the yellow/blue boundary in Fig. 3(a).

Next, we discuss the above four scenarios separately.

MMOs are robust to increasing C1

Increasing C1 slows down the fast variable V1 and hence moves the three timescale (1F, 2S, 1SS) problem closer to (3S, 1SS) separation. As a result, the critical manifold MS and the folded singularities become less relevant with increased C1. Nonetheless, this does not affect the existence of MMOs, as they arise from the upper Hopf bifurcation. As C1 is increased, the upper HB moves to smaller V2 values [see Fig. 10(a)]. This change allows the trajectory to travel a longer distance along the stable part of MSS during phase ◯3 and generate more SAOs, as shown in Fig. 14. To simplify the presentation, we omit singular orbits and focus only on MSS 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 Fig. 9(c)]. Increasing C1 causes the upper HB in (V2,V1)-space (Fig. 14, red circle) to move to the left and become further away from the maximum of V2 (green square), which remains unchanged. As a result, the trajectory with larger C1 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 V2 [compare Figs. 9(c) and 14]. As C1 is increased to 11, the trajectory passes over the HB to MSSr and experiences a delay before jumping away (see Fig. 14, 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 V2 at the green triangle.

Decreasing C1 leads to three MMOs/non-MMOs transitions

Contrary to the preservation of MMOs with increasing C1, MMOs appear to be sensitive to the decrease of C1. Specifically, speeding up the fast variable V1 by decreasing C1 leads to transitions from MMOs (e.g., at C1=8) to non-MMOs (e.g., C1=7), back to MMOs at C1∈(1.46,4.2), and then to non-MMOs for C1<1.46.

The initial decrease of C1 [e.g., from C1=8 to C1=7, see Fig. 15(a)] 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 C1. Interestingly, a further decrease of C1 to a range of C1∈(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 Figs. 15(b) and 15(c)]. This is because, for C1∈(1.46,4.2) and V2 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 MSS. When C1 is reduced to be below 1.46, the voltage V1 exhibits rapid spikes that occur immediately after the green square, consistent with the singular orbits shown in Figs. 6(b) and 6(c). As a result, there is no more MMOs [see Fig. 15(d), C1=1.4].

To better understand why the SAOs for C1∈(1.46,4.2) grow in amplitude [see Figs. 15(b) and 15(c)], we project the trajectory when C1=2 onto (V1,w1,V2) space (see Fig. 16). The blue triangle denotes a saddle-focus of the (V1,w1,V2) subsystem near the green square. During phase ◯2 from the green circle to the green square, the trajectory travels toward the saddle-focus (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 upward 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^48^ and a singular Hopf bifurcation in two-timescale settings.^1^

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 Fig. 13(a), where the amplitude and the number of SAOs exhibit an alternation between increase and decrease.

Smaller perturbations should cause the solution trajectory to follow MSS 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 Fig. 13(a) (top three rows) and also in Fig. 9 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 Fig. 13(a), 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 Fig. 13(a)].

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 MSS 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 MSSa 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 MSSa 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.

Increasing ϕ2 leads to five MMOs/non-MMOs transitions

When ϕ2 increases, the superslow variable w2 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(ε13)]. 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 Fig. 13(b), 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 Figs. 13(a) and 13(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.

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 Sec. IV A that both canard and DHB mechanisms coexist and interact to produce MMOs that exhibit significant robustness to timescale variations as shown in Fig. 3(b). We explain the robust occurrence of MMOs in Sec. 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.

Relation of the trajectory to MS and MSS

The solution of (2) for gsyn=4.4 in Fig. 1 is projected onto the (V1,w1,V2)-space, together with the singular orbit from Fig. 6(d) (green curve), critical manifold MS (blue surface), fold Ls (blue curve), and the superslow manifold MSS (red curves) (see Fig. 17). As the coupling strength gsyn increases from 4.3 to 4.4, two new features regarding the upper MSS emerge. First, the stability of the upper MSS changes at a fold point Lss1 (the green triangle) rather than a DHB. Specifically, as V2 increases, the equilibrium along the upper MSS 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 Ls now intersects the upper branch of MSS at a CDH singularity (the blue diamond), which is a folded node at the default parameter values given in Table I. As proved in Sec. 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 Fig. 18, which shows that the upper FSN1 and HB point converge to the same CDH point on the upper Ls in the double singular limits (ε,δ)→(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 Fig. 17). During phase ◯3, the singular orbit (green) traces MSSa 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 MSL. 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 Figs. 22 and 23), the delay is more substantial, but the small oscillations are very small and difficult to observe due to the stronger attraction to MSSa during phase ◯3. We omit explaining the solution dynamics of (2) for gsyn=4.4 during other phases as they are similar to those observed when gsyn=4.3. In particular, the absence of SAOs near the lower CDH point is due to the same mechanism as discussed in Sec. III C.

![FIG. 22.: Time traces of the solutions of (2) (left panels) and their projections onto the corresponding bifurcation diagrams in (V2,V1)-space (right panels) for gsyn=4.4 and (a) C1=0.8

(ε=0.01), (b) C1=40

(ε=0.5), and (c) C1=80 ( ε=1). Increasing C1 moves the upper DHB (red circle) to the lower V2 and manifests DHB-like features. Color coding and symbols of the bifurcation diagrams on the right panel are the same as in Fig. 17.](CHAOEH-000034-053119_1-g022.jpg)

At default parameters, δ=O(ε). The emergence of SAOs for gsyn=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 Ls, as seen in examples considered in Refs. 26 and 27. To understand why folded node singularities play an important role in the occurrence of SAOs, we view the trajectory and folded singularities projected onto (V1,V2,w2)-space, as illustrated in Fig. 19, as well as draw the funnel volume associated with the folded node singularity curve on the upper fold (see Fig. 20). On the other hand, to grasp the importance of the DHB toward the generation of SAOs, we examine the trajectory on the (V1,V2) projection, which includes the periodic orbit branches born at the upper Hopf bifurcation (see Fig. 21). For simplicity, we omit depicting singular orbits in the following figures and focus only on MS, MSS and their relations to the full trajectory.

MMOs are organized by both canard and DHB mechanisms

Canard dynamics

Figure 19 shows the solution trajectory of (2) projected onto (V1,V2,w2)-space when gsyn=4.4. The upper folded singularity curve M comprises folded singularities of two folded nodes (solid green) and folded saddles (dashed green). The blue star denotes FSN1, which occurs O(δ) away from the upper CDH at the blue diamond. Different from the previous case ( gsyn=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 Ls (see Fig. 20). As discussed in Sec. II B, the funnel for a folded node point represents a two-dimensional trapping region on MS. The 2D funnels for all folded nodes on the upper fold together form a three-dimensional funnel volume. In the (V1,V2,w2) projection (see the top panel of Fig. 20), the funnel volume is bounded by the singular strong canard surface (shown in magenta) and the upper fold surface (in blue). The bottom panel shows the projection onto the (V1,V2)-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 Ls at a regular jump point and there are no SAOs. These observations suggest that the MMOs for gsyn=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 21 shows the projection of the solution trajectory and the bifurcation diagram of the slow layer problem (13) onto the (V2,V1)-plane. Starting at the green square, the trajectory exhibits SAOs as it follows the upper branch of MSS toward the left. As the trajectory passes through the attracting region of MSS, the oscillation amplitude decays in a typical DHB fashion. Moreover, the orbit experiences a delay along the repelling MSS for an amount of time as it passes through the HB. It is worth noting that after the trajectory enters the unstable part of MSS, there are no symmetric oscillations with respect to the DHB, i.e., there are no SAOs with increasing amplitudes. This is due to the fact that V2 switches from a superslow to a slow timescale at the green triangle, as discussed earlier when gsyn=4.3.

Thus, for δ=O(ε), the MMOs for gsyn=4.4 exhibit characteristics of both the canard and DHB mechanisms. To further confirm this, we performed two perturbations on the system. First, we increased ε by raising the value of C1 to make the canard mechanism less relevant. This is because increasing C1 slows down the evolution speed of V1, which in turn drives the three-timescale system (2) closer to (3S, 1SS) splitting. We observed that MMOs persisted for C1 as large as 80 ( ε=O(1)), at which the folded singularities were no longer relevant (see Fig. 22). Second, 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.

Effects of varying

ε and δ on MMOs

Unlike the sensitivity of MMOs at gsyn=4.3 to timescale variations, the interaction of canard and DHB mechanisms due to the existence of the upper CDH when gsyn=4.4 makes MMOs much more robust, as discussed above and illustrated in Fig. 3(b). Specifically, MMOs persist over biologically relevant ranges of C1∈(0.1,80) and ϕ2∈(10−4,0.01) (also see Figs. 22 and 23).

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 Figs. 22(a) and 23(a)]. 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 DHB-like features, whereas by tuning δ≥O(ε), DHB-like features diminish and the canard mechanism dominates.

Below we summarize the effects of increasing the singular

∙For fixed ϕ2=0.001, increasing ε (i.e., increasing C1) enhances the DHB-like features of the SAOs (see Fig. 22). ∙For fixed C1=8, increasing δ (i.e., increasing ϕ2) makes the canard mechanism dominant and causes DHB-like features to gradually vanish (see Fig. 23).

Increasing C1 makes DHB dominate

Increasing C1 drives the three-timescale system (2) closer to (3S, 1SS) splitting. As a result, the critical manifold Ms and folded singularities including the CDH point become less meaningful and eventually irrelevant for sufficiently large C1 (e.g., C1=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 C1 increases, while the influence of the canard mechanism becomes less significant.

Figure 22 shows the effect of increasing C1 on voltage traces of the full system and the bifurcation diagrams of the fast subsystem (13). As C1 increases from panel (A2) to panel (C2), the upper DHB (red circle) moves away from the upper CDH (blue diamond) to smaller V2 values, whereas the CDH points and slow/superslow timescale transitions denoted by magenta and green symbols all remain unaffected by C1. As a result, the trajectory with larger C1 begins small oscillations at a larger distance from the Hopf bifurcation point, similar to what we observed in the case of gsyn=4.3. Furthermore, we have noticed that trajectories for larger C1 exhibit more pronounced DHB-like characteristics, including SAOs with decreasing amplitudes and a more extended travel distance along the unstable branch of MSS. As discussed before, due to a switch of the V2 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 w2 and hence drives the three-timescale system (2) closer to (1F, 3S) splitting. As a result, the superslow manifold MSS 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 23 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 MSS 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., Figs. 19 and 23(a)] exhibit both canard and DHB characteristics. When δ≈12ε [Fig. 23(b)], there are still some DHB-like features. Further increasing δ to δ=0.006≈ε or δ=0.01≈ε14 [Figs. 23(c) and 23(d)], the trajectories no longer closely follow MSS 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 MSS closer together and amplifies the DHB characteristics of the sustained MMOs (see Fig. 23).

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.^1,19–24,49,50^ In contrast, progress on MMOs in three-timescale problems has been made only in the recent past (see, e.g., Refs. 30, 27, 26, 41, 11, 39, 38, and 40). In this work, 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 ( gsyn=4.3 and 4.4mS/cm2). Applying geometric singular perturbation theory and bifurcation analysis,^18,35^ we have revealed that the two MMOs exhibit different mechanisms, leading to remarkably different sensitivities to variations in timescales (see Fig. 3). Specifically, when gsyn=4.3, only the δ→0 singular limit provides a reliable prediction for the observed MMOs, demonstrating that the fast subsystem delayed Hopf (DHB) is the only MMO mechanism. Conversely, for gsyn=4.4, both ε→0 and δ→0 singular limits yield faithful predictions for the MMOs. Therefore, both the canard and DHB mechanisms exist and interact to produce MMOs that are significantly more robust than MMOs with only a single mechanism.

The existence of three distinct timescales leads to two important the critical (or slow) manifold MS and the superslow manifold MSS. The point where a fold of the critical manifold Ls intersects the superslow manifold MSS is referred to as the canard-delayed-Hopf (CDH) singularity, which naturally arises in problems that involve three different timescales. Reference 26 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.^26,27^

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 CDH singularity as documented in Ref. 26 with its center manifold transverse to Ls at the CDH point (see our proof in Appendix C), SAOs in our system are not guaranteed to occur in the neighborhood of a CDH singularity. Specifically, we observed that CDH singularities on the lower Ls did not support SAOs in either the case of gsyn=4.3 or gsyn=4.4. We have explained in Sec. III C why neither of the canard and DHB mechanisms near the lower CDH gives rise to SAOs. Our analysis suggests that the absence of SAOs near the lower CDH might be due to the proximity of this CDH to the actual fold of the superslow manifold MSS. Further analytical work is still required to confirm this observation and should be considered for future work. Second, we have explored the conditions underlying the robust occurrence of MMOs in a three-timescale setting. Our analysis has revealed that the existence of CDH singularities critically determines whether or not the two MMO mechanisms (canard and DHB) can coexist and interact, which in turn greatly impacts the robustness of MMOs. In particular, we have found that MMOs occurring near a CDH singularity are much more robust than MMOs when CDH singularities are absent.

It is worth highlighting that the existence of MMOs in (2), even in the presence of an upper CDH singularity, is still contingent upon the trajectory closely approaching the CDH, as discussed in Sec. I and further demonstrated in Sec. IV (see also Fig. 19). Although the upper CDH persists for gsyn>4.3212 (see Fig. 2), increasing gsyn causes a shift of the CDH toward lower V2 values. This shift leads to numerous transitions between MMOs and non-MMOs due to the changing distance between the trajectory and the CDH near the end of large V1 spikes (data not shown). With higher gsyn values, the transition from V1 spikes to V1 plateau or SAOs can happen during phase ◯1 instead of phase ◯3. In this case, whether the trajectory can approach the CDH to produce MMOs is influenced by a mechanism analogous to the spike-adding like phenomenon that governs the transitions between MMOs and non-MMOs as δ increases when gsyn=4.3 (see Sec. III). Ultimately, a sufficiently large value of gsyn will bring the upper CDH below the minimal V2 of the trajectory, resulting in a complete loss of MMOs. A comprehensive investigation of these complex transitions induced by variations in gsyn is left for future work.

In addition to uncovering the relationship between the upper CDH singularity and robustness of MMOs, we have also provided a detailed investigation on how the features and mechanisms of MMOs without or with CDH singularities vary with respect to timescale variations. Table II outlines a summary of the different mechanisms and robustness properties of MMOs as we vary timescales. When gsyn=4.3, where no CDH was found near the small oscillations, the only mechanism for the MMOs is the DHB mechanism as we justified in Sec. III. In this case, speeding up the fast variable V1 via decreasing ε led to a total of three transitions between MMOs and non-MMO solutions due to its effect on the upper DHB point [Fig. 10(a)]. Initially, MMOs disappeared as the decrease of ε moved the DHB closer to the green square where the SAO phase began, resulting in insufficient time for generating small oscillations [see Fig. 15(a)]. However, 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 MSS [see Figs. 15(b), 15(c), and 16]. Eventually, SAOs disappeared entirely when V1 became so rapid that it failed to remain in proximity of MSS to generate small oscillations.

As one would expect, MMOs when gsyn=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 Figs. 3 and 13(b)]. 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 gsyn=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 Figs. 12(a) and 14]. 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 Fig. 13(a)]. 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 gsyn=4.3, an upper CDH occurs near the SAOs at gsyn=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 gsyn=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 gsyn=4.3, MMOs at gsyn=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 gsyn=4.3 [Fig. 10(a)], decreasing ε at gsyn=4.4 brings the DHB point closer to the CDH (Fig. 18, right panel) and, consequently, closer to the midst of small oscillations (Fig. 22). Moreover, this timescale change manifests the canard mechanism by moving the system closer to the ε→0 singular limit. Hence, MMOs at gsyn=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 gsyn=4.4 due to the presence of a nearby upper CDH, leads to significantly enhanced robustness against timescale variations compared with gsyn=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 gsyn=4.4, unlike what we observed in gsyn=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.