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

Ngoc Anh Phan; Yangyang Wang

doi:10.1063/5.0181308

Chaos. 2024 May; 34(5): 053119.

Published online 2024 May 7. doi: 10.1063/5.0181308

PMCID: PMC11087137

PMID: 38717416

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

Ngoc Anh Phan^1,^a) and Yangyang Wang^2,^b)

Author information Article notes Copyright and License information PMC Disclaimer

Associated Data

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

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.

I. INTRODUCTION

Mixed-mode oscillations (MMOs) are frequently perceived in the dynamical systems involving multiple timescales;¹ 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);¹⁸ see Ref. ¹ 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.²⁸ 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.³⁵ 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. ¹⁰^,³⁰^,¹¹^,²⁷^,²⁵^,²⁶, and ^36–42 and references therein). Initial approaches were to consider three-dimensional systems,

\begin{matrix} ε \frac{d x}{d t} = f (x, y, z), \\ \frac{d y}{d t} = g (x, y, z), \\ \frac{d z}{d t} = δ h (x, y, z), \end{matrix}

(1)

with special cases $δ = ε$ or $δ = \sqrt{ε}$ .^{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. ²⁶ 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. ³⁵. The model equations are given by

\begin{matrix} \frac{d V_{1}}{d t} = (I_{1} - g_{Ca} m_{\infty} (V_{1}) (V_{1} - V_{Ca}) - g_{K} w_{1} (V_{1} - V_{K}) \\ - g_{L} (V_{1} - V_{L}) - g_{syn} S (V_{2}) (V_{1} - V_{syn})) / C_{1}, \\ \frac{d w_{1}}{d t} = ϕ_{1} (w_{\infty} (V_{1}) - w_{1}) / τ_{w} (V_{1}), \\ \frac{d V_{2}}{d t} = (I_{2} - g_{Ca} m_{\infty} (V_{2}) (V_{2} - V_{Ca}) \\ - g_{K} w_{2} (V_{2} - V_{K}) - g_{L} (V_{2} - V_{L})) / C_{2}, \\ \frac{d w_{2}}{d t} = ϕ_{2} (w_{\infty} (V_{2}) - w_{2}) / τ_{w} (V_{2}), \end{matrix}

(2)

with

\begin{matrix} S (V_{i}) = α (V_{i}) / (α (V_{i}) + β), \\ α (V_{i}) = 1 / (1 + \exp (- (V_{i} - θ_{s}) / σ_{s})), \\ m_{\infty} (V_{i}) = 0.5 (1 + \tanh ((V_{i} - K_{1}) / K_{2})), \\ w_{\infty} (V_{i}) = 0.5 (1 + \tanh ((V_{i} - K_{3}) / K_{4})), \\ τ_{w} (V_{i}) = 1 / \cosh ((V_{i} - K_{3}) / 2 K_{4}) . \end{matrix}

(3)

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

TABLE I.

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

Parameter values
C ₁	8 μF/cm²	I ₁	0 μA/cm²	ϕ ₁	0.01
C ₂	100 μF/cm²	I ₂	60 μA/cm²	ϕ ₂	0.001
V _Ca	120 mV	g _Ca	4 mS/cm²	K ₁	−1.2 mV
V _K	−84 mV	g _K	8 mS/cm²	K ₂	18 mV
V _L	−60 mV	g _L	2 mS/cm²	K ₃	12 mV
V _syn	30 mV	θ _s	−20 mV	K ₄	17.4 mV
β	0.5 ms⁻¹	σ _s	10 mV

Open in a separate window

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. ³⁵ 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. ³⁵ (e.g., $g_{syn} = 4.1$ and $5.1$ in Fig. 1).

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g001.jpg

Open in a separate window

FIG. 1.

Time traces of the model (2) for different values of $g_{syn}$ . MMOs are observed for $g_{syn} = 4.3$ and $g_{syn} = 4.4$ .

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. ²⁶ (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 \in {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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g002.jpg

Open in a separate window

FIG. 2.

Bifurcation curves of DHB (red) and CDH singularities (blue) for (2) with respect to $g_{syn}$ . Specifically, these curves represent the upper DHB and upper CDH, corresponding to larger $V_{1}$ and $V_{2}$ values. The lower CDH and DHB, associated with smaller $V_{i}$ values, are not presented here. The two vertical asymptotes are given by $g_{syn} \approx 4.2628$ and $g_{syn} \approx 4.3213$ .

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g004.jpg

Open in a separate window

FIG. 4.

Projections to $(V_{1}, V_{2}, w_{2})$ -space of the critical manifold fold surfaces $L_{s}$ (blue surface) for [(a) and (c)] $g_{s y n} = 4.3$ and [(b) and (d)] $g_{s y n} = 4.4$ . Also shown are the curves of folded singularities $M$ including folded node (solid green), folded saddle (dashed green), and two types of folded saddle-nodes $FSN$ (blue star: ${FSN}^{1}$ ; cyan star: ${FSN}^{2}$ ). The yellow curve consists of mostly folded foci points and small segments of other singularities (e.g., folded node, folded saddle) that are barely visible and hence are not displayed here. In the top two panels [(a) and (b)] when $δ \neq 0$ , an ${FSN}^{1}$ point (blue star) is $O (δ)$ close to a CDH singularity (blue diamond), whereas an ${FSN}^{2}$ (cyan star) is far away from any CDH. In the lower two panels [(c) and (d)] at the singular limit $δ = 0$ , the ${FSN}^{1}$ point becomes a CDH singularity (blue star overlapping with blue diamond). The center subspace of an ${FSN}^{1}$ (respectively, an ${FSN}^{2}$ ) is denoted by a pink plane (respectively, a yellow plane). It follows that the center manifolds of both ${FSN}^{1}$ and ${FSN}^{2}$ are transverse to $L_{s}$ .

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. ³⁵), 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. ³⁵. 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.

In this paper, we focus only on MMOs at $g_{syn} = 4.3$ and $g_{syn} = 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 problem:

\begin{matrix} ε \frac{d V_{1}}{d t_{s}} = f_{1} (V_{1}, w_{1}, V_{2}), \\ \frac{d w_{1}}{d t_{s}} = g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{s}} = f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{s}} = δ g_{2} (V_{2}, w_{2}), \end{matrix}

(4)

where $ε = 0.1$ , $δ = 0.053$ , $t_{s}$ is the slow dimensionless time variable, $f_{1}$ , $f_{2}$ , $g_{1}$ , and $g_{2}$ are $O (1)$ functions specified in (A1) in Appendix A, which include details of the nondimensionalization procedure. For simplicity, we did not rescale $V_{1}$ and $V_{2}$ in (4) as the scalings of voltage have no influence on the timescales.

We call system (4) that is described over the slow timescale the slow system in which $V_{1}$ evolves on a timescale of $O (ε^{- 1})$ , $(w_{1}, V_{2})$ on $O (1)$ and $w_{2}$ on $O (δ)$ . Introducing a superslow time $t_{s s} = δ t_{s}$ yields an equivalent description of dynamics,

\begin{matrix} ε δ \frac{d V_{1}}{d t_{s s}} = f_{1} (V_{1}, V_{2}, w_{1}), \\ δ \frac{d w_{1}}{d t_{s s}} = g_{1} (V_{1}, w_{1}), \\ δ \frac{d V_{2}}{d t_{s s}} = f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{s s}} = g_{2} (V_{2}, w_{2}), \end{matrix}

(5)

which evolves on the superslow timescale and is called the superslow system. Alternatively, defining a fast time $t_{f} = t_{s} / ε$ , we obtain the following fast system:

\begin{matrix} \frac{d V_{1}}{d t_{f}} = f_{1} (V_{1}, V_{2}, w_{1}), \\ \frac{d w_{1}}{d t_{f}} = ε g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{f}} = ε f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{f}} = ε δ g_{2} (V_{2}, w_{2}), \end{matrix}

(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 $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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g003.jpg

Open in a separate window

FIG. 3.

Regions of MMOs (yellow) and non-MMOs (blue) of the full system (2) in $(ϕ_{2}, C_{1})$ -space for (A) $g_{syn} = 4.3$ and (B) $g_{syn} = 4.4$ . Increasing $C_{1}$ slows down the fast variable $V_{1}$ , whereas increasing $ϕ_{2}$ speeds up the superslow variable $w_{2}$ . The timescales of $w_{1}$ and $V_{2}$ remain unaffected. The red star marks the default parameter values of $C_{1}$ and $ϕ_{2}$ as given in Table I. (a) $g_{syn} = 4.3$ . While MMOs are robust to increasing $C_{1}$ and decreasing $ϕ_{2}$ , decreasing $C_{1}$ or increasing $ϕ_{2}$ leads to multiple transitions between MMOs and non-MMOs (crossings between the dashed lines with the yellow/blue boundaries). (b) $g_{syn} = 4.4$ . MMOs are robust to changes of both $C_{1}$ and $ϕ_{2}$ over the ranges of $0.1 \leq C_{1} \leq 20$ and $0.1 e - 3 \leq ϕ_{2} \leq 8 e - 3$ . Note that the MMOs at $g_{syn} = 4.4$ will eventually vanish for $C_{1}$ and $ϕ_{2}$ large enough at which there is no more timescale separation (data not shown).

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 $δ \to 0$ singular limit provides a faithful prediction for the observed MMOs, whereas the $ε \to 0$ limit does not. This observation pinpoints the DHB from the $δ \to 0$ 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 $δ \leq O (ε)$ . Increasing $δ$ via increasing $ϕ_{2}$ or decreasing $ε$ via decreasing $C_{1}$ to a degree where $δ > O (ε)$ leads to multiple MMOs/non-MMOs transitions. MMOs are completely lost when $δ \geq O (ε^{\frac{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 $ε \to 0$ and $δ \to 0$ 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 $δ = O (ε)$ , MMOs exhibit both canard and DHB features. Upon tuning $δ \geq O (\sqrt{ε})$ , 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. ³⁹ 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 $(ε, δ) \to (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 $ε$ or $δ$ [i.e., varying $C_{1}$ and $ϕ_{2}$ in (2)]. We explain the transitions between MMO and non-MMO dynamics as we vary one perturbation parameter while keeping the other fixed at its default value, as illustrated by the two lines in 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 $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 $ε$ 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. ³⁵, 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. ³⁵. 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. $ε \to 0$ Singular limit

Fixing $δ > 0$ and taking $ε \to 0$ in the fast system (6) yields the one-dimensional (1D) fast layer problem, a system that describes the dynamics of the fast variable, $V_{1}$ , for fixed values of the other variables,

\frac{d V_{1}}{d t_{f}} = f_{1} (V_{1}, w_{1}, V_{2}) .

(7)

The set of equilibrium points of the fast layer problem is called the critical manifold and is denoted as $M_{S}$ ,

M_{S} := {(V_{1}, w_{1}, V_{2}, w_{2}) : f_{1} (V_{1}, w_{1}, V_{2}) = 0} .

(8)

Although $M_{S}$ is a three-dimensional (3D) manifold in $R^{4}$ space, it does not depend on $w_{2}$ . We can solve $f_{1} (V_{1}, w_{1}, V_{2}) = 0$ for $w_{1}$ as a function of $V_{1}$ and $V_{2}$ and can, therefore, represent $M_{S}$ as

w_{1} = F_{1} (V_{1}, V_{2})

(9)

for a function $F_{1}$ . It is well known that, for sufficiently small $ε > 0$ , normally hyperbolic parts of $M_{S}$ each perturb to a locally invariant manifold called a slow manifold, on which $w_{1}$ is given by an $O (ε)$ -perturbation of $F_{1}$ (Ref. ¹⁸); we simply use $M_{S}$ as a convenient numerical approximation of these slow manifolds.

$M_{S}$ is a 3D folded manifold with two-dimensional (2D) fold surface, $L_{s}$ , given by

L_{s} := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M_{S} : F_{1 V_{1}} = 0}

(10)

or equivalently

L_{s} := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M_{S} : f_{1 V_{1}} = 0},

(11)

where $F_{1 V_{1}}$ and $f_{1 V_{1}}$ denote the partial derivatives of $F_{1}$ and $f_{1}$ with respect to $V_{1}$ . The fold surface divides the critical manifold $M_{S}$ into attracting upper $(M_{S}^{U})$ and lower $(M_{S}^{L})$ branches where $F_{1 V_{1}} < 0$ and repelling middle branch $(M_{S}^{M})$ where $F_{1 V_{1}} > 0$ .

Taking the same limit, i.e., $ε \to 0$ with $δ > 0$ , in the slow system (4) yields the 3D slow reduced problem, a system that describes the dynamics of $w_{1}, V_{2}, w_{2}$ along $M_{S}$ ,

\begin{matrix} \frac{d w_{1}}{d t_{s}} = g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{s}} = f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{s}} = δ g_{2} (V_{2}, w_{2}), \end{matrix}

(12)

where $f_{1} = 0$ .

2. δ → 0 Singular limit

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

\begin{matrix} ε \frac{d V_{1}}{d t_{s}} = f_{1} (V_{1}, w_{1}, V_{2}), \\ \frac{d w_{1}}{d t_{s}} = g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{s}} = f_{2} (V_{2}, w_{2}), \end{matrix}

(13)

where the superslow variable $w_{2}$ is a parameter.

The set of equilibrium points of the slow layer problem (13) is defined to be the superslow manifold and is denoted as $M_{S S}$ ,

\begin{aligned} M_{S S} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \\ : f_{1} (V_{1}, w_{1}, V_{2}) = g_{1} (V_{1}, w_{1}) = f_{2} (V_{2}, w_{2}) = 0} . \end{aligned}

(14)

$M_{S S}$ is a 1D subset of $M_{S}$ . Similar to $M_{S}$ , the normally hyperbolic parts of $M_{S S}$ 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 $M_{S S}$ where Fenichel’s theory (GSPT) breaks down.

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

\frac{d w_{2}}{d t_{s s}} = g_{2} (V_{2}, w_{2}),

(15)

where $f_{1} = g_{1} = f_{2} = 0$ . The superslow motions of trajectories of (15) are slaved to $M_{S S}$ until nonhyperbolic points are reached.

3. $ε \to 0, δ \to 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 $δ \to 0$ in (12) or taking the limit $ε \to 0$ in (13) yields the same slow reduced layer problem

\begin{matrix} 0 = f_{1} (V_{1}, V_{2}, w_{1}), \\ \frac{d w_{1}}{d t_{s}} = g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{s}} = f_{2} (V_{2}, w_{2}), \end{matrix}

(16)

which describes the slow motion along $M_{S}$ and the superslow variable $w_{2}$ is fixed as a constant.

It follows that the double singular limits lead to three subsystems: the fast layer problem (7), the slow reduced layer problem (16) and the superslow reduced problem (15). In addition to the naturally expected fast/slow transitions and slow/superslow transitions, transitions directly from superslow to fast dynamics and from superslow to fast–slow relaxation oscillations have also been observed in Ref. ³⁵.

B. Slow reduced problem and canard dynamics

To investigate canard dynamics, we project the slow reduced problem (12) onto $(V_{1}, V_{2}, w_{2})$ to obtain a complete description of the dynamics along $M_{S}$ . To this end, we differentiate the graph representation of $M_{S}$ given by $w_{1} = F_{1} (V_{1}, V_{2})$ to obtain

\begin{matrix} F_{1 V_{1}} \frac{d V_{1}}{d t_{s}} = g_{1} (V_{1}, w_{1}) - F_{1 V_{2}} f_{2} (V_{2}, w_{2}), \\ \frac{d V_{2}}{d t_{s}} = f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{s}} = δ g_{2} (V_{2}, w_{2}), \end{matrix}

(17)

where $F_{1 V_{2}} := \frac{\partial F_{1}}{\partial V_{2}}$ . Note that the reduced system (17) is singular at the fold surfaces $L_{s}$ (10). Nonetheless, this singular term can be removed by a time rescaling $t_{s} = - F_{1 V_{1}} t_{d}$ , and we obtain the following desingularized system:

\begin{matrix} \frac{d V_{1}}{d t_{d}} = F_{1 V_{2}} f_{2} (V_{2}, w_{2}) - g_{1} (V_{1}, w_{1}) := F (V_{1}, w_{1}, V_{2}, w_{2}), \\ \frac{d V_{2}}{d t_{d}} = - F_{1 V_{1}} f_{2} (V_{2}, w_{2}) := G (V_{1}, V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{d}} = - δ F_{1 V_{1}} g_{2} (V_{2}, w_{2}) := H (V_{1}, V_{2}, w_{2}, δ) . \end{matrix}

(18)

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

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

\begin{aligned} E & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M_{S} : g_{1} (V_{1}, w_{1}) \\ = f_{2} (V_{2}, w_{2}) = g_{2} (V_{2}, w_{2}) = 0} . \end{aligned}

(19)

For the chosen parameter set in Table I, $E$ always lies on the repelling branch of $M_{S}$ and hence is unstable. In contrast to the ordinary singularities, the folded singularities are not equilibria of the full system. They lie on one-dimensional curves along the fold surface $L_{s}$ defined by

M := {(V_{1}, w_{1}, V_{2}, w_{2}) \in L_{s} : F (V_{1}, w_{1}, V_{2}, w_{2}) = 0} .

(20)

Folded singularities are special points that allow trajectories of (17) to cross the fold $L_{s}$ with nonzero speed. Such solutions are called singular canards.¹⁹ Note that when projecting to $(V_{1}, w_{1}, V_{2})$ -space, the condition $F = 0$ is redundant and the fold surfaces $L_{s}$ become curves that overlap with the folded singularity curves $M$ .

Since there is a curve of folded singularities, the Jacobian of (18) evaluated along $M$ (denoted as $J_{D}$ , see Appendix B) always has a zero eigenvalue and the eigenvector corresponding to this zero eigenvalue is tangent to the curve of folded singularities. Generically, the other two eigenvalues ( $λ_{w}, λ_{s}$ ) where $| λ_{w} | < | λ_{s} |$ have nonzero real part and are used to classify the folded singularities. Folded singularities with two real eigenvalues with the same sign (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 $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)

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

\begin{aligned} FSN & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : f_{2} Q (V_{1}, V_{2}, w_{2}) \\ = δ P (V_{1}, V_{2}, w_{2})}, \end{aligned}

(21)

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

\begin{matrix} {FSN}^{1} := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : f_{2} = δ K_{1} (V_{1}, V_{2}, w_{2})} \\ = {(V_{1}, w_{1}, V_{2}, w_{2}) \in L_{s} : \\ f_{2} = δ K_{1} (V_{1}, V_{2}, w_{2}), g_{1} = δ K_{2} (V_{1}, V_{2}, w_{2})}, \end{matrix}

(22)

or

\begin{aligned} {FSN}^{2} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : Q (V_{1}, V_{2}, w_{2}) \\ = δ K_{3} (V_{1}, V_{2}, w_{2})}, \end{aligned}

(23)

where $K_{1}, K_{2}$ , and $K_{3}$ are defined in Appendix B. We discuss below in Remark II.1 that both ${FSN}^{1}$ and ${FSN}^{2}$ are novel types of $FSN$ .²⁶

Remark II.1

In our parameter regime, the ordinary singularity point lies in the middle branch of critical manifold and is not involved in any bifurcations of the folded singularities. Hence, the $FSN$ singularities (21) 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. ²⁶, we identify ${FSN}^{1}$ and ${FSN}^{2}$ as novel types of $FSN$ .

Remark II.2

The condition (22) suggests that an ${FSN}^{1}$ is $O (δ)$ 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. ²⁶. 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 $δ \neq 0$ (top panels) and the singular limit $δ = 0$ (bottom panels). Recall the bifurcation of the upper CDH with respect to $g_{syn}$ is shown in 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 $δ \to 0$ singular limit, which exhibits delayed Hopf bifurcations that allow for interesting dynamics.

Let $J_{SL}$ denote the Jacobian matrix of (13) evaluated along the superslow manifold $M_{S S}$ , which is given by

J_{SL} = (\begin{matrix} \frac{1}{ε} f_{1 V_{1}} & \frac{1}{ε} f_{1 w_{1}} & \frac{1}{ε} f_{1 V_{2}} \\ g_{1 V_{1}} & g_{1 w_{1}} & 0 \\ 0 & 0 & f_{2 V_{2}} \end{matrix}),

(24)

where the nonzero entries denote partial derivatives.

The eigenvalues of $J_{SL}$ are given by $f_{2 V_{2}}$ and the eigenvalues of

J = (\begin{matrix} \frac{1}{ε} f_{1 V_{1}} & \frac{1}{ε} f_{1 w_{1}} \\ g_{1 V_{1}} & g_{1 w_{1}} \end{matrix}) .

Thus, the Hopf bifurcation points on $M_{S S}$ are given by $tr (J) = \frac{1}{ε} f_{1 V_{1}} + g_{1 w_{1}} = 0$ and $det J > 0$ . The former defining condition can be rewritten as

M_{S S}^{H} := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M_{S S} : f_{1 V_{1}} = - ε g_{1 w_{1}}} .

(25)

Remark II.3

It follows from (25) that an $M_{S S}^{H}$ is $O (ε)$ 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).

The isolated fold bifurcation points on $M_{S S}$ are located by letting $det J_{SL} = 0$ . That is,

\begin{aligned} L_{s s} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M_{S S} : f_{2 V_{2}} \\ = 0 or f_{1 V_{1}} g_{1 w_{1}} - f_{1 w_{1}} g_{1 V_{1}} = 0} . \end{aligned}

(26)

The fold points $L_{s s}$ that satisfy the former condition (denoted as $L_{s s}^{1}$ ) are the folds of the $V_{2}$ -nullcline (see Fig. 5, green circle and green triangle), which correspond to the transition between superslow dynamics along $M_{S S}$ and slow jumps. $L_{s s}$ given by the latter condition (denoted as $L_{s s}^{2}$ ) corresponds to the actual fold of $M_{S S}$ when projected to $(V_{1}, w_{1}, V_{2})$ -space. Since $g_{1 w_{1}} < 0$ , $f_{1 w_{1}} < 0$ and $g_{1 V_{1}} > 0$ , it follows that $L_{s s}^{2}$ lies on the middle branch of $M_{S}$ ( $f_{1 V_{1}} = f_{1 w_{1}} g_{1 V_{1}} / g_{1 w_{1}} > 0$ ) and hence will not play a role in dynamics. At the double singular limit $(ε, δ) \to (0, 0)$ , the $L_{s s}^{2}$ fold point of $M_{S S}$ 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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g005.jpg

Open in a separate window

FIG. 5.

Projection of the singular orbit (green) and the solution trajectory $Γ_{(ε, δ)}$ (black) of the full system (2) onto the phase plane of ( $V_{2}, w_{2}$ ) system with parameters given in Table I. The red curve is the $V_{2}$ -nullcline, which is the projection of the superslow manifold $M_{S S}$ . Yellow and green symbols mark the key transitions between the slow and superslow sections of the singular orbit and the perturbed oscillation, respectively: the star and square indicate the transitions from the slow to the superslow motions, and the circle and triangle mark the transitions from superslow to slow sections at the fold of the $V_{2}$ -nullcline. The circled numbers indicate four phases of the oscillations: superslow excursions along $M_{S S}$ during $◯ 1$ and $◯ 3$ and slow jumps at the fold of $M_{S S}$ during $◯ 2$ and $◯ 4$ .

D. Interaction between canard and delayed Hopf

To investigate the interaction between the canard and the delayed Hopf mechanisms in the double limit case ( $ε \to 0, δ \to 0$ ), we need to examine the slow reduced layer problem (16). The corresponding desingularized system is given by

\begin{matrix} \frac{d V_{1}}{d t_{d}} = F (V_{1}, w_{1}, V_{2}, w_{2}), \\ \frac{d V_{2}}{d t_{d}} = G (V_{1}, V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{d}} = 0, \end{matrix}

(27)

where $F$ and $G$ are defined in (18). Note that (27) is the $δ \to 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 $M_{S S}$ . The FSN condition at the double singular limit can be obtained from letting $δ \to 0$ in the ${FSN}^{1}$ condition (22) or the ${FSN}^{2}$ condition (23). This implies that, at the double singular limit, an FSN¹ becomes a CDH singularity [see Figs. 4(c) and 4(d)],

\begin{aligned} {FSN}_{(ε, δ) \to (0, 0)}^{1} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : f_{2} = 0} \\ = {(V_{1}, w_{1}, V_{2}, w_{2}) \in L_{s} : f_{2} = g_{1} = 0} \\ = M_{S S} \cap L_{s}, \end{aligned}

(28)

whereas an ${FSN}^{2}$ singularity is characterized by

\begin{aligned} {FSN}_{(ε, δ) \to (0, 0)}^{2} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : \\ Q (V_{1}, V_{2}, w_{2}) = 0} . \end{aligned}

(29)

According to Remarks II.2 and II.3, an ${FSN}^{1}$ singularity from the $ε$ viewpoint converges to a CDH as $δ \to 0$ and a DHB $M_{S S}^{H}$ from the $δ$ viewpoint converges to a CDH as $ε \to 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. ²⁶. 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. ²⁶.

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. ²⁶, 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,¹⁸ 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 $ε \to 0$ singular limit. Let $Γ_{(ε, 0)}^{S}$ and $Γ_{(ε, 0)}^{S S}$ denote trajectories of the slow layer problem (13) and the superslow reduced problem (15) obtained from the $δ \to 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 \in {F, S, S S}$ 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 $ε$ . The singular orbit (Fig. 5, green trajectory) consists of a superslow excursion along the left branch of $V_{2}$ -nullcline ( $Γ_{(ε, 0)}^{S S}$ , phase $◯ 1$ ), a slow jump at the lower fold of $V_{2}$ -nullcline up to its right branch ( $Γ_{(ε, 0)}^{S}$ , phase $◯ 2$ ), a superslow excursion through the active phase ( $Γ_{(ε, 0)}^{S S}$ , 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 $(V_{2}, w_{2})$ -space, the $V_{2}$ -nullcline (red) corresponds to the superslow manifold $M_{S S}$ and the singular orbit $Γ_{(ε, 0)}^{S} \cup Γ_{(ε, 0)}^{S S}$ overlaps with the singular orbit from the double singular limits due to its independence on $ε$ . For sufficiently small perturbation $δ$ , $Γ_{(ε, 0)}^{S} \cup Γ_{(ε, 0)}^{S S}$ perturbs to the full system trajectory $Γ_{(ε, δ)}$ shown by the black curve.

Figure 6(a) illustrates the singular orbit $Γ_{(ε, 0)}^{S} \cup Γ_{(ε, 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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g006.jpg

Open in a separate window

FIG. 6.

Projection of singular periodic orbit (green curve) of system (2) for $g_{syn} = 4.3$ [(a)–(c)] and $g_{syn} = 4.4$ (d) to $(V_{1}, w_{1}, V_{2})$ -space. Specifically, panels (a) and (c) show the singular orbits for $(ε, δ) = (0.1, 0)$ , and $(ε, δ) = (0, 0.053)$ , respectively. Panels (b) and (d) show the singular orbits at the double singular limit $(ε, δ) = (0, 0)$ . Also shown are the superslow manifold $M_{S S}$ (red curves), the critical manifold $M_{S}$ (blue surface), and folds of $M_{S}$ (blue curves). The solid (respectively, dashed) red curves are the branches of $M_{S S}$ that are attracting (respectively, repelling) under the superslow flow. Yellow symbols represent transition points between slow and superslow flow for the $(V_{2}, w_{2})$ oscillation as shown in Fig. 5. The red circle denotes the fast subsystem DHB $M_{S S}^{H}$ . ${FSN}^{1}$ and CDH are denoted by the blue star and diamond, respectively.

Starting from the yellow star in panel (a), the singular orbit is in phase $◯ 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 $◯ 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 $◯ 2$ ), a few more spikes occur before the slow flow $Γ_{(ε, 0)}^{S}$ travels to the yellow square on $M_{S S}$ , at which phase $◯ 3$ begins. After that, the superslow reduced problem takes over until the singular orbit $Γ_{(ε, 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 $◯ 1$ . As phase $◯ 4$ begins at the yellow triangle, several additional $V_{1}$ 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 $M_{S S}$ . Instead, we observe a continuum of $V_{1}$ spikes throughout phase $◯ 3$ . This is because the upper $M_{S S}$ becomes unstable for all $V_{2}$ values as $ε \to 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 $V_{1}$ jumps from $L_{s}$ , $Γ_{(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 $Γ_{(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 $◯ 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 $g_{syn} = 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 $M_{S S}$ differs between the two panels, can be viewed as an $O (ε)$ perturbation of the singular orbit in panel (b). Hence, when discussing the fast–slow oscillations in the full system, we still refer to different segments of them as being governed by (7), which describes fast $V_{1}$ jumps and (16), which describes the slow motion along $M_{S}$ .

For completeness, we also plot the singular orbit for $ε = 0$ but $δ \neq 0$ in Fig. 6(c). Instead of a continuum of $V_{1}$ spikes, we observe a finite number of $V_{1}$ spikes during phases $◯ 1$ and $◯ 3$ since $δ \neq 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)}^{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, $Γ_{(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 $◯ 4$ due to its unique starting position. Nonetheless, the construction of singular orbits during phase $◯ 2$ remains non-unique as discussed in the $g_{syn} = 4.3$ case. Later in 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²

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 $ε$ 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.

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¹ as $δ \to 0$ as discussed in Sec. II D. Moreover, the nearby HB bifurcation (red circle in the lower inset) is $O (ε)$ close to this CDH.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g007.jpg

Open in a separate window

FIG. 7.

Projection of an attracting MMO solution trajectory (black curve) of system (2) for $g_{syn} = 4.3$ from Fig. 1 to $(V_{1}, w_{1}, V_{2})$ -space. Also shown are parts of the singular orbit from Fig. 6(a) (green curves), the critical manifold $M_{S}$ (blue surface), folds of $M_{S}$ (blue curves), and the superslow manifold $M_{S S}$ (red curves). The upper inset shows the upper branch of $M_{S S}$ is always below the upper fold of $M_{S}$ , indicating the absence of an upper CDH. The black circle near the upper fold is the true equilibrium of the full system (2), which is unstable. The lower inset shows a magnified view around the lower CDH (blue diamond) at which the lower branch of $M_{S S}$ intersects the lower fold of $M_{S}$ . The lower HB bifurcation (red circle) is $O (ε)$ close to the CDH singularity. Other color coding and symbols have the same meaning as in Figs. 5 and and66.

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)}^{S S}$ perturbs to $Γ_{(ε, δ)}$ , 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 $◯ 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 $◯ 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 $Γ_{(ε, 0)}^{S S}$ 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 $(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).

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g008.jpg

Open in a separate window

FIG. 8.

The solution from Fig. 7 projected onto ( $V_{1}, V_{2}, w_{2}$ )-space. Also shown are the projections of the fold surface $L_{s}$ (blue surface), superslow manifold $M_{S S}$ (red curve), and folded singularities $M$ (green and yellow curves). The cyan star denotes ${FSN}^{2}$ defined in (23). Color codings of the folded singularity curves are the same as in Fig. 4. Other symbols have the same meanings as in Fig. 7.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g009.jpg

Open in a separate window

FIG. 9.

Solutions of (2) for $g_{syn} = 4.3$ , $C_{1} = 8$ and various values of $ϕ_{2}$ and bifurcation diagrams of the slow layer problem (13), projected onto $(V_{2}, V_{1})$ -plane. From (a) to (b) to (c), $ϕ_{2}$ increases from $0.0008$ to $0.0009$ to its default value of $0.001$ . In addition to the trajectory (black), the singular orbit (green) and $M_{S S}$ (red curve), also shown is the periodic orbit (PO) branches (solid cyan: stable; dashed cyan: unstable) born at the upper HB bifurcation. The upper two magenta circles denoting the saddle-node bifurcations of $M_{S S}$ exhibit the same $V_{2}$ values as the folds of the $V_{2}$ -nullcline (green circle and green triangle). The lower magenta circle (see the inset) represents the actual fold of $M_{S S}$ (denoted as $L_{s s}^{2}$ ) and is not a fold of the $V_{2}$ -nullcline. Other symbols have the same meanings as in Fig. 7. (d) Real part of the eigenvalues of the upper triangular block of $J_{SL}$ (24), the Jacobian matrix of the slow layer problem (13) along the superslow manifold $M_{S S}$ . The eigenvalues along $M_{S S}$ are real when there are two branches of $Re (λ)$ (blue curves) and complex when there is a single branch of $Re (λ)$ (black curve).

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 $δ$ 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 $δ \to 0$ singular orbits from Fig. 6(a) perturb to MMOs for various perturbations $δ$ . With small perturbations [Fig. 9(a), $δ \approx 0.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 $Γ_{(ε, δ)}$ under the default perturbation $δ \approx 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 $δ \approx 0.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 $δ \approx 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. ³⁵, 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 $δ$ on SAOs is non-monotonic. In certain cases, smaller perturbations can lead to fewer SAOs with even larger amplitudes. This effect is related to how the slow flow during phase $◯ 2$ approaches $M_{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 $◯ 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}$ .

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g013.jpg

Open in a separate window

FIG. 13.

Effect of varying $ϕ_{2}$ (or $δ$ ) on time traces of solutions of (2) for $g_{syn} = 4.3, C_{1} = 8$ and other parameters given in Table I. (a) Decreasing $ϕ_{2}$ preserves MMOs, but the amplitude and the number of the small oscillations alternates between increasing and decreasing with $ϕ_{2}$ . (b) Increasing $ϕ_{2}$ leads to a total of five transitions between MMOs and non-MMOs. MMOs exist for $ϕ_{2} < 0.0016$ [i.e., $δ \leq O (ε)$ ] and are completely lost for $ϕ_{2} > 0.007$ [i.e., $δ \geq O (ε^{\frac{1}{3}})]$ .

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.

B. Effects of varying $ε$ and δ on DHB and FSN points

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

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g010.jpg

Open in a separate window

FIG. 10.

2-parameter bifurcation curves for (13) projected onto ( $V_{2}, C_{1}$ )-space when $g_{syn} = 4.3$ . (a) The bifurcation curve of the upper HB (Fig. 7, upper red circle), with a horizontal asymptote at $C_{1} = 2.9097$ . (b) The bifurcation curves of the lower HB (red curve) and the lower fold $L_{s s}^{2}$ of $M_{S S}$ (magenta curve), which meet at a Bogdanov–Takens (BT) bifurcation. As $C_{1} \to 0$ (i.e., $ε \to 0$ ), the lower HB converges to the lower CDH (blue).

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g011.jpg

Open in a separate window

FIG. 11.

Relation of ${FSN}^{1}$ (22), ${FSN}^{2}$ (23), and the CDH point on the lower fold for $g_{syn} = 4.3$ and various values of $ϕ_{2}$ . (a) ${FSN}^{1}$ (blue star) converges to the $CDH$ (blue diamond) as $ϕ_{2} \to 0$ (equivalently, $δ \to 0$ ). (b) ${FSN}^{2}$ (cyan star) are located far away from the $CDH$ and converge to the leftmost cyan star at the intersection of the folded singularity curve $M$ and $Q (V_{1}, V_{2}, w_{2}) = 0$ [see (29)]. The yellow curve denotes the lower folded singularity curve without showing stability and types.

1. Effect of $ε$ on DHB

The effects of $ε$ 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 $ε$ decreases), the upper Hopf moves to larger $V_{2}$ values and eventually vanishes for $ε$ small enough [see Fig. 10(a)]. This explains why the upper singular orbit $Γ_{(ε, 0)}^{S S}$ in Fig. 6(a) for $ε \neq 0$ switches to a continuum of spikes in Fig. 6(b) as $ε \to 0$ . On the other hand, since there is a CDH on the lower $L_{s}$ , the lower Hopf will converge to that CDH as $ε \to 0$ [see Fig. 10(b) and recall Remark II.3]. When $ε$ 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 $δ$ on FSN singularities on the lower $L_{s}$ . Figure 11(a) shows as $δ \to 0$ (or equivalently, $ϕ_{2} \to 0$ ), the ${FSN}^{1}$ singularity converges to the lower CDH as demonstrated in our analysis [see (28)]. On the other hand, the ${FSN}^{2}$ singularity is significantly distant from the CDH [see Fig. 11(b) and recall the condition (29)].

C. 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 ${FSN}^{1}$ and $O (ε)$ close to a DHB. One would naturally expect to observe SAOs arising from Canard and/or DHB mechanisms near this lower CDH. Nonetheless, there is no SAOs near the lower fold regardless of $ε$ and $δ$ values. Below we explain why none of the two mechanisms produces SAOs.

As discussed earlier, the reason why there is no 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 $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 $δ$ is not guaranteed.

On the other hand, with default $ε$ and $δ$ values, the trajectory jumps away from the lower fold before reaching the Hopf bifurcation and hence there is no DHB-induced small oscillations. Increasing $δ$ will make the DHB less relevant, whereas increasing $ε$ will move the DHB further away from the lower $L_{s}$ and eventually vanish upon crossing the actual fold of $M_{S S}$ [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 $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 $λ$ 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 λ$ 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 $λ$ becomes real again shortly after. As a result, the trajectory has insufficient time to oscillate and we do not observe any small-amplitude oscillations before the trajectory jumps up to the outer periodic orbit branch.

D. Effects of varying $ε$ and δ on MMOs

Recall when $g_{syn} = 4.3$ , the only mechanism available for MMOs is the DHB mechanism. While there exist folded node singularities, they do not play any significant role in generating MMOs. In this subsection, we explore the effects of $ε$ and $δ$ on the dynamics of the full system (4) by mainly examining their effects on the upper DHB around which SAOs are observed. As before, we will vary $ε$ by changing $C_{1}$ and vary $δ$ by changing $ϕ_{2}$ . Recall that increasing (respectively, decreasing) $C_{1}$ or $ε$ slows down (respectively, speeds up) the fast variable $V_{1}$ , whereas increasing (respectively, decreasing) $ϕ_{2}$ or $δ$ 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}, ϕ_{2})$ on MMOs when $g_{syn} = 4.3$ . Our findings suggest that MMOs with only the DHB mechanism are robust to changes that slow down either the fast variable or the superslow variable, but they are vulnerable to perturbations that speed up either timescale to a degree where $δ$ is approximately greater than $O (ε)$ . Specifically, we observe the following:

$∙$
For fixed $δ = 0.053$ at $ϕ_{2} = 0.001$ , slowing down the fast variable $V_{1}$ by increasing $C_{1}$ from its default value $C_{1} = 8$ [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.
Open in a separate window
FIG. 12.
Effect of varying $C_{1}$ (or $ε$ ) on time traces of solutions of (2) for $g_{syn} = 4.3, ϕ_{2} = 0.001$ and other parameters given in Table I. (a) Increasing $C_{1}$ from its default value $8$ (equivalently, increasing $ε$ from $0.1$ ) preserves MMOs. The number of small oscillations increases with $C_{1}$ . (b) Decreasing $C_{1}$ from the default value leads to a transition from MMOs ( $C_{1} = 8$ ) to non-MMOs ( $C_{1} = 7$ ) to MMOs again ( $C_{1} = 3$ ) to non-MMOs ( $C_{1} = 1.4$ ).
$∙$
For fixed $δ$ at $ϕ_{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 $ε = 0.1$ at $C_{1} = 8$ , slowing down the superslow variable $w_{2}$ 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 $C_{1} = 8$ , speeding up the superslow variable $w_{2}$ by increasing $ϕ_{2}$ leads to a total of five transitions between MMOs and non-MMOs [see 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₁

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 $◯ 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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g014.jpg

Open in a separate window

FIG. 14.

Solutions of (2) for $g_{syn} = 4.3$ and various values of $C_{1}$ and bifurcation diagrams of the slow layer problem (13), projected to $(V_{2}, V_{1})$ -space. From left to right, $C_{1} = 9$ ( $ε = 0.1125$ ), $C_{1} = 10$ ( $ε = 0.125$ ), $C_{1} = 11$ ( $ε = 0.1375$ ). Increasing $C_{1}$ moves the upper HB bifurcation (red circle) to lower $V_{2}$ values and hence increases the number of small oscillations surrounding $M_{s s}$ . Color codings and symbols are the same as in Figs. 9(a)–9(c).

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 $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 and14].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₁ 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} \in (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} \in (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} \in (1.46, 4.2)$ and $V_{2}$ near the green square, (13) exhibits either unstable periodic orbits with negligible amplitudes or saddle-foci equilibrium characterized by a negative real eigenvalue and complex eigenvalues with positive real parts. As a result, the SAOs during phase $◯ 3$ grow in amplitude as the trajectory spirals away from $M_{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$ ].

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g015.jpg

Open in a separate window

FIG. 15.

Projections of the full system solutions and bifurcation diagrams of (13) for $g_{syn} = 4.3$ and (a) $C_{1} = 7$ ( $ε = 0.0875$ ), (b) $C_{1} = 3$ ( $ε = 0.0375$ ), (c) $C_{1} = 2$ ( $ε = 0.025$ ), (d) $C_{1} = 1.4$ ( $ε = 0.0175$ ). The upper DHB (red circle) moves to larger $V_{2}$ values with decreased $C_{1}$ and eventually vanishes for $C_{1} < 2.9$ . Color codings and symbols have the same meaning as in Figs. 9(a)–9(c).

To better understand why the SAOs for $C_{1} \in (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 $◯ 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⁴⁸ and a singular Hopf bifurcation in two-timescale settings.¹

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g016.jpg

Open in a separate window

FIG. 16.

The solution (black curve) from Fig. 15(c) projected onto $(V_{1}, w_{1}, V_{2})$ space, together with $M_{S S}$ (red curve). The blue triangle denotes a saddle-focus equilibrium on $M_{S S}$ with maximum $V_{2}$ , which has one negative real eigenvalue and a pair of complex eigenvalues with positive real parts ( $0.0027 \pm 0.31 i$ ). A stable manifold branch associated with the saddle-focus at the blue triangle is denoted by the magenta curve.

3. Decreasing ϕ₂ 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 $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 $ϕ_{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 $M_{S S}$ at the beginning of phase $◯ 3$ , leading to fewer SAOs with larger amplitudes. Hence, the amplitude and number of SAOs in the full system are not only determined by the size of the perturbation but also by how the flow approaches $M_{S S}^{a}$ during phase $◯ 2$ . As previously discussed, infinitely many singular orbit segments can be constructed during this phase. This leads to different ways for the full trajectory to reach $M_{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 $ϕ_{2}$ causes non-monotonic effects on SAOs.

4. Increasing ϕ₂ leads to five MMOs/non-MMOs transitions

When $ϕ_{2}$ increases, the superslow variable $w_{2}$ speeds up, moving the system closer to the 1-fast/3-slow splitting and making the DHB mechanism less relevant. Since there is no canard mechanism, MMOs should be eliminated for $ϕ_{2}$ large enough. Indeed, we observe a total of five transitions between MMOs and non-MMOs as $ϕ_{2}$ increases, and eventually, MMOs are lost for $ϕ_{2} > 0.007$ [i.e., $δ \geq O (ε^{\frac{1}{3}})]$ . The mechanism driving these MMOs/non-MMOs transitions over the increase of $ϕ_{2}$ is similar to the mechanism underlying the non-monotonic effects on SAOs when $ϕ_{2}$ is decreased. Specifically, if no additional big (full) spike before phase $◯ 3$ is lost with the increase of $ϕ_{2}$ , there will be fewer SAOs with larger amplitudes or no MMOs as one would naturally anticipate [e.g., when $ϕ_{2}$ increases from $0.0012$ to $0.0016$ , see 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} \geq 0.008$ for which the HB is no longer relevant.

For a more detailed discussion, we refer the readers to Appendix E.

IV. ANALYSIS OF MMOS WHEN g_syn = 4.4 mS/cm²

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.

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 (δ)$ close to a folded saddle-node singularity FSN¹ (upper blue star) and $O (ε)$ 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 $(ε, δ) \to (0, 0)$ .

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g017.jpg

Open in a separate window

FIG. 17.

Projection of a solution trajectory (black curve) of system (2) for $g_{syn} = 4.4$ from Fig. 1(b) to $(V_{1}, w_{1}, V_{2})$ -space. Also shown are parts of the singular orbit $Γ_{(0, 0)}^{F} \cup Γ_{(0, 0)}^{S} \cup Γ_{(0, 0)}^{S S}$ from Fig. 6(d) (green curve), the critical manifold $M_{S}$ (blue surface), folds $L_{s}$ of $M_{S}$ (blue curves), and the superslow manifold $M_{S S}$ (red curves). The solid (respectively, dashed) red curves represent the attracting (respectively, repelling) branches of $M_{S S}$ . Yellow and green symbols mark the transitions between slow and superslow pieces of the $(V_{2}, w_{2})$ oscillations as in Fig. 5. The blue diamonds denote CDH singularities—the intersection of $L_{s}$ and $M_{S S}$ . The upper CDH is a folded node and the lower CDH (almost overlapping with the lower DHB denoted by a red circle) is a folded focus. The upper CDH point is $O (δ)$ close to the folded saddle-node singularity ${FSN}^{1}$ (blue star) and $O (ε)$ close to the upper DHB (red circle). The black circle denotes the isolated ordinary singularity of the full system, whose type is a saddle-focus.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g018.jpg

Open in a separate window

FIG. 18.

The upper $CDH$ represents the interaction of canard and DHB mechanisms when $g_{syn} = 4.4$ . (Left) The upper ${FSN}^{1}$ singularity (blue star) converges to the upper $CDH$ (blue diamond) as $ϕ_{2} \to 0$ (or equivalently $δ \to 0$ ). (Right) As $C_{1}$ (or $ε$ ) decreases, the upper HB (red curve) moves to larger $V_{2}$ values and approaches the upper $CDH$ (blue curve) in the singular limit $ε \to 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 $M_{S S}^{a}$ and jumps down at the fold point. The full trajectory $Γ_{(ε, δ)}$ does not immediately jump at the fold. Instead, there is a delay before $Γ_{(ε, δ)}$ jumps to the lower attracting manifold $M_{S}^{L}$ . 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 and23),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 $◯ 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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g022.jpg

Open in a separate window

FIG. 22.

Time traces of the solutions of (2) (left panels) and their projections onto the corresponding bifurcation diagrams in $(V_{2}, V_{1})$ -space (right panels) for $g_{syn} = 4.4$ and (a) $C_{1} = 0.8$ $(ε = 0.01)$ , (b) $C_{1} = 40$ $(ε = 0.5)$ , and (c) $C_{1} = 80$ ( $ε = 1$ ). Increasing $C_{1}$ moves the upper DHB (red circle) to the lower $V_{2}$ and manifests DHB-like features. Color coding and symbols of the bifurcation diagrams on the right panel are the same as in Fig. 17.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g023.jpg

Open in a separate window

FIG. 23.

Time traces of the solutions of (2) (left panels) and their projections (black curves) onto $(w_{2}, V_{2}, V_{1})$ -space (right panels) for $g_{syn} = 4.4$ and (a) $ϕ_{2} = 0.0005$ (i.e., $δ = 0.0265 = \frac{1}{3} ε$ ), (b) $ϕ_{2} = 0.003$ (i.e., $δ = 0.159 = \frac{1}{2} \sqrt{ε}$ ), (c) $ϕ_{2} = 0.006$ (i.e., $δ = 0.318 = \sqrt{ε}$ ), and (d) $ϕ_{2} = 0.01$ (i.e., $δ = 0.56 = ε^{\frac{1}{4}}$ ). Increasing $ϕ_{2}$ preserves MMOs with more canard-like feature. Codings and symbols on the right panels are the same as in Fig. 20.

At default parameters, $δ = O (ε)$ . The emergence of SAOs for $g_{syn} = 4.4$ is governed by both the canard dynamics due to folded node singularities and the slow passage effects associated with the DHB on the upper fold $L_{s}$ , as seen in examples considered in Refs. ²⁶ and ²⁷. 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.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g019.jpg

Open in a separate window

FIG. 19.

Projection of the solution from Fig. 17 onto ( $w_{2}, V_{2}, V_{1}$ )-space. Also shown are the projections of the fold surface $L_{s}$ (blue surfaces), superslow manifold $M_{S S}$ , and folded singularities $M$ (curves of green and yellow). Other color coding and symbols have the same meaning as in Figs. 4 and and1717.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g020.jpg

Open in a separate window

FIG. 20.

Projections of the solution from Fig. 17 (black) and the singular 3D funnel volume corresponding to the curve of folded nodes onto $(w_{2}, V_{2}, V_{1})$ - space (top) and $(V_{2}, V_{1})$ - plane (bottom). The 3D funnel volume is bounded between the singular canard surface (magenta surface) and the fold surface $L_{s}$ . The black, cyan, and yellow curves denote solutions of (2) with different initial conditions. Other color coding and symbols have the same meaning as in Fig. 19.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g021.jpg

Open in a separate window

FIG. 21.

Projection of the solution trajectory and superslow manifold $M_{S S}$ (red curve) from Fig. 17 onto ( $V_{2}, V_{1}$ )-plane. The solid (respectively, dashed) cyan curves denote stable (respectively, unstable) periodic orbit branches. Magenta circles represent saddle-node bifurcations of $M_{S S}$ , in which the upper two have the same $V_{2}$ values as the folds of the $V_{2}$ -nullcline at green circle and triangle. Other colors and symbols are the same as in Fig. 17.

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

C. Effects of varying $ε$ and δ on MMOs

Unlike the sensitivity of MMOs at $g_{syn} = 4.3$ to timescale variations, the interaction of canard and DHB mechanisms due to the existence of the upper CDH when $g_{syn} = 4.4$ makes MMOs much more robust, as discussed above and illustrated in Fig. 3(b). Specifically, MMOs persist over biologically relevant ranges of $C_{1} \in (0.1, 80)$ and $ϕ_{2} \in (10^{- 4}, 0.01)$ (also see Figs. 22 and and2323).

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 $δ \geq O (\sqrt{ε})$ , DHB-like features diminish and the canard mechanism dominates.

Below we summarize the effects of increasing the singular perturbations:

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

1. Increasing C₁ 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 ϕ₂ makes the canard mechanism dominate

Increasing $ϕ_{2}$ speeds up the superslow variable $w_{2}$ and hence drives the three-timescale system (2) closer to (1F, 3S) splitting. As a result, the superslow manifold $M_{S S}$ 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 $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 $δ = O (ε)$ [e.g., Figs. 19 and 23(a)] exhibit both canard and DHB characteristics. When $δ \approx \frac{1}{2} \sqrt{ε}$ [Fig. 23(b)], there are still some DHB-like features. Further increasing $δ$ to $δ = 0.006 \approx \sqrt{ε}$ or $δ = 0.01 \approx ε^{\frac{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 $ϕ_{2}$ makes the canard mechanism dominate and the SAOs exhibit fewer DHB-like features. Conversely, decreasing $ϕ_{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. ³⁰^,²⁷^,²⁶^,⁴¹^,¹¹^,³⁹^,³⁸, and ⁴⁰). 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 $δ \to 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 $g_{syn} = 4.4$ , both $ε \to 0$ and $δ \to 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 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 ²⁶ 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. ²⁶ 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 $◯ 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 $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 $ε$ 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 $M_{S S}$ [see Figs. 15(b), 15(c), and and16].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.

TABLE II.

Effects of varying $ε$ and δ on MMOs.

	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₁ (or $ε$ )	MMOs are preserved with more DHB characteristics.	MMOs are preserved with more DHB characteristics.
Decreasing C₁ (or $ε$ )	Two MMOs/non-MMOs transitions are observed before a complete loss of MMOs.	MMOs are preserved and organized by both mechanisms.
Increasing ϕ₂ (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 ϕ₂ (or δ)	MMOs are preserved, but a non-monotonic effect on the small oscillations is observed.	MMOs are preserved with both DHB and canard characteristics.

Open in a separate window

As one would expect, MMOs when $g_{syn} = 4.3$ is also sensitive to increasing $δ$ , which speeds up the superslow variable and thus makes the DHB mechanism less relevant. Interestingly, however, increasing $δ$ does not just simply eliminate MMOs. Instead, it led to a total of five transitions between MMOs and non-MMO states [see 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 $g_{syn} = 4.3$ show strong robustness to increasing $ε$ or decreasing $δ$ . This is not surprising as both of these changes manifest the DHB mechanism by moving the three timescale problem closer to (3S, 1SS) separation. As a result, we observed more DHB characteristic in the MMOs as demonstrated in the case of increasing $ε$ [see Figs. 12(a) and and14].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 $g_{syn} = 4.3$ , an upper CDH occurs near the SAOs at $g_{syn} = 4.4$ . In this case, we showed that this CDH allowed for the coexistence and interaction of canard and DHB mechanisms, resulting in MMOs with strong robustness against timescale variations. Since both MMO solutions for $g_{syn} = 4.3$ and $4.4$ exhibit DHB mechanisms, they show similar responses and robustness to increasing $ε$ and decreasing $δ$ , both of which lead to more DHB characteristic in the MMOs. In contrast to $g_{syn} = 4.3$ , MMOs at $g_{syn} = 4.4$ are also robust against changes that speed up fast or superslow variables. This is due to the existence of the canard mechanism in addition to the DHB. Instead of eliminating the upper DHB in the case of $g_{syn} = 4.3$ [Fig. 10(a)], decreasing $ε$ 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 $ε \to 0$ singular limit. Hence, MMOs at $g_{syn} = 4.4$ persist as $ε$ decreases and are organized by both mechanisms. On the other hand, increasing $δ$ diminishes the relevance of the DHB mechanism and moves the trajectory further away from the superslow manifold, resulting in MMOs with more canard-like features.

To summarize, the coexistence of canard and DHB mechanisms for MMOs at $g_{syn} = 4.4$ due to the presence of a nearby upper CDH, leads to significantly enhanced robustness against timescale variations compared with $g_{syn} = 4.3$ , where only one mechanism is present. While we have not examined the case when only a canard mechanism is present, based on our analysis, we expect that MMOs with only canard mechanism would show more sensitivities to timescale variations that reduces the relevance of the canard mechanism such as increasing $ε$ or decreasing $δ$ . Similarly, such MMOs should show stronger robustness to decreasing $ε$ or increasing $δ$ . It would be of interest to explore such a scenario for future investigation. Furthermore, we did not notice any non-monotonic effects on the features of SAOs as $ϕ_{2}$ is decreased when $g_{syn} = 4.4$ , unlike what we observed in $g_{syn} = 4.3$ . This difference is likely attributed to the presence of an additional canard mechanism which may have hindered the occurrence of complex non-monotonic behaviors. A complete analysis of this phenomenon could be investigated in future work.

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.

APPENDIX A: DIMENSIONAL ANALYSIS

Here, we present a dimensional analysis of (2) to reveal its timescales. To this end, we introduce a dimensionless timescale $t_{s} = t / Q_{t}$ with reference timescale $Q_{t} = min (τ_{w} (V)) / ϕ_{1} = 18.86 ms = O (10) ms$ . This transforms (2) to the dimensionless system (4) given in the Introduction, namely,

\begin{matrix} \frac{C_{1}}{g_{\max} Q_{t}} \frac{d V_{1}}{d t_{s}} := ε \frac{d V_{1}}{d t_{s}} = f_{1} (V_{1}, w_{1}, V_{2}), \\ \frac{d w_{1}}{d t_{s}} = Q_{t} ϕ_{1} (w_{\infty} (V_{1}) - w_{1}) / τ_{w} (V_{1}) := g_{1} (V_{1}, w_{1}), \\ \frac{d V_{2}}{d t_{s}} = \frac{Q_{t} g_{\max}}{C_{2}} {\bar{f}}_{2} (V_{2}, w_{2}) := f_{2} (V_{2}, w_{2}), \\ \frac{d w_{2}}{d t_{s}} = Q_{t} ϕ_{2} (w_{\infty} (V_{2}) - w_{2}) / τ_{w} (V_{2}) := δ g_{2} (V_{2}, w_{2}), \end{matrix}

(A1)

where $ε = \frac{C_{1}}{g_{\max} Q_{t}} = 0.1$ , $δ = Q_{t} ϕ_{2} / \min (τ_{w} (V_{2})) = 0.053$ , $g_{\max} = 8 {mS / cm}^{2}$ ,

f_{1} (V_{1}, w_{1}, V_{2}) = \frac{(I_{1} - g_{Ca} m_{\infty} (V_{1}) (V_{1} - V_{Ca}) - g_{K} w_{1} (V_{1} - V_{K}) - g_{L} (V_{1} - V_{L}) - g_{syn} S (V_{2}) (V_{1} - V_{syn}))}{g_{\max}},

and

{\bar{f}}_{2} (V_{2}, w_{2}) = \frac{(I_{2} - g_{Ca} m_{\infty} (V_{2}) (V_{2} - V_{Ca}) - g_{K} w_{2} (V_{2} - V_{K}) - g_{L} (V_{2} - V_{L})}{g_{\max}} .

From system (A1), we can see that the voltage $V_{1}$ evolves on a fast timescale $(\frac{C_{1}}{g_{\max}} = 1 ms)$ , $(w_{1}, V_{2})$ evolve on a slow timescale $(\frac{\min (τ_{w} (V_{1}))}{ϕ_{1}} \approx \frac{C_{2}}{g_{\max}} = O (10) ms)$ , whereas $w_{2}$ evolves on a superslow timescale $(\frac{\min (τ_{w} (V_{2}))}{ϕ_{2}} = O (100) ms)$ .

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

To derive the conditions for FSN singularities, we note that the eigenvalues $λ$ of the folded singularities satisfy the following algebraic equation:

λ^{3} - tr (J_{D}) λ^{2} + \frac{1}{2} ({(tr (J_{D}))}^{2} - tr (J_{D}^{2})) λ - det (J_{D}) = 0,

(B1)

where $tr (J_{D})$ and $det (J_{D})$ denote the trace and determinant of the Jacobian matrix $J_{D}$ of the desingularized system (18), respectively. As discussed before, $det (J_{D}) \equiv 0$ along the folded singularity curve $M$ . This also directly follows from the following calculations:

\begin{aligned} det (J_{D}) & = det (\begin{matrix} F_{V_{1}} & F_{V_{2}} & F_{w_{2}} \\ G_{V_{1}} & G_{V_{2}} & G_{w_{2}} \\ H_{V_{1}} & H_{V_{2}} & H_{w_{2}} \end{matrix}) |_{M} \\ = det (\begin{matrix} F_{V_{1}} & F_{V_{2}} & F_{w_{2}} \\ G_{V_{1}} & G_{V_{2}} & 0 \\ H_{V_{1}} & H_{V_{2}} & 0 \end{matrix}) |_{M} \\ = F_{w_{2}} (G_{V_{1}} H_{V_{2}} - G_{V_{2}} H_{V_{1}}) \equiv 0, \end{aligned}

where the entries in $J_{D}$ denote partial derivatives. The last equality in the above equations holds because $G_{V_{1}} H_{V_{2}} \equiv G_{V_{2}} H_{V_{1}} = δ f_{2} g_{2} F_{1 V_{1} V_{2}} F_{1 V_{1} V_{1}}$ [see (18)], where $F_{1 V_{1} V_{1}} = \frac{\partial^{2} F_{1}}{\partial V_{1}^{2}}$ and $F_{1 V_{1} V_{2}} = \frac{\partial^{2} F_{1}}{\partial V_{1} \partial V_{2}}$ . Hence, the remaining two nontrivial eigenvalues $λ_{w}$ and $λ_{s}$ satisfy

λ^{2} - tr (J_{D}) λ + \frac{1}{2} ({(tr (J_{D}))}^{2} - tr (J_{D}^{2})) = 0.

It follows that an $FSN$ is given by

\frac{1}{2} {(tr (J_{D}))}^{2} - tr (J_{D}^{2}) = F_{V_{1}} G_{V_{2}} - F_{V_{2}} G_{V_{1}} - F_{w_{2}} H_{V_{1}} = 0.

Plugging in functions $F$ , $G$ , and $H$ from (18), we can rewrite the above condition as Eq. (21) in Sec. II B, namely,

\begin{aligned} FSN & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : f_{2} Q (V_{1}, V_{2}, w_{2}) \\ = δ P (V_{1}, V_{2}, w_{2})}, \end{aligned}

(B2)

where $Q (V_{1}, V_{2}, w_{2}) := F_{V_{1}} F_{1 V_{1} V_{2}} - F_{V_{2}} F_{1 V_{1} V_{1}}$ and $P (V_{1}, V_{2}, w_{2}) := g_{2} F_{1 V_{2}} f_{2 w_{2}} F_{1 V_{1} V_{1}}$ .

Next, we prove there are two ways an $FSN$ can occur. In case 1, suppose $Q \neq 0$ . An $FSN$ can occur when $f_{2} = δ K_{1} (V_{1}, V_{2}, w_{2})$ , where $K_{1} (V_{1}, V_{2}, w_{2}) := P / Q$ . Note that $FSN$ is a folded singularity point (20), that is, $F = g_{1} - F_{1 V_{2}} f_{2} = 0$ . This implies $g_{1} = F_{1 V_{2}} f_{2} = δ K_{2} (V_{1}, V_{2}, w_{2})$ , where $K_{2} (V_{1}, V_{2}, w_{2}) := F_{1 V_{2}} K_{1} (V_{1}, V_{2}, w_{2})$ . Thus, the first way that an $FSN$ can occur is described as the following:

\begin{aligned} {FSN}^{1} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : f_{2} = δ K_{1} (V_{1}, V_{2}, w_{2})} \\ = {(V_{1}, w_{1}, V_{2}, w_{2}) \in L_{s} : \\ f_{2} = δ K_{1} (V_{1}, V_{2}, w_{2}), g_{1} = δ K_{2} (V_{1}, V_{2}, w_{2})}, \end{aligned}

which is Eq. (22) in Sec. II B. This implies that an ${FSN}^{1}$ point becomes a CDH when $δ = 0$ .

In case 2, suppose $f_{2} \neq 0$ . An $FSN$ can occur when $Q (V_{1}, V_{2}, w_{2}) = δ K_{3} (V_{1}, V_{2}, w_{2})$ where $K_{3} (V_{1}, V_{2}, w_{2}) = P / f_{2}$ . As a result, the second way that an $FSN$ can occur is defined as

\begin{aligned} {FSN}^{2} & := {(V_{1}, w_{1}, V_{2}, w_{2}) \in M : Q (V_{1}, V_{2}, w_{2}) \\ = δ K_{3} (V_{1}, V_{2}, w_{2})}, \end{aligned}

which is Eq. (23) in Sec. II B. This implies that an ${FSN}^{2}$ is $O (δ)$ close to the intersection of $M$ and $Q (V_{1}, V_{2}, w_{2}) = 0$ .

APPENDIX C: CDH AT DOUBLE SINGULAR LIMIT EXHIBITS TWO LINEARLY INDEPENDENT CRITICAL EIGENVECTORS

Recall the Jacobian matrix $J_{D}$ at an ${FSN}^{1}$ , which becomes a CDH at the double singular limit, has two zero eigenvalues. We prove that the center subspace associated with the two zero eigenvalues is two dimensional and is transverse to the fold surface of the critical manifold (see Sec. II D and Fig. 4).

The Jacobian matrix $J_{D}$ of the desingularized system at the double singular limit ( $ε \to 0, δ \to 0)$ is given by

J_{D} = (\begin{matrix} F_{V_{1}} & F_{V_{2}} & F_{w_{2}} \\ G_{V_{1}} & G_{V_{2}} & 0 \\ 0 & 0 & 0 \end{matrix}),

in which $G_{V_{1}} = - F_{1 V_{1} V_{1}} f_{2}$ and $G_{V_{2}} = - F_{1 V_{1} V_{2}} f_{2}$ . At a CDH, $f_{2} = g_{1} = 0$ . It follows that $G_{V_{1}} = G_{V_{2}} = 0$ . Thus, the Jacobian matrix $J_{D}$ at a CDH singularity becomes

J_{D} = (\begin{matrix} F_{V_{1}} & F_{V_{2}} & F_{w_{2}} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}),

which has a nullity of 2, implying the existence of two linearly independent critical eigenvectors associated with zero eigenvalues. Moreover, the center subspace at the CDH, given by the plane $F_{V_{1}} V_{1} + F_{V_{2}} V_{2} + F_{w_{2}} w_{2} = 0$ [see pink planes in Figs. 4(c) and 4(d)], is transverse to the fold surface $L_{s}$ of the critical manifold. This is because the normal vector of $L_{s}$ at the CDH, which is given by $n_{f} = \frac{n_{s}}{| n_{s} |}$ where $n_{s} = (F_{1 V_{1} V_{1}}, F_{1 V_{1} V_{2}}, 0)$ is not parallel to the normal vector of the center subspace given by $(F_{V_{1}}, F_{V_{2}}, F_{w_{2}})$ . Thus, the CDH singularity (an ${FSN}^{1}$ at the double singular limit) considered in system (2) is a folded saddle-node singularity, with the center manifold of the ${FSN}^{1}$ transverse to the fold of the critical manifold.

Finally, we verify a previous claim that we made in Sec. II D that the eigenvector associated with the first trivial zero eigenvalue of $J_{D}$ (denoted by $v_{0}$ ) is always tangent to $L_{s}$ . Through direct computations, we obtain $v_{0} = \frac{u_{0}}{| u_{0} |}$ with

\begin{aligned} u_{0} & = (F_{V_{2}} F_{w_{2}} F_{1 V_{1} V_{2}}, - F_{V_{2}} F_{w_{2}} F_{1 V_{1} V_{1}}, \\ {(F_{V_{2}})}^{2} F_{1 V_{1} V_{1}} - F_{V_{1}} F_{V_{2}} F_{1 V_{1} V_{2}}) . \end{aligned}

The dot product of $v_{0}$ and the unit normal vector of $L_{s}$ at the CDH is given by

\begin{aligned} v_{0} \cdot n_{f} & = \frac{(F_{V_{2}} F_{w_{2}} F_{1 V_{1} V_{2}} F_{1 V_{1} V_{1}} - F_{V_{2}} F_{w_{2}} F_{1 V_{1} V_{1}} F_{1 V_{1} V_{2}})}{| u_{0} | | n_{s} |} = 0. \end{aligned}

It follows directly that $v_{0}$ is tangent to $L_{s}$ as expected.

APPENDIX D: DECREASING $ϕ_{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 $ϕ_{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 $ϕ_{2}$ , the number of SAOs will decrease and their amplitudes will increase (e.g., when $ϕ_{2}$ decreases from $0.0008$ to $0.0007$ ); (2) if no additional big spike is generated during the process, the changes to the SAOs will be reversed (e.g., there are more SAOs with smaller amplitudes when $ϕ_{2}$ decreases from $0.0007$ to $0.0006$ ). Below, we examine the two opposite effects using the projections of the corresponding solutions and bifurcation diagrams onto $(V_{2}, V_{1})$ space (see Fig. 24).

FIG. 24.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g024.jpg

Open in a separate window

Projections of solutions of (2) with $g_{syn} = 4.3$ , $C_{1} = 8$ , and different $ϕ_{2}$ values onto $(V_{2}, V_{1})$ -plane. (a) Projections of trajectories with $ϕ_{2} = 0.0008$ (black) and $ϕ_{2} = 0.0007$ (dark yellow). (b) Projections of trajectories with $ϕ_{2} = 0.0008$ (black) and $ϕ_{2} = 0.0006$ (green). The black and dark yellow stars, circles, squares, and triangles have the same meaning as the green symbols in Fig. 5. Other symbols and curves have the same meaning as in Fig. 9.

Figure 24(a) shows that the amplitude of SAOs with $ϕ_{2} = 0.0008$ (black) is smaller than SAOs with $ϕ_{2} = 0.0007$ (dark yellow). There are also more black SAOs when $ϕ_{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 $◯ 1$ . Thus, one extra full spike occurs for the yellow trajectory compared to the black trajectory during the jump at phase $◯ 2$ . As a result, the black trajectory, after its last big spike, approaches the black square at its maximum $V_{2}$ along $M_{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 $ϕ_{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 $◯ 1$ , similar to the case when $ϕ_{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 $ϕ_{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 $ϕ_{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 $ϕ_{2}$ is decreased (see Appendix D).

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

FIG. 25.

An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g025.jpg

Open in a separate window

Projections of solutions (2) with $g_{syn} = 4.3$ , $C_{1} = 8$ ( $ε = 0.1$ ) and different values of $ϕ_{2}$ onto the $(V_{2}, V_{1})$ -space. (a)–(c) Projections of trajectories with $ϕ_{2} = 0.0012$ , $ϕ_{2} = 0.0016$ and $ϕ_{2} = 0.0018$ are shown by green, dark yellow, and blue curves, respectively. Other symbols and color codings have the same meaning as in Fig. 9.

In contrast, increasing $ϕ_{2}$ from $0.0012$ to $0.0016$ does not change the number of large spikes before phase $◯ 3$ [compare the green and dark yellow trajectories in Figs. 25(a) and 25(b)]. However, speeding up of $ϕ_{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.

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.

REFERENCES

1. Desroches M., Guckenheimer J., Krauskopf B., Kuehn C., Osinga H. M., and Wechselberger M., “Mixed-mode oscillations with multiple time scales,” SIAM Rev. 54(2), 211–288 (2012). 10.1137/100791233 [CrossRef] [Google Scholar]

2. Awal N. M., Epstein I. R., Kaper T. J., and Vo T., “Symmetry-breaking rhythms in coupled, identical fast–slow oscillators,” Chaos 33(1), 011102 (2023). 10.1063/5.0131305 [PubMed] [CrossRef] [Google Scholar]

3. Hudson J. L., Hart M., and Marinko D., “An experimental study of multiple peak periodic and nonperiodic oscillations in the Belousov–Zhabotinskii reaction,” J. Chem. Phys. 71(4), 1601–1606 (1979). 10.1063/1.438487 [CrossRef] [Google Scholar]

4. Battaglin S. and Pedersen M. G., “Geometric analysis of mixed-mode oscillations in a model of electrical activity in human beta-cells,” Nonlinear Dyn. 104(4), 4445–4457 (2021). 10.1007/s11071-021-06514-z [CrossRef] [Google Scholar]

5. Curtu R., “Singular Hopf bifurcation and mixed-mode oscillations in a two-cell inhibitory neural network,” Physica D 239(9), 504–514 (2010). 10.1016/j.physd.2009.12.010 [CrossRef] [Google Scholar]

6. Curtu R. and Rubin J., “Interaction of canard and singular Hopf mechanisms in a neural model,” SIAM J. Appl. Dyn. Syst. 10(4), 1443–1479 (2011). 10.1137/110823171 [CrossRef] [Google Scholar]

7. Harvey E., Kirk V., Wechselberger M., and Sneyd J., “Multiple timescales, mixed mode oscillations and canards in models of intracellular calcium dynamics,” J. Nonlinear Sci. 21, 639–683 (2011). 10.1007/s00332-011-9096-z [CrossRef] [Google Scholar]

8. Kimrey J., Vo T., and Bertram R., “Big ducks in the heart: Canard analysis can explain large early afterdepolarizations in cardiomyocytes,” SIAM J. Appl. Dyn. Syst. 19(3), 1701–1735 (2020). 10.1137/19M1300777 [CrossRef] [Google Scholar]

9. Kimrey J., Vo T., and Bertram R., “Canard analysis reveals why a large

{Ca}^{2 +}

window current promotes early afterdepolarizations in cardiac myocytes,” PLoS Comput. Biol. 16(11), e1008341 (2020). 10.1371/journal.pcbi.1008341 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

10. Krupa M., Popović N., Kopell N., and Rotstein H. G., “Mixed-mode oscillations in a three time-scale model for the dopaminergic neuron,” Chaos 18(1), 015106 (2008). 10.1063/1.2779859 [PubMed] [CrossRef] [Google Scholar]

11. Krupa M., Vidal A., Desroches M., and Clément F., “Mixed-mode oscillations in a multiple time scale phantom bursting system,” SIAM J. Appl. Dyn. Syst. 11(4), 1458–1498 (2012). 10.1137/110860136 [CrossRef] [Google Scholar]

12. Kügler P., Erhardt A. H., and Bulelzai M. A. K., “Early afterdepolarizations in cardiac action potentials as mixed mode oscillations due to a folded node singularity,” PLoS One 13(12), e0209498 (2018). 10.1371/journal.pone.0209498 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

13. Pavlidis E., Campillo F., Goldbeter A., and Desroches M., “Multiple-timescale dynamics, mixed mode oscillations and mixed affective states in a model of bipolar disorder,” Cognit. Neurodyn. (published online) (2022). 10.1007/s11571-022-09900-4 [CrossRef] [Google Scholar]

14. Teka W., Tabak J., and Bertram R., “The relationship between two fast/slow analysis techniques for bursting oscillations,” Chaos 22(4), 043117 (2012). 10.1063/1.4766943 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

15. Vo T., Bertram R., Tabak J., and Wechselberger M., “Mixed mode oscillations as a mechanism for pseudo-plateau bursting,” J. Comput. Neurosci. 28, 443–458 (2010). 10.1007/s10827-010-0226-7 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

16. Vo T., Tabak J., Bertram R., and Wechselberger M., “A geometric understanding of how fast activating potassium channels promote bursting in pituitary cells,” J. Comput. Neurosci. 36, 259–278 (2014). 10.1007/s10827-013-0470-8 [PubMed] [CrossRef] [Google Scholar]

17. Yu N., Kuske R., and Li Y. X., “Stochastic phase dynamics and noise-induced mixed-mode oscillations in coupled oscillators,” Chaos 18(1), 015112 (2008). 10.1063/1.2790369 [PubMed] [CrossRef] [Google Scholar]

18. Fenichel N., “Geometric singular perturbation theory for ordinary differential equations,” J. Differ. Equ. 31(1), 53–98 (1979). 10.1016/0022-0396(79)90152-9 [CrossRef] [Google Scholar]

19. Szmolyan P. and Wechselberger M., “Canards in

R^{3}

,” J. Differ. Equ. 177(2), 419–453 (2001). 10.1006/jdeq.2001.4001 [CrossRef] [Google Scholar]

20. Wechselberger M., “Existence and bifurcation of canards in

R^{3}

in the case of a folded node,” SIAM J. Appl. Dyn. Syst. 4(1), 101–139 (2005). 10.1137/030601995 [CrossRef] [Google Scholar]

21. Baer S. M., Erneux T., and Rinzel J., “The slow passage through a Hopf bifurcation: Delay, memory effects, and resonance,” SIAM J. Appl. Math. 49(1), 55–71 (1989). 10.1137/0149003 [CrossRef] [Google Scholar]

22. Hayes M. G., Kaper T. J., Szmolyan P., and Wechselberger M., “Geometric desingularization of degenerate singularities in the presence of fast rotation: A new proof of known results for slow passage through Hopf bifurcations,” Indag. Math. 27(5), 1184–1203 (2016). 10.1016/j.indag.2015.11.005 [CrossRef] [Google Scholar]

23. Neishtadt A., “On delayed stability loss under dynamical bifurcations I,” Differ. Equ. 23, 1385–1390 (1987). [Google Scholar]

24. Neishtadt A., “On delayed stability loss under dynamical bifurcations II,” Differ. Equ. 24, 171–176 (1988). [Google Scholar]

25. De Maesschalck P., Kutafina E., and Popović N., “Three time-scales in an extended Bonhoeffer–van der Pol oscillator,” J. Dyn. Differ. Equ. 26, 955–987 (2014). 10.1007/s10884-014-9356-3 [CrossRef] [Google Scholar]

26. Letson B., Rubin J. E., and Vo T., “Analysis of interacting local oscillation mechanisms in three-timescale systems,” SIAM J. Appl. Dyn. Syst. 77(3), 1020–1046 (2017). 10.1137/16M1088429 [CrossRef] [Google Scholar]

27. Vo T., Bertram R., and Wechselberger M., “Multiple geometric viewpoints of mixed mode dynamics associated with pseudo-plateau bursting,” SIAM J. Appl. Dyn. Syst. 12(2), 789–830 (2013). 10.1137/120892842 [CrossRef] [Google Scholar]

28. Baldemir H., Avitabile D., and Tsaneva-Atanasova K., “Pseudo-plateau bursting and mixed-mode oscillations in a model of developing inner hair cells,” Commun. Nonlinear Sci. Numer. Simul. 80, 104979 (2020). 10.1016/j.cnsns.2019.104979 [CrossRef] [Google Scholar]

29. Chumakov G. A., Chumakova N. A., and Lashina E. A., “Modeling the complex dynamics of heterogeneous catalytic reactions with fast, intermediate, and slow variables,” Chem. Eng. J. 282, 11–19 (2015). 10.1016/j.cej.2015.03.017 [CrossRef] [Google Scholar]

30. Jalics J., Krupa M., and Rotstein H. G., “Mixed-mode oscillations in a three time-scale system of ODEs motivated by a neuronal model,” Dyn. Syst. 25(4), 445–482 (2010). 10.1080/14689360903535760 [CrossRef] [Google Scholar]

31. Perryman C. and Wieczorek S., “Adapting to a changing environment: Non-obvious thresholds in multi-scale systems,” Proc. R. Soc. A 470(2170), 20140226 (2014). 10.1098/rspa.2014.0226 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

32. Wang Y. and Rubin J. E., “Multiple timescale mixed bursting dynamics in a respiratory neuron model,” J. Comput. Neurosci. 41, 245–268 (2016). 10.1007/s10827-016-0616-6 [PubMed] [CrossRef] [Google Scholar]

33. Wang Y. and Rubin J. E., “Timescales and mechanisms of sigh-like bursting and spiking in models of rhythmic respiratory neurons,” J. Math. Neurosci. 7, 1–39 (2017). 10.1186/s13408-017-0045-5 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

34. Wang Y. and Rubin J. E., “Complex bursting dynamics in an embryonic respiratory neuron model,” Chaos 30(4), 043127 (2020). 10.1063/1.5138993 [PubMed] [CrossRef] [Google Scholar]

35. Nan P., Wang Y., Kirk V., and Rubin J. E., “Understanding and distinguishing three-time-scale oscillations: Case study in a coupled Morris-Lecar system,” SIAM J. Appl. Dyn. Syst. 14(3), 1518–1557 (2015). 10.1137/140985494 [CrossRef] [Google Scholar]

36. De Maesschalck P., Kutafina E., and Popović N., “Sector-delayed-Hopf-type mixed-mode oscillations in a prototypical three-time-scale model,” Appl. Math. Comput. 273, 337–352 (2016). 10.1016/j.amc.2015.09.083 [CrossRef] [Google Scholar]

37. Desroches M. and Kirk V., “Spike-adding in a canonical three-time-scale model: Superslow explosion and folded-saddle canards,” SIAM J. Appl. Dyn. Syst. 17(3), 1989–2017 (2018). 10.1137/17M1143411 [CrossRef] [Google Scholar]

38. Kaklamanos P. and Popović N., “Complex oscillatory dynamics in a three-timescale El Niño Southern Oscillation model,” Physica D 449, 133740 (2023). 10.1016/j.physd.2023.133740 [CrossRef] [Google Scholar]

39. Kaklamanos P., Popović N., and Kristiansen K. U., “Bifurcations of mixed-mode oscillations in three-timescale systems: An extended prototypical example,” Chaos 32(1), 013108 (2022). 10.1063/5.0073353 [PubMed] [CrossRef] [Google Scholar]

40. Kaklamanos P., Popović N., and Kristiansen K. U., “Geometric singular perturbation analysis of the multiple-timescale Hodgkin-Huxley equations,” SIAM J. Appl. Dyn. Syst. 22(3), 1552–1589 (2023). 10.1137/22M1477477 [CrossRef] [Google Scholar]

41. Krupa M., Popović N., and Kopell N., “Mixed-mode oscillations in three time-scale systems: A prototypical example,” SIAM J. Appl. Dyn. Syst. 7(2), 361–420 (2008). 10.1137/070688912 [CrossRef] [Google Scholar]

42. Sadhu S., “Complex oscillatory patterns in a three-timescale model of a generalist predator and a specialist predator competing for a common prey,” Discrete Contin. Dyn. Syst. Ser. B 28(5), 3014–3051 (2023). 10.3934/dcdsb.2022202 [CrossRef] [Google Scholar]

43. Morris C. and Lecar H., “Voltage oscillations in the barnacle giant muscle fiber,” Biophys. J. 35(1), 193–213 (1981). 10.1016/S0006-3495(81)84782-0 [PMC free article] [PubMed] [CrossRef] [Google Scholar]

44. Rinzel J. and Ermentrout G. B., “Analysis of neural excitability and oscillations,” Methods Neuronal Model. 2, 251–292 (1998). [Google Scholar]

45. Krupa M. and Wechselberger M., “Local analysis near a folded saddle-node singularity,” J. Differ. Equ. 248(12), 2841–2888 (2010). 10.1016/j.jde.2010.02.006 [CrossRef] [Google Scholar]

46. Roberts K. L., Rubin J. E., and Wechselberger M., “Averaging, folded singularities, and torus canards: Explaining transitions between bursting and spiking in a coupled neuron model,” SIAM J. Appl. Dyn. Syst. 14(4), 1808–1844 (2015). 10.1137/140981770 [CrossRef] [Google Scholar]

47. Vo T. and Wechselberger M., “Canards of folded saddle-node type I,” SIAM J. Math. Anal. 47(4), 3235–3283 (2015). 10.1137/140965818 [CrossRef] [Google Scholar]

48. Guckenheimer J. and Willms A. R., “Asymptotic analysis of subcritical Hopf-homoclinic bifurcation,” Physica D 139(3-4), 195–216 (2000). 10.1016/S0167-2789(99)00225-0 [CrossRef] [Google Scholar]

49. Brøns M., Krupa M., and Wechselberger M., “Mixed mode oscillations due to the generalized canard phenomenon,” Fields Inst. Commun. 49, 39–63 (2006). [Google Scholar]

50. Rubin J. E. and Wechselberger M., “The selection of mixed-mode oscillations in a Hodgkin-Huxley model with multiple timescales,” Chaos 18(1), 015105 (2008). 10.1063/1.2789564 [PubMed] [CrossRef] [Google Scholar]

Articles from Chaos are provided here courtesy of American Institute of Physics

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

Ngoc Anh Phan

Yangyang Wang

Associated Data

Abstract

I. INTRODUCTION

TABLE I.

II. GEOMETRIC SINGULAR PERTURBATION ANALYSIS

A. Singular limits

1. ε→0 Singular limit

2. δ → 0 Singular limit

3. ε→0,δ→0 Double singular limits

B. Slow reduced problem and canard dynamics

1. Folded saddle node (FSN)

C. Slow layer problem and delayed Hopf bifurcations

D. Interaction between canard and delayed Hopf

E. Singular orbit construction

III. ANALYSIS OF MMOs WHEN gsyn = 4.3 mS/cm2

A. Relation of the trajectory to MS and MSS

1. Canard mechanism does not contribute to MMOs

2. MMOs arise from the delayed Hopf mechanism

B. Effects of varying ε and δ on DHB and FSN points

1. Effect of ε on DHB

2. Effect of δ on FSN points

C. Why there are no SAOs near the lower CDH

D. Effects of varying ε and δ on MMOs

1. MMOs are robust to increasing C1

2. Decreasing C1 leads to three MMOs/non-MMOs transitions

3. Decreasing ϕ2 preserves MMOs but causes non-monotonic effects on SAOs

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

IV. ANALYSIS OF MMOS WHEN gsyn = 4.4 mS/cm2

A. Relation of the trajectory to MS and MSS

B. MMOs are organized by both canard and DHB mechanisms

1. Canard dynamics

2. Delayed Hopf bifurcation mechanism

C. Effects of varying ε and δ on MMOs

1. Increasing C1 makes DHB dominate

2. Increasing ϕ2 makes the canard mechanism dominate

V. DISCUSSION

TABLE II.

ACKNOWLEDGMENTS

APPENDIX A: DIMENSIONAL ANALYSIS

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

APPENDIX C: CDH AT DOUBLE SINGULAR LIMIT EXHIBITS TWO LINEARLY INDEPENDENT CRITICAL EIGENVECTORS

APPENDIX D: DECREASING ϕ2 WHEN gsyn=4.3 CAUSES NON-MONOTONIC EFFECTS ON SAOS

FIG. 24.

APPENDIX E: INCREASING ϕ2 WHEN gsyn=4.3 LEADS TO MULTIPLE MMOS/NON-MMOS TRANSITIONS

FIG. 25.

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

REFERENCES

1. $ε \to 0$ Singular limit

3. $ε \to 0, δ \to 0$ Double singular limits

III. ANALYSIS OF MMOs WHEN g_syn = 4.3 mS/cm²

A. Relation of the trajectory to M_S and M_SS

B. Effects of varying $ε$ and δ on DHB and FSN points

1. Effect of $ε$ on DHB

D. Effects of varying $ε$ and δ on MMOs

1. MMOs are robust to increasing C₁

2. Decreasing C₁ leads to three MMOs/non-MMOs transitions

3. Decreasing ϕ₂ preserves MMOs but causes non-monotonic effects on SAOs

4. Increasing ϕ₂ leads to five MMOs/non-MMOs transitions

IV. ANALYSIS OF MMOS WHEN g_syn = 4.4 mS/cm²

A. Relation of the trajectory to M_S and M_SS

C. Effects of varying $ε$ and δ on MMOs

1. Increasing C₁ makes DHB dominate

2. Increasing ϕ₂ makes the canard mechanism dominate

APPENDIX D: DECREASING $ϕ_{2}$ WHEN $g_{syn} = 4.3$ CAUSES NON-MONOTONIC EFFECTS ON SAOS

APPENDIX E: INCREASING $ϕ_{2}$ WHEN $g_{syn} = 4.3$ LEADS TO MULTIPLE MMOS/NON-MMOS TRANSITIONS