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.

## I. 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.

^{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 $\epsilon $ and $\delta $, 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.

^{43,44}that was introduced by Ref. 35. The model equations are given by

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., Refs. 32 and 33). For the physiological description of functions in (2) and (3), we refer readers to Ref. 35 for details.

Parameter values . | |||||
---|---|---|---|---|---|

C_{1} | 8 μF/cm^{2} | I_{1} | 0 μA/cm^{2} | ϕ_{1} | 0.01 |

C_{2} | 100 μF/cm^{2} | I_{2} | 60 μA/cm^{2} | ϕ_{2} | 0.001 |

V_{Ca} | 120 mV | g_{Ca} | 4 mS/cm^{2} | K_{1} | −1.2 mV |

V_{K} | −84 mV | g_{K} | 8 mS/cm^{2} | K_{2} | 18 mV |

V_{L} | −60 mV | g_{L} | 2 mS/cm^{2} | K_{3} | 12 mV |

V_{syn} | 30 mV | θ_{s} | −20 mV | K_{4} | 17.4 mV |

β | 0.5 ms^{−1} | σ_{s} | 10 mV |

Parameter values . | |||||
---|---|---|---|---|---|

C_{1} | 8 μF/cm^{2} | I_{1} | 0 μA/cm^{2} | ϕ_{1} | 0.01 |

C_{2} | 100 μF/cm^{2} | I_{2} | 60 μA/cm^{2} | ϕ_{2} | 0.001 |

V_{Ca} | 120 mV | g_{Ca} | 4 mS/cm^{2} | K_{1} | −1.2 mV |

V_{K} | −84 mV | g_{K} | 8 mS/cm^{2} | K_{2} | 18 mV |

V_{L} | −60 mV | g_{L} | 2 mS/cm^{2} | K_{3} | 12 mV |

V_{syn} | 30 mV | θ_{s} | −20 mV | K_{4} | 17.4 mV |

β | 0.5 ms^{−1} | σ_{s} | 10 mV |

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 s y n$ and the dynamics of $( V 1, w 1$). 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 $ 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. 35 (e.g., $ g syn=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 $ g syn=4.3$ and $ g syn=4.4$ (with the unit of $ mS / cm 2$), as highlighted in Fig. 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 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 $ V i$ values, $i\u2208{1,2}$ (denoted as upper CDH) and the other at smaller $ V i$ (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 $( 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 and facilitate the existence of MMOs. When $ g syn<4.2628$ (e.g., $ g syn=1$ and $4.1$ as considered in Ref. 35), the system (2) does not produce MMOs. When $ g syn>4.3213$, MMOs may arise through mechanisms similar to those we will thoroughly examine for $ g syn=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 $ g syn=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 $ V 1$ spikes to a $ V 1$ plateau, as detailed in Ref. 35. Therefore, as $ g syn$ 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 $ g syn$ values (e.g., $ g syn=80$), where both CDH and DHB points fall to low $ V 2$ values that lie beyond the system’s trajectory.

^{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 problem:

*slow timescale*the

*slow system*in which $ V 1$ evolves on a timescale of $O( \epsilon \u2212 1)$, $( w 1, V 2)$ on $O(1)$ and $ w 2$ on $O(\delta )$. Introducing a superslow time $ t s s=\delta t s$ yields an equivalent description of dynamics,

*superslow timescale*and is called the

*superslow system*. Alternatively, defining a fast time $ t f= t s/\epsilon $, we obtain the following

*fast system*:

*fast timescale*.

The presence of two independent singular perturbation parameters, $\epsilon $ and $\delta $, 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 $ g syn=4.3$ and $ g syn=4.4$ arise from distinct mechanisms, resulting in remarkably different sensitivities to variations in timescales, as illustrated in Fig. 3.

When $ g syn=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 $\delta \u21920$ singular limit provides a faithful prediction for the observed MMOs, whereas the $\epsilon \u21920$ limit does not. This observation pinpoints the DHB from the $\delta \u21920$ limit as the only mechanism for the MMOs at $ g syn=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 $\delta \u2264 O(\epsilon )$. Increasing $\delta $ via increasing $ \varphi 2$ or decreasing $\epsilon $ via decreasing $ C 1$ to a degree where $\delta > O(\epsilon )$ leads to multiple MMOs/non-MMOs transitions. MMOs are completely lost when $\delta \u2265 O( \epsilon 1 3)$ for which the DHB is no longer present.

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 co-modulate properties of the local oscillatory behavior. In this case, both the $\epsilon \u21920$ and $\delta \u21920$ singular limits contribute faithful predictions for the observed dynamics, resulting in MMOs with significantly stronger robustness than $ g syn=4.3$ [Fig. 3(b)]. We demonstrate that when $\delta = O(\epsilon )$, MMOs exhibit both canard and DHB features. Upon tuning $\delta \u2265O( \epsilon )$, 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 $\epsilon $ and $\delta $.

The paper is organized as follows. In Sec. II, we perform a geometric singular perturbation analysis of the 3-timescale problem (2) by treating $\epsilon $ as the principal perturbation parameter while keeping $\delta $ fixed, treating $\delta $ as the principal perturbation parameter while keeping $\epsilon $ fixed, and by treating $\epsilon $ and $\delta $ as two independent perturbation parameters. We review both mechanisms for MMOs and discuss their interaction at the double singular limit $(\epsilon ,\delta )\u2192(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 $ g syn=4.3$ to variations in perturbation sizes $\epsilon $ or $\delta $ [i.e., varying $ C 1$ and $ \varphi 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 $\epsilon $ and $\delta $ 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 $ 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 the 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 Fig. 3(b)]. Finally, we conclude in Sec. V with a discussion.

## II. 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 $\epsilon $ as the only singular perturbation parameter,^{19,20} treating $\delta $ as the only singular perturbation parameter,^{21–24} and finally treating $\epsilon $ and $\delta $ 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.

### A. Singular limits

#### 1. $\epsilon \u21920$ Singular limit

*fast layer problem*, a system that describes the dynamics of the fast variable, $ V 1$, for fixed values of the other variables,

*critical manifold*and is denoted as $ M S$,

*slow manifold*, on which $ w 1$ is given by an $O(\epsilon )$-perturbation of $ F 1$ (Ref. 18); we simply use $ M S$ as a convenient numerical approximation of these slow manifolds.

*slow reduced problem*, a system that describes the dynamics of $ w 1, V 2, w 2$ along $ M S$,

#### 2. *δ* → 0 Singular limit

*δ*→

*slow layer problem*in the form

*superslow manifold*and is denoted as $ M S S$,

#### 3. $\epsilon \u21920,\delta \u21920$ Double singular limits

*slow reduced layer problem*

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 Ref. 35.

### B. Slow reduced problem and canard dynamics

^{19}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 ( $ \lambda w, \lambda s$) where $ | \lambda w |< | \lambda 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> \lambda w> \lambda 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.

#### 1. Folded saddle node (FSN)

^{26}

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) that mark the boundary between folded nodes and folded saddles are neither of type II nor type III.^{26,45–47} They are also not type I because the center manifolds associated with $ FSN 1$ and $ FSN 2$ are transverse to the fold surface at the singularity points [see Figs. 4(c) and 4(d), yellow and pink planes]. Following Ref. 26, we identify $ FSN 1$ and $ FSN 2$ as novel types of $ FSN$.

The condition (22) suggests that an $ FSN 1$ is $O(\delta )$ close to the intersection point of the superslow manifold $ M s s$ and the fold surface $ L s$, which was defined as the canard-delayed-Hopf (CDH) singularity in Ref. 26. 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 $\delta \u22600$ (top panels) and the singular limit $\delta =0$ (bottom panels). Recall the bifurcation of the upper CDH with respect to $ g syn$ is shown in Fig. 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 $\delta \u21920$ singular limit, which exhibits delayed Hopf bifurcations that allow for interesting dynamics.

It follows from (25) that an $ M S S H$ is $O(\epsilon )$ close to the intersection of $ M S S$ and $ L s$, i.e., a CDH singularity. The subsystem HB bifurcation $ M S S H$ is also known as delayed Hopf bifurcation (DHB).

### D. Interaction between canard and delayed Hopf

According to Remarks II.2 and II.3, an $ FSN 1$ singularity from the $\epsilon $ viewpoint converges to a CDH as $\delta \u21920$ and a DHB $ M S S H$ from the $\delta $ viewpoint converges to a CDH as $\epsilon \u21920$. 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 $ L s$ [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 $ FSN 1$ 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 $ L s$ 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 $ g s y n=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 $ g syn=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 $ L s$. However, MMOs are not observed in the neighborhood of any lower CDH (see Sec. III for more discussions).

### E. 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 $ \Gamma ( \epsilon , \delta )$. We let $ \Gamma ( 0 , \delta ) F$ and $ \Gamma ( 0 , \delta ) S$ denote trajectories of the fast layer problem (7) and the slow reduced problem (12) obtained from the $\epsilon \u21920$ singular limit. Let $ \Gamma ( \epsilon , 0 ) S$ and $ \Gamma ( \epsilon , 0 ) S S$ denote trajectories of the slow layer problem (13) and the superslow reduced problem (15) obtained from the $\delta \u21920$ 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 $ \Gamma ( 0 , 0 ) x$, $x\u2208{ 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 $( V 2, w 2)$-subsystem (Fig. 5), which is a relaxation oscillation that is independent of $ g syn$, $( V 1, w 1)$, and $\epsilon $. The singular orbit (Fig. 5, green trajectory) consists of a superslow excursion along the left branch of $ V 2$-nullcline ( $ \Gamma ( \epsilon , 0 ) S S$, phase $\u25ef 1$), a slow jump at the lower fold of $ V 2$-nullcline up to its right branch ( $ \Gamma ( \epsilon , 0 ) S$, phase $\u25ef 2$), a superslow excursion through the active phase ( $ \Gamma ( \epsilon , 0 ) S S$, phase $\u25ef 3$), and a slow jump back to its left branch ( $ \Gamma ( \epsilon , 0 ) S$, phase $\u25ef 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 $( V 2, w 2)$-space, the $ V 2$-nullcline (red) corresponds to the superslow manifold $ M S S$ and the singular orbit $ \Gamma ( \epsilon , 0 ) S\u222a \Gamma ( \epsilon , 0 ) S S$ overlaps with the singular orbit from the double singular limits due to its independence on $\epsilon $. For sufficiently small perturbation $\delta $, $ \Gamma ( \epsilon , 0 ) S\u222a \Gamma ( \epsilon , 0 ) S S$ perturbs to the full system trajectory $ \Gamma ( \epsilon , \delta )$ shown by the black curve.

Figure 6(a) illustrates the singular orbit $ \Gamma ( \epsilon , 0 ) S\u222a \Gamma ( \epsilon , 0 ) S S$ for $ g syn=4.3$ in $( V 1, V 2, w 1)$-space, together with the superslow manifold $ M S S$ (red solid curve $ M S S a$: attracting; red dashed curve $ M S S r$: 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 S U$) and lower ( $ M S L$) branches are stable and the middle branch ( $ M S M$) is unstable.

Starting from the yellow star in panel (a), the singular orbit is in phase $\u25ef 1$ and evolves along the lower branch of $ M S S$ under the superslow reduced problem (15). After hitting the DHB where the stability of $ M S S$ changes, the slow layer problem (13) takes over but with $ V 2$ still evolving on a superslow timescale (see Fig. 5, phase $\u25ef 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 yellow circle (beginning of phase $\u25ef 2$), a few more spikes occur before the slow flow $ \Gamma ( \epsilon , 0 ) S$ travels to the yellow square on $ M S S$, at which phase $\u25ef 3$ begins. After that, the superslow reduced problem takes over until the singular orbit $ \Gamma ( \epsilon , 0 ) S S$ reaches the upper DHB at which $ M S S$ becomes unstable. As such, a family of $( V 1, w 1)$ oscillations emerges as we observed during phase $\u25ef 1$. As phase $\u25ef 4$ begins at the yellow triangle, several additional $ V 1$ spikes occur before the slow flow $ \Gamma ( \epsilon , 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 $ M S S$. Instead, we observe a continuum of $ V 1$ spikes throughout phase $\u25ef 3$. This is because the upper $ M S S$ becomes unstable for all $ V 2$ values as $\epsilon \u21920$, which will be discussed further in Sec. III B 1. In (b), the singular orbit consists of $ \Gamma ( 0 , 0 ) F$ (triple arrow) that are fast $ V 1$ jumps from $ L s$, $ \Gamma ( 0 , 0 ) S$ (double arrow) which travels along stable branches of $ M S$ or the intersection of $ M S$ and the $ V 2$-nullcline, and $ \Gamma ( 0 , 0 ) S S$ (single arrow) that follows the stable branch of $ M S S$.

While the $ V 2$ value at the slow/superslow transition (yellow circle or yellow 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 $\u25ef 2$ and $\u25ef 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 $ \Gamma ( \epsilon , 0 ) S$ during phases $\u25ef 2$ and $\u25ef 4$.

Comparing (a) and (b) indicates that to obtain MMO dynamics in the perturbed system, $\epsilon $ cannot be too small. To explain the MMO dynamics for $ g syn=4.3$ in Sec. III, we mainly make use of the singular orbit in panel (a) with $0<\epsilon \u226a1$. This orbit, aside from the upper superslow segment where the stability of $ M S S$ differs between the two panels, can be viewed as an $ O(\epsilon )$ 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$.

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

Finally, the singular orbit for $ g syn=4.4$ at the double singular limits is shown in Fig. 6(d). There are two major differences between $ g syn=4.3$ and $ g syn=4.4$: First, in contrast to the $ g syn=4.3$ case where the upper DHB vanishes, it now converges to the upper CDH at the double singular limits in (d). Second, $ \Gamma ( 0 , 0 ) S S$ in panel (d) follows the upper stable branch of $ M S S$ until reaching a saddle-node bifurcation at the yellow triangle where $ M S S$ becomes unstable. This results in a unique construction of singular orbit segment during phase $\u25ef 4$ due to its unique starting position. Nonetheless, the construction of singular orbits during phase $\u25ef 2$ remains non-unique as discussed in the $ g syn=4.3$ case. Later in Sec. IV, we show how this singular orbit perturbs to MMOs at $ g syn=4.4$ for various perturbations.

## III. ANALYSIS OF MMOs WHEN *g*_{syn} **=** 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 $ g syn=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 $\epsilon $ or $\delta $ in Sec. III D [see Fig. 3(a)]. We also discuss why there is no MMOs near the lower CDH in Sec. III C.

### A. Relation of the trajectory to *M*_{S} and *M*_{SS}

The MMO solution of (2) for $ g syn=4.3$ from Fig. 1 is projected onto the $( V 1, w 1, V 2)$-space in Fig. 7. The full system equilibrium (black circle) lies on $ M S M$ and is unstable. The stability of the superslow manifold $ M S S$ changes at the two DHBs (red circles). In particular, as $ V 2$ decreases, the upper branch of $ M S S$ 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 S S$ so there is no upper CDH, whereas the lower fold intersects the lower branch of $ M S S$ at a CDH singularity (blue diamond in the lower inset). This CDH is a folded focus and will become an FSN $ 1$ as $\delta \u21920$ as discussed in Sec. II D. Moreover, the nearby HB bifurcation (red circle in the lower inset) is $O(\epsilon )$ 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 $\u25ef 3$, the upper superslow segment $ \Gamma ( \epsilon , 0 ) S S$ perturbs to $ \Gamma ( \epsilon , \delta )$, which displays two SAOs around the stable branch of $ M S S$. These SAOs soon transition to large-amplitude oscillations before crossing the DHB at the red circle to reach the unstable branch of $ M S S$. As phase $\u25ef 4$ begins at the green triangle, the slow jump down of $ V 2$ brings the trajectory to a region of $ M S L$ where there is a nearby stable $ M S S$ that attracts the trajectory. From the green star, the trajectory follows $ M S S$ on the superslow timescale to the lower fold of $ M S$ (see Fig. 7, lower inset), where it jumps up to $ M S U$ 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 $\u25ef 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 $ \Gamma ( \epsilon , 0 ) S S$ switches to a continuum of big spikes when the upper DHB vanishes at the singular limit $\epsilon =0$, as shown in Fig. 6(b). To understand why canard dynamics are not involved, we view the trajectory in $( V 1, V 2, w 2)$-space (see Fig. 8). Unlike Fig. 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 Fig. 9).

#### 1. Canard mechanism does not contribute to MMOs

In Fig. 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 S U$. 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.

#### 2. 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 $\delta $ as the main singular perturbation parameter as discussed inSec. II C. The bifurcation diagrams of the slow layer problem (13) are projected onto the $( V 2, V 1)$-plane, along with singular and perturbed orbits [see Figs. 9(a)–9(c)]. Both the upper and lower DHBs (red circle) are subcritical. As $ V 2$ 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 $ V 2$ as it coalesces with a third family of unstable periodic solutions born in the upper DHB.

Figures 9(a)–9(c) illustrate how the $\delta \u21920$ singular orbits from Fig. 6(a) perturb to MMOs for various perturbations $\delta $. With small perturbations [Fig. 9(a), $\delta \u22480.042$], the trajectory closely follows the attracting side of $ M S S$ (denoted as $ M S S a$), passes over the upper DHB (red circle) to the repelling side of $ M S S$ (denoted as $ M S S r$), and undergoes a delay in which the trajectory traces $ M S S r$ before it jumps away. As the perturbation size increases, the delay is less substantial and the observed small oscillations around $ M S S a$ become larger and fewer [see Fig. 9(b)]. Panel (c) illustrates the projection of the full system trajectory $ \Gamma ( \epsilon , \delta )$ under the default perturbation $\delta \u22480.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 $\delta \u22480.042$ in panel (a), the delay is no longer present and the trajectory jumps away from $ M S S a$ before reaching the upper DHB. While this might initially suggest that the default perturbation $\delta \u22480.053$ is distant from the singular limit, we show below that this is not the case. Moreover, this specific value of $\delta $ aligns with the perturbation used in Ref. 35, where GSPT analysis has been successfully employed to elucidate the dynamics across various $ g syn$ values.

It is worth emphasizing two interesting points: First, the impact of decreasing $\delta $ 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 $\u25ef 2$ approaches $ M S S a$, 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 $ M S S a$ during phase $\u25ef 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 $ 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 S S r$ is significantly shorter than that near $ M S S a$.

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 $\epsilon $ and $\delta $ 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.

### B. Effects of varying $\epsilon $ and *δ* on DHB and FSN points

*δ*

#### 1. Effect of $\epsilon $ on DHB

The effects of $\epsilon $ on the DHBs $ M S S H$ are summarized by the two-parameter bifurcation diagrams of (13) projected onto $( V 2, C 1)$-space (see Fig. 10). As $ C 1$ decreases (or equivalently, as $\epsilon $ decreases), the upper Hopf moves to larger $ V 2$ values and eventually vanishes for $\epsilon $ small enough [see Fig. 10(a)]. This explains why the upper singular orbit $ \Gamma ( \epsilon , 0 ) S S$ in Fig. 6(a) for $\epsilon \u22600$ switches to a continuum of spikes in Fig. 6(b) as $\epsilon \u21920$. On the other hand, since there is a CDH on the lower $ L s$, the lower Hopf will converge to that CDH as $\epsilon \u21920$ [see Fig. 10(b) and recall Remark II.3]. When $\epsilon $ increases, the lower Hopf and the lower fold of $ M S S$ 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 S S$ [also see Figs. 9(a) and 9(c), lower magenta circle].

#### 2. Effect of *δ* on FSN points

*δ*

There is no CDH or $ FSN 1$ on the upper fold, hence, we only examine the effect of $\delta $ on FSN singularities on the lower $ L s$. Figure 11(a) shows as $\delta \u21920$ (or equivalently, $ \varphi 2\u21920$), 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 Fig. 11(b) and recall the condition (29)].

### C. Why there are no SAOs near the lower CDH

Before examining the effects of varying $\epsilon $ and $\delta $ 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(\delta )$ close to an $ FSN 1$ and $O(\epsilon )$ 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 $\epsilon $ and $\delta $ 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 $\epsilon $ and $\delta $ 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 $\epsilon $ varies or as $\delta $ increases. While as $\delta $ decreases, the trajectory will follow $ M S S$ 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 $\delta $ is not guaranteed.

On the other hand, with default $\epsilon $ and $\delta $ values, the trajectory jumps away from the lower fold before reaching the Hopf bifurcation and hence there is no DHB-induced small oscillations. Increasing $\delta $ will make the DHB less relevant, whereas increasing $\epsilon $ will move the DHB further away from the lower $ L s$ and eventually vanish upon crossing the actual fold of $ M S S$ [see Fig. 10(b)]. Thus, we do not expect to detect MMOs with increased $\epsilon $ or $\delta $. As $\epsilon $ or $\delta $ decreases, the trajectory should pass closer to the lower DHB point. This is because decreasing $\epsilon $ moves the Hopf bifurcation closer to the CDH singularity [see Fig. 10(b)] and reducing $\delta $ pushes the trajectory to travel along $ M S S$ 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 S S$, i.e., close to a double zero eigenvalue at a BT bifurcation of the $( V 1, w 1, V 2)$ 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 $ M S S$ along which the Jacobian matrix $ J SL$ (24) of the slow layer problem (13) has complex eigenvalues. Figure 9(d) shows the real part of the first two eigenvalues $\lambda $ of (24) evaluated along $ M S S$, excluding the third eigenvalue given by $ f 2 V 2$. The eigenvalues are real when there are two branches of curves for $ Re\lambda $ and become complex when the curves coalesce and become a single branch. In panel (d), the eigenvalues on the stable lower branch of $ M S S$ are initially real and negative. That is, the trajectory approaches the attracting $ M S S$ 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 $\lambda $ 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 $\epsilon $ 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 $\epsilon $ and $\delta $ 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 $\epsilon $ by changing $ C 1$ and vary $\delta $ by changing $ \varphi 2$. Recall that increasing (respectively, decreasing) $ C 1$ or $\epsilon $ slows down (respectively, speeds up) the fast variable $ V 1$, whereas increasing (respectively, decreasing) $ \varphi 2$ or $\delta $ speeds up (respectively, slows down) the superslow variable $ w 2$. Other (slow) variables are not affected.

Figure 3(a) summarizes the effects of $( C 1, \varphi 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 $\delta $ is approximately greater than $O(\epsilon )$. Specifically, we observe the following:

For fixed $\delta =0.053$ at $ \varphi 2=0.001$, slowing down the fast variable $ V 1$ by increasing $ C 1$ from its default value $ C 1=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 $\epsilon $ increases.

For fixed $\delta $ at $ \varphi 2=0.001$, speeding up the fast variable $ V 1$ via decreasing $ C 1$ 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 $\epsilon =0.1$ at $ C 1=8$, slowing down the superslow variable $ w 2$ by decreasing $ \varphi 2$ from the default value $ \varphi 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 $ \varphi 2$ decreases. Additionally, the amplitudes of the small oscillations also display a similar non-monotonic behavior [see Fig. 13(a)].

For fixed $\epsilon $ at $ C 1=8$, speeding up the superslow variable $ w 2$ by increasing $ \varphi 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.

#### 1. MMOs are robust to increasing *C*_{1}

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 Fig. 10(a)]. This change allows the trajectory to travel a longer distance along the stable part of $ M S S$ during phase $\u25ef 3$ and generate more SAOs, as shown in Fig. 14. To simplify the presentation, we omit singular orbits and focus only on $ M S S$ and bifurcation diagrams that are essential for organizing SAOs in the full model, as demonstrated earlier.

Recall that the small oscillations occurring during phase $\u25ef 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 $ C 1$ causes the upper HB in $( V 2, V 1)$-space (Fig. 14, 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 Figs. 9(c) and 14]. As $ C 1$ is increased to $11$, the trajectory passes over the HB to $ M S S r$ 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 $ V 2$ at the green triangle.

#### 2. Decreasing *C*_{1} leads to three MMOs/non-MMOs transitions

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

The initial decrease of $ C 1$ [e.g., from $ C 1=8$ to $ C 1=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 $ C 1$. Interestingly, a further decrease of $ C 1$ to a range of $ C 1\u2208(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 $ C 1\u2208(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 $\u25ef 3$ grow in amplitude as the trajectory spirals away from $ M S S$. 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 Figs. 6(b) and 6(c). As a result, there is no more MMOs [see Fig. 15(d), $ C 1=1.4$].

To better understand why the SAOs for $ C 1\u2208(1.46,4.2)$ grow in amplitude [see Figs. 15(b) and 15(c)], we project the trajectory when $ C 1=2$ onto $( V 1, w 1, V 2)$ space (see Fig. 16). The blue triangle denotes a saddle-focus of the $( V 1, w 1, V 2)$ subsystem near the green square. During phase $\u25ef 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 $\u25ef 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}

#### 3. Decreasing *ϕ*_{2} preserves MMOs but causes non-monotonic effects on SAOs

*ϕ*

_{2}

Decreasing the perturbation parameter $ \varphi 2$ (i.e., reducing $\delta $) 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 $ M S S$ 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 $ \varphi 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 $ \varphi 2=0.0007$, the MMOs exhibit less SAOs with larger amplitudes than the SAOs at $ \varphi 2=0.0008$ [see Fig. 13(a), the third and the fourth row]. The number of SAOs increases and their amplitudes decrease again as $ \varphi 2$ decreases from $0.0007$ to $0.0006$. This alternating pattern of changes in the number and amplitude of SAOs repeats as $ \varphi 2$ continues to decrease [see Fig. 13(a)].

Our analysis reveals that as $ \varphi 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 $ \varphi 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 S S$ at the beginning of phase $\u25ef 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 S S a$ during phase $\u25ef 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 S S a$ 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 $ \varphi 2$ causes non-monotonic effects on SAOs.

#### 4. Increasing *ϕ*_{2} leads to five MMOs/non-MMOs transitions

*ϕ*

_{2}

When $ \varphi 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 $ \varphi 2$ large enough. Indeed, we observe a total of five transitions between MMOs and non-MMOs as $ \varphi 2$ increases, and eventually, MMOs are lost for $ \varphi 2>0.007$ [i.e., $\delta \u2265O( \epsilon 1 3)]$. The mechanism driving these MMOs/non-MMOs transitions over the increase of $ \varphi 2$ is similar to the mechanism underlying the non-monotonic effects on SAOs when $ \varphi 2$ is decreased. Specifically, if no additional big (full) spike before phase $\u25ef 3$ is lost with the increase of $ \varphi 2$, there will be fewer SAOs with larger amplitudes or no MMOs as one would naturally anticipate [e.g., when $ \varphi 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 $ \varphi 2$, changes to the SAOs will be reversed such that there will be more SAOs with smaller amplitude [e.g., when $ \varphi 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 $ \varphi 2\u22650.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 *g*_{syn} **=** 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 $\epsilon $ and $\delta $. Specifically, increasing $\epsilon $ manifests the DHB-like characteristics, while an increase in $\delta $ leads to dominance of the canard mechanism.

### A. Relation of the trajectory to *M*_{S} and *M*_{SS}

The solution of (2) for $ g syn=4.4$ in Fig. 1 is projected onto the $( V 1, w 1, V 2)$-space, together with the singular orbit from Fig. 6(d) (green curve), critical manifold $ M S$ (blue surface), fold $ L s$ (blue curve), and the superslow manifold $ M S S$ (red curves) (see Fig. 17). As the coupling strength $ g syn$ increases from $4.3$ to $4.4$, two new features regarding the upper $ M S S$ emerge. First, the stability of the upper $ M S S$ changes at a fold point $ L s s 1$ (the green triangle) rather than a DHB. Specifically, as $ V 2$ increases, the equilibrium along the upper $ M S S$ 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 S S$ 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(\delta )$ close to a folded saddle-node singularity FSN $ 1$ (upper blue star) and $O(\epsilon )$ close to a HB (upper red circle). This is further confirmed in Fig. 18, 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 limits $(\epsilon ,\delta )\u2192(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 $\u25ef 3$, the singular orbit (green) traces $ M S S a$ and jumps down at the fold point. The full trajectory $ \Gamma ( \epsilon , \delta )$ does not immediately jump at the fold. Instead, there is a delay before $ \Gamma ( \epsilon , \delta )$ jumps to the lower attracting manifold $ M S L$. Small amplitude oscillations are observed as the trajectory passes through the neighborhood of the canard and DHB points. For smaller $\epsilon $ or $\delta $ 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 $ M S S a$ during phase $\u25ef 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 Sec. III C.

At default parameters, $\delta =O(\epsilon )$. 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 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 $( V 1, V 2, w 2)$-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 $( V 1, V 2)$ 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 $ M S$, $ M S S$ and their relations to the full trajectory.

### B. MMOs are organized by both canard and DHB mechanisms

#### 1. Canard dynamics

Figure 19 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 $ FSN 1$, which occurs $O(\delta )$ 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 Fig. 20). As discussed in Sec. 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 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 $( 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.

#### 2. 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 $( V 2, V 1)$-plane. Starting at the green square, the trajectory exhibits SAOs as it follows the upper branch of $ M S S$ toward the left. As the trajectory passes through the attracting region of $ M S S$, the oscillation amplitude decays in a typical DHB fashion. Moreover, the orbit experiences a delay along the repelling $ M S S$ 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 S S$, 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 $ V 2$ switches from a superslow to a slow timescale at the green triangle, as discussed earlier when $ g syn=4.3$.

Thus, for $\delta = O(\epsilon )$, 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. First, we increased $\epsilon $ 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 ( $\epsilon =O(1)$), at which the folded singularities were no longer relevant (see Fig. 22). Second, we increased $\delta $ by raising $ \varphi 2$ to drive the system closer to (1F, 3S) splitting, and observed that SAOs persisted for $ \varphi 2$ as large as 0.01 ( $\delta =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 $\epsilon $ 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 Fig. 3(b). Specifically, MMOs persist over biologically relevant ranges of $ C 1\u2208(0.1,80)$ and $ \varphi 2\u2208( 10 \u2212 4,0.01)$ (also see Figs. 22 and 23).

The robustness of MMOs or SAOs to decreasing $\epsilon $ or $\delta $ 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 $\epsilon $ or $\delta $ on the features of small oscillations within MMOs, as decreasing them yields analogous effects but reversed. We find that MMOs with $\delta = O(\epsilon )$ exhibit both canard- and DHB-like features, whereas by tuning $\delta \u2265 O( \epsilon )$, DHB-like features diminish and the canard mechanism dominates.

Below we summarize the effects of increasing the singular perturbations:

For fixed $ \varphi 2=0.001$, increasing $\epsilon $ (i.e., increasing $ C 1$) enhances the DHB-like features of the SAOs (see Fig. 22).

For fixed $ C 1=8$, increasing $\delta $ (i.e., increasing $ \varphi 2$) makes the canard mechanism dominant and causes DHB-like features to gradually vanish (see Fig. 23).

#### 1. Increasing *C*_{1} 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 22 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 magenta 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 S S$. 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.

#### 2. Increasing *ϕ*_{2} makes the canard mechanism dominate

*ϕ*

_{2}

Increasing $ \varphi 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 S S$ and the DHB points become less relevant and eventually no longer meaningful for sufficiently large $ \varphi 2$ (e.g., $ \varphi 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 $ \varphi 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 $ \varphi 2$ increases, small oscillations tend to pull away from $ M S S$ and lose their DHB-like features (i.e., oscillations with decaying amplitude and a delay after passing the HB). Specifically, MMOs with $\delta =O(\epsilon )$ [e.g., Figs. 19 and 23(a)] exhibit both canard and DHB characteristics. When $\delta \u2248 1 2 \epsilon $ [Fig. 23(b)], there are still some DHB-like features. Further increasing $\delta $ to $\delta =0.006\u2248 \epsilon $ or $\delta =0.01\u2248 \epsilon 1 4$ [Figs. 23(c) and 23(d)], the trajectories no longer closely follow $ M S S$ and the amplitudes of small oscillations become almost constant, which reflects the absence of DHB-like features.

In summary, increasing $ \varphi 2$ makes the canard mechanism dominate and the SAOs exhibit fewer DHB-like features. Conversely, decreasing $ \varphi 2$ brings the solution and $ M S S$ closer together and amplifies the DHB characteristics of the sustained MMOs (see Fig. 23).

## 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.^{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 ( $ g syn=4.3$ and $4.4 mS / cm 2$). 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 $ g syn=4.3$, only the $\delta \u21920$ 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 $ g syn=4.4$, both $\epsilon \u21920$ and $\delta \u21920$ 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 subjects: the critical (or slow) manifold $ M S$ and the superslow manifold $ M S S$. The point where a fold of the critical manifold $ L s$ intersects the superslow manifold $ M S S$ 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 $\epsilon $ and $\delta $. 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 $ L s$ 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 $ L s$ did not support SAOs in either the case of $ g syn=4.3$ or $ g syn=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 $ M S S$. 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 $ g syn>4.3212$ (see Fig. 2), increasing $ g syn$ causes a shift of the CDH toward lower $ V 2$ 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 $ V 1$ spikes (data not shown). With higher $ g syn$ values, the transition from $ V 1$ spikes to $ V 1$ plateau or SAOs can happen during phase $\u25ef 1$ instead of phase $\u25ef 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 $\delta $ increases when $ g syn=4.3$ (see Sec. III). Ultimately, a sufficiently large value of $ g syn$ will bring the upper CDH below the minimal $ V 2$ of the trajectory, resulting in a complete loss of MMOs. A comprehensive investigation of these complex transitions induced by variations in $ g syn$ 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 $ g syn=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 $ V 1$ via decreasing $\epsilon $ 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 $\epsilon $ 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 $\epsilon $ 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 S S$ [see Figs. 15(b), 15(c), and 16]. Eventually, SAOs disappeared entirely when $ V 1$ became so rapid that it failed to remain in proximity of $ M S S$ to generate small oscillations.

. | g_{syn} = 4.3
. | g_{syn} = 4.4
. |
---|---|---|

Mechanisms | No upper CDH; Only DHB mechanism | Upper CDH exists; Canard and DHB mechanisms coexist |

Increasing C_{1} (or $\epsilon $) | MMOs are preserved with more DHB characteristics. | MMOs are preserved with more DHB characteristics. |

Decreasing C_{1} (or $\epsilon $) | Two MMOs/non-MMOs transitions are observed before a complete loss of MMOs. | MMOs are preserved and organized by both mechanisms. |

Increasing ϕ_{2} (or δ) | Four MMOs/non-MMOs transitions are observed before MMOs are entirely lost. | MMOs are preserved with more canard-like and less DHB features. |

Decreasing ϕ_{2} (or δ) | MMOs are preserved, but a non-monotonic effect on the small oscillations is observed. | MMOs are preserved with both DHB and canard characteristics. |

. | g_{syn} = 4.3
. | g_{syn} = 4.4
. |
---|---|---|

Mechanisms | No upper CDH; Only DHB mechanism | Upper CDH exists; Canard and DHB mechanisms coexist |

Increasing C_{1} (or $\epsilon $) | MMOs are preserved with more DHB characteristics. | MMOs are preserved with more DHB characteristics. |

Decreasing C_{1} (or $\epsilon $) | Two MMOs/non-MMOs transitions are observed before a complete loss of MMOs. | MMOs are preserved and organized by both mechanisms. |

Increasing ϕ_{2} (or δ) | Four MMOs/non-MMOs transitions are observed before MMOs are entirely lost. | MMOs are preserved with more canard-like and less DHB features. |

Decreasing ϕ_{2} (or δ) | MMOs are preserved, but a non-monotonic effect on the small oscillations is observed. | MMOs are preserved with both DHB and canard characteristics. |

As one would expect, MMOs when $ g syn=4.3$ is also sensitive to increasing $\delta $, which speeds up the superslow variable and thus makes the DHB mechanism less relevant. Interestingly, however, increasing $\delta $ 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 $\delta $, 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 $\delta $ 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 $\epsilon $ or decreasing $\delta $. 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 $\epsilon $ [see Figs. 12(a) and 14]. Nonetheless, decreasing $\delta $ 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 $\delta $ is reduced [see Fig. 13(a)]. Our analysis showed that the mechanism underlying such non-monotonic effects on SAOs over the decrease of $\delta $ is similar to the mechanism that drives multiple MMOs/non-MMOs transitions as $\delta $ 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 $\epsilon $ and decreasing $\delta $, 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$ [Fig. 10(a)], decreasing $\epsilon $ at $ g syn=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 $\epsilon \u21920$ singular limit. Hence, MMOs at $ g syn=4.4$ persist as $\epsilon $ decreases and are organized by both mechanisms. On the other hand, increasing $\delta $ 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 $\epsilon $ or decreasing $\delta $. Similarly, such MMOs should show stronger robustness to decreasing $\epsilon $ or increasing $\delta $. 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 $ \varphi 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.

## ACKNOWLEDGMENTS

This work was supported, in part, by NIH/NIDA R01DA057767 to Y.W., as part of the NSF/NIH/DOE/ANR/BMBF/BSF/NICT/AEI/ISCIII Collaborative Research in Computational Neuroscience Program.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Ngoc Anh Phan:** Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Yangyang Wang:** Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Software (equal); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

### APPENDIX A: DIMENSIONAL ANALYSIS

From system (A1), we can see that the voltage $ V 1$ evolves on a fast timescale $ ( C 1 g max = 1 ms )$, $( w 1, V 2)$ evolve on a slow timescale $ ( min ( \tau w ( V 1 ) ) \varphi 1 \u2248 C 2 g max = O ( 10 ) ms )$, whereas $ w 2$ evolves on a superslow timescale $ ( min ( \tau w ( V 2 ) ) \varphi 2 = O ( 100 ) ms )$.

### APPENDIX B: FOLDED SADDLE-NODE (FSN) SINGULARITIES OF (18)

### 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 Sec. II D and Fig. 4).

### APPENDIX D: DECREASING $ \varphi 2$ WHEN $ g syn=4.3$ CAUSES NON-MONOTONIC EFFECTS ON SAOS

In this subsection, we explain the mechanism underlying the non-monotonic effects of decreasing $ \varphi 2$ on SAOs [see Fig. 13(a) in Sec. III D]. We claim that (1) if an additional (full) big spike is generated during the process of decreasing $ \varphi 2$, the number of SAOs will decrease and their amplitudes will increase (e.g., when $ \varphi 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 $ \varphi 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 Fig. 24).

Figure 24(a) shows that the amplitude of SAOs with $ \varphi 2=0.0008$ (black) is smaller than SAOs with $ \varphi 2=0.0007$ (dark yellow). There are also more black SAOs when $ \varphi 2=0.0008$ [see Fig. 13(a)], 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 yellow one has a slower evolving rate of $ V 2$, which is slaved to $ w 2$, during phase $\u25ef 1$. Thus, one extra full spike occurs for the yellow trajectory compared to the black trajectory during the jump at phase $\u25ef 2$. As a result, the black trajectory, after its last big spike, approaches the black square at its maximum $ V 2$ along $ M S S$. In contrast, the yellow 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 yellow trajectory that is being pushed further away from $ M S S$ exhibits fewer small oscillations with larger amplitudes compared with the black trajectory, which remains close to $ M S S$.

When $ \varphi 2$ is reduced from $0.0007$ to $0.0006$, the solution trajectory changes to the green curve in Fig. 24(b). Slowing down of $ w 2$ even further keeps the number of full spikes between the yellow and green trajectories the same, but now all five full green spikes occur within phase $\u25ef 1$, similar to the case when $ \varphi 2=0.0008$. As a result, the green and black SAOs exhibit similar characteristics as shown in Fig. 24(b), resulting in more SAOs with smaller amplitudes than the yellow solution in Fig. 24(a).

### APPENDIX E: INCREASING $ \varphi 2$ WHEN $ g syn=4.3$ LEADS TO MULTIPLE MMOS/NON-MMOS TRANSITIONS

In this subsection, we explain the mechanism underlying the MMOs/non-MMOs transitions induced by an increase of $ \varphi 2$ in the case of $ g syn=4.3$ (see Fig. 13 in Sec. III D), which is similar to the mechanism that causes the non-monotonic effects on SAOs when $ \varphi 2$ is decreased (see Appendix D).

During the increase of $ \varphi 2$ which speeds up $ w 2$, the big spikes produced during phase $\u25ef 1$ will gradually decrease. If one big (full) spike before phase $\u25ef 3$ is lost during this process, there will be more SAOs with smaller amplitudes [e.g., when $ \varphi 2$ increases from $0.001$ to $0.0012$ in Fig. 25(a)] or the earlier lost MMOs will reappear [e.g., when $ \varphi 2$ increases from 0.0016 to 0.0018 in Fig. 25(c)]. As $ \varphi 2$ increases from $0.001$ to $0.0012$ in Fig. 25(a), the number of full spikes during phases $\u25ef 1$ and $\u25ef 2$ decreases by one. Consequently, all three green full spikes occur within phase $\u25ef 1$, while there is a large black spike during phase $\u25ef 2$. As explained in the previous subsection, this causes the green trajectory to approach the maximum $ V 2$ along $ M S S$ 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 $ \varphi 2$ from $0.0016$ to $0.0018$ leads to a decrease in the number of big spikes that occur before phase $\u25ef 3$, more SAOs are observed and the transition from non-MMOs to MMOs occurs [see Fig. 25(c)].

In contrast, increasing $ \varphi 2$ from $0.0012$ to $0.0016$ does not change the number of large spikes before phase $\u25ef 3$ [compare the green and dark yellow trajectories in Figs. 25(a) and 25(b)]. However, speeding up of $ \varphi 2$ pushes the third big yellow spike to occur during the jump, similar to the black trajectory. This causes the trajectory to be pushed away from $ M S S$, resulting in fewer SAOs with larger amplitudes. In fact, in this case, SAOs in the yellow trajectory are eliminated, and, therefore, a transition from MMOs to non-MMOs results.