Learn more: PMC Disclaimer | PMC Copyright Notice
Mixed-mode oscillations in a three-timescale coupled Morris–Lecar system
Associated Data
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;1 these are complex oscillatory dynamics characterized by the concatenation of small-amplitude oscillations (SAOs) and large-amplitude excursions in each periodic cycle. Such phenomena have been recognized in many branches of sciences including physics, chemistry,2,3 and particularly life sciences.4–17
Theoretical analysis of MMOs in systems with two distinct timescales has been well developed with the implementation of the geometric singular perturbation theory (GSPT);18 see Ref. 1 for review. Two common mechanisms leading to the occurrence of MMOs in multiple timescale problems are canard dynamics associated with the twisting of slow manifolds due to folded singularities19,20 and a slow passage through the delayed Andronov–Hopf bifurcation (DHB) of the fast subsystem.21–24 While in two-timescale settings, these two mechanisms remain separated, they can coexist and interact in three-timescale regime.14,25–27
Compared with the extremely well-studied MMOs in two-timescale problems, the theory of MMOs in the three-timescale settings has been less well developed. Traditionally, three-timescale problems are simplified to two-timescale problems, which is the natural setting for geometric singular perturbation theory.28 However, many real-world systems have more than two timescales.25,27,29–34 It has also been established that a two-timescale decomposition fails to capture certain aspects of the system’s dynamics.35 Therefore, classifying three timescales into two groups is not a sufficient approach for modeling and analysis.
MMOs in three-timescale systems have been studied before (see, e.g., Refs. 10, 30, 11, 27, 25, 26, and 36–42 and references therein). Initial approaches were to consider three-dimensional systems,
with special cases or .10,25,30,36,41 MMOs were shown to emerge through an effect analogous to a slow passage through a canard explosion.25,30,41 More recently, there has been a growing interest in MMOs with independent singular perturbation parameters and , as explored in various three-dimensional models.26,38–40 In particular, Ref. 26 centered on a novel singularity type denoted as canard-delayed-Hopf (CDH) singularity, which naturally arises in three-timescale settings when the two mechanisms for MMOs (the fast subsystem Hopf and a folded node) coexist and interact. The authors investigated the existence and properties of MMOs near the CDH singularity.
In this paper, we contribute to the investigation of MMOs in three-timescale settings by considering a model of four-dimensional coupled Morris–Lecar neurons43,44 that was introduced by Ref. 35. The model equations are given by
with
Table I lists the parameter values for the model chosen to ensure that (2) exhibits three distinct timescales, where is fast, are slow, and is superslow. In a more biologically realistic model for calcium and voltage interactions, might represent membrane potential, while might represent intracellular calcium concentration with appropriate adjustments to parameter units and functional terms (see, e.g., Refs. 32 and 33). For the physiological description of functions in (2) and (3), we refer readers to Ref. 35 for details.
TABLE I.
Parameter values | |||||
---|---|---|---|---|---|
C 1 | 8 μF/cm2 | I 1 | 0 μA/cm2 | ϕ 1 | 0.01 |
C 2 | 100 μF/cm2 | I 2 | 60 μA/cm2 | ϕ 2 | 0.001 |
V Ca | 120 mV | g Ca | 4 mS/cm2 | K 1 | −1.2 mV |
V K | −84 mV | g K | 8 mS/cm2 | K 2 | 18 mV |
V L | −60 mV | g L | 2 mS/cm2 | K 3 | 12 mV |
V syn | 30 mV | θ s | −20 mV | K 4 | 17.4 mV |
β | 0.5 ms−1 | σ s | 10 mV |
In the absence of coupling , is excitable with an attracting critical point at relatively low value, whereas is oscillatory with an attracting limit cycle independent of the value of and the dynamics of ). To analyze the three-timescale coupled Morris–Lecar neurons, Ref. 35 extended two approaches previously developed in the context of GSPT for the analysis of two-timescale systems to the three-timescale setting and showed these two approaches complemented each other nicely. By varying in system (2), the authors identified various solution features that truly require three timescales, thus demonstrating the functional relevance of three timescales in the model. While system (2) exhibits both the fast subsystem Hopf and folded nodes that can support MMOs, MMOs were not observed within the parameter regime examined by Ref. 35 (e.g., and in Fig. 1).
The goal of this work is the analysis of MMOs and their robustness in three-timescale systems by focusing on a coupled Morris–Lecar system (2). Based on our simulations, we have selected two MMO solutions on which to focus our analysis. Specifically, we consider and (with the unit of ), as highlighted in Fig. 1. To provide further insight into our choice of values, we perform a bifurcation analysis to explore the effect of on a singularity called canard-delayed-Hopf (CDH) that was first introduced by Ref. 26 (see Fig. 2, blue). As noted above, this singularity plays a crucial role in organizing MMOs within the three-timescale setting. We show in Sec. II that the full system (2) may exhibit two CDH points—one at larger values, (denoted as upper CDH) and the other at smaller (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.
Bifurcation curves of DHB (red) and CDH singularities (blue) for (2) with respect to . Specifically, these curves represent the upper DHB and upper CDH, corresponding to larger and values. The lower CDH and DHB, associated with smaller values, are not presented here. The two vertical asymptotes are given by and .
Projections to -space of the critical manifold fold surfaces (blue surface) for [(a) and (c)] and [(b) and (d)] . Also shown are the curves of folded singularities including folded node (solid green), folded saddle (dashed green), and two types of folded saddle-nodes (blue star: ; cyan star: ). 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 , an point (blue star) is close to a CDH singularity (blue diamond), whereas an (cyan star) is far away from any CDH. In the lower two panels [(c) and (d)] at the singular limit , the point becomes a CDH singularity (blue star overlapping with blue diamond). The center subspace of an (respectively, an ) is denoted by a pink plane (respectively, a yellow plane). It follows that the center manifolds of both and are transverse to .
In Fig. 2, we plot the bifurcation curves of the upper CDH and the upper DHB in the plane, both of which approach vertical asymptotes as decreases. Our first choice of represents values between the two asymptotes, at which the upper DHB exists but there is no upper CDH. In contrast, 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 (e.g., and as considered in Ref. 35), the system (2) does not produce MMOs. When , MMOs may arise through mechanisms similar to those we will thoroughly examine for 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 , 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 spikes to a plateau, as detailed in Ref. 35. Therefore, as 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 values (e.g., ), where both CDH and DHB points fall to low values that lie beyond the system’s trajectory.
In this paper, we focus only on MMOs at and . 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:
where , , is the slow dimensionless time variable, , , , and are functions specified in (A1) in Appendix A, which include details of the nondimensionalization procedure. For simplicity, we did not rescale and 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 evolves on a timescale of , on and on . Introducing a superslow time yields an equivalent description of dynamics,
which evolves on the superslow timescale and is called the superslow system. Alternatively, defining a fast time , we obtain the following fast system:
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 and 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 An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g003.jpg](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11087137/bin/CHAOEH-000034-053119_1-g003.jpg)
Regions of MMOs (yellow) and non-MMOs (blue) of the full system (2) in -space for (A) and (B) . Increasing slows down the fast variable , whereas increasing speeds up the superslow variable . The timescales of and remain unaffected. The red star marks the default parameter values of and as given in Table I. (a) . While MMOs are robust to increasing and decreasing , decreasing or increasing leads to multiple transitions between MMOs and non-MMOs (crossings between the dashed lines with the yellow/blue boundaries). (b) . MMOs are robust to changes of both and over the ranges of and . Note that the MMOs at will eventually vanish for and large enough at which there is no more timescale separation (data not shown).
When , 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 singular limit provides a faithful prediction for the observed MMOs, whereas the limit does not. This observation pinpoints the DHB from the limit as the only mechanism for the MMOs at . As a result, these MMO dynamics are sensitive to variations in timescales [see Fig. 3(a)]. Specifically, we show that MMOs persist for . Increasing via increasing or decreasing via decreasing to a degree where leads to multiple MMOs/non-MMOs transitions. MMOs are completely lost when for which the DHB is no longer present.
In contrast, there exists a CDH in the middle of the SAOs when 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 and singular limits contribute faithful predictions for the observed dynamics, resulting in MMOs with significantly stronger robustness than [Fig. 3(b)]. We demonstrate that when , MMOs exhibit both canard and DHB features. Upon tuning , DHB-like features disappear and the canard mechanism dominates. In summary, our findings reveal that MMOs near a CDH exhibit stronger robustness compared to those governed by a single mechanism, and that the CDH singularity is a determining factor in whether the two distinct MMO mechanisms can interact or not. However, it is essential to note that not all CDH singularities can support local MMOs. Specifically, we demonstrate that the lower CDH does not produce any local oscillations.
Our work is novel in two main aspects. First, to the best of our knowledge, our study is the first to investigate the geometric conditions that lead to robust occurrences of MMOs in three-timescale systems. It is worth noting that while Ref. 39 also considered the robustness of MMOs in a three-timescale system, their focus was specifically on MMOs with double epochs of SAOs. Second, we discovered that the CDH singularities do not always enable the two MMO mechanisms to interact and produce MMOs. This is different from past studies26,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 . 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 to variations in perturbation sizes or [i.e., varying and 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 . 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 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 analysis18,35 to the three-timesale coupled Morris–Lecar system (4) by treating as the only singular perturbation parameter,19,20 treating as the only singular perturbation parameter,21–24 and finally treating and as two independent singular perturbation parameters.26,27,35
Although the detailed GSPT analysis and derivation of subsystems have been previously presented in Ref. 35, we provide a brief overview in this paper for the sake of completeness. However, the focus of our current work is on the investigation of MMOs, which is distinct from the emphasis of Ref. 35. Specifically, we concentrate on reviewing and discussing the canard mechanism in Sec. II B, delayed Hopf bifurcation in Sec. II C, and their interactions in Sec. II D.
A. Singular limits
1. Singular limit
Fixing and taking in the fast system (6) yields the one-dimensional (1D) fast layer problem, a system that describes the dynamics of the fast variable, , for fixed values of the other variables,
The set of equilibrium points of the fast layer problem is called the critical manifold and is denoted as ,
Although is a three-dimensional (3D) manifold in space, it does not depend on . We can solve for as a function of and and can, therefore, represent as
for a function . It is well known that, for sufficiently small , normally hyperbolic parts of each perturb to a locally invariant manifold called a slow manifold, on which is given by an -perturbation of (Ref. 18); we simply use as a convenient numerical approximation of these slow manifolds.
is a 3D folded manifold with two-dimensional (2D) fold surface, , given by
or equivalently
where and denote the partial derivatives of and with respect to . The fold surface divides the critical manifold into attracting upper and lower branches where and repelling middle branch where .
Taking the same limit, i.e., with , in the slow system (4) yields the 3D slow reduced problem, a system that describes the dynamics of along ,
where .
2. δ → 0 Singular limit
Alternatively, fixing and taking in the slow system (4) yields the 3D slow layer problem in the form
where the superslow variable 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 ,
is a 1D subset of . Similar to , the normally hyperbolic parts of 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 where Fenichel’s theory (GSPT) breaks down.
Taking the same singular limit in the superslow system (5) leads to the 1D superslow reduced problem
where . The superslow motions of trajectories of (15) are slaved to until nonhyperbolic points are reached.
3. Double singular limits
Both the slow reduced problem (12) and the slow layer problem (13) still include two distinct timescales. Further taking the limit in (12) or taking the limit in (13) yields the same slow reduced layer problem
which describes the slow motion along and the superslow variable 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. 35.
B. Slow reduced problem and canard dynamics
To investigate canard dynamics, we project the slow reduced problem (12) onto to obtain a complete description of the dynamics along . To this end, we differentiate the graph representation of given by to obtain
where . Note that the reduced system (17) is singular at the fold surfaces (10). Nonetheless, this singular term can be removed by a time rescaling , and we obtain the following desingularized system:
We observe that the desingularized system (18) is equivalent to (17) on the attracting branch, i.e., for , but has the opposite orientation on the repelling branch, i.e., for .
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
For the chosen parameter set in Table I, always lies on the repelling branch of 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 defined by
Folded singularities are special points that allow trajectories of (17) to cross the fold with nonzero speed. Such solutions are called singular canards.19 Note that when projecting to -space, the condition is redundant and the fold surfaces become curves that overlap with the folded singularity curves .
Since there is a curve of folded singularities, the Jacobian of (18) evaluated along (denoted as , 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 ( ) where 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 . 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 form a two-dimensional trapping region (the funnel) on the attracting branch of such that all solutions in the funnel converge to that folded node. The funnel family of all folded nodes of and the fold surface 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 to a repelling 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, , becomes zero. This folded singularity is referred to as a “folded saddle node” ( ) and is characterized by the condition
where and are defined in (B2) in Appendix B, which contains a detailed derivation of the FSN condition. Similar to Ref. 27, we demonstrate in Appendix B that our system can exhibit an in two different ways:
or
where , and are defined in Appendix B. We discuss below in Remark II.1 that both and are novel types of .26
In our parameter regime, the ordinary singularity point lies in the middle branch of critical manifold and is not involved in any bifurcations of the folded singularities. Hence, the 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 and are transverse to the fold surface at the singularity points [see Figs. 4(c) and 4(d), yellow and pink planes]. Following Ref. 26, we identify and as novel types of .
The condition (22) suggests that an is close to the intersection point of the superslow manifold and the fold surface , which was defined as the canard-delayed-Hopf (CDH) singularity in Ref. 26. In contrast, is always far away from a CDH. Figure 4 shows the positions of (blue star), (cyan star), and CDH (blue diamond) in -space, for (top panels) and the singular limit (bottom panels). Recall the bifurcation of the upper CDH with respect to is shown in Fig. 2, which explains why there is no upper CDH for . It is worth noting that a CDH point of (2) is always a folded singularity because the critical manifold does not depend on the superslow variable .
C. Slow layer problem and delayed Hopf bifurcations
In this subsection, we turn to the slow layer problem (13) resulting from the singular limit, which exhibits delayed Hopf bifurcations that allow for interesting dynamics.
Let denote the Jacobian matrix of (13) evaluated along the superslow manifold , which is given by
where the nonzero entries denote partial derivatives.
The eigenvalues of are given by and the eigenvalues of
Thus, the Hopf bifurcation points on are given by and . The former defining condition can be rewritten as
It follows from (25) that an is close to the intersection of and , i.e., a CDH singularity. The subsystem HB bifurcation is also known as delayed Hopf bifurcation (DHB).
The isolated fold bifurcation points on are located by letting . That is,
The fold points that satisfy the former condition (denoted as ) are the folds of the -nullcline (see Fig. 5, green circle and green triangle), which correspond to the transition between superslow dynamics along and slow jumps. given by the latter condition (denoted as ) corresponds to the actual fold of when projected to -space. Since , and , it follows that lies on the middle branch of ( ) and hence will not play a role in dynamics. At the double singular limit , the fold point of 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.
Projection of the singular orbit (green) and the solution trajectory (black) of the full system (2) onto the phase plane of ( ) system with parameters given in Table I. The red curve is the -nullcline, which is the projection of the superslow manifold . 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 -nullcline. The circled numbers indicate four phases of the oscillations: superslow excursions along during and and slow jumps at the fold of during and .
D. Interaction between canard and delayed Hopf
To investigate the interaction between the canard and the delayed Hopf mechanisms in the double limit case ( ), we need to examine the slow reduced layer problem (16). The corresponding desingularized system is given by
where and are defined in (18). Note that (27) is the limit of the desingularized system (18) from Sec. II B. The folded singularities of (27) are exactly the same as given by (20), whereas the ordinary singularities of (27) are relaxed to be . The FSN condition at the double singular limit can be obtained from letting in the condition (22) or the condition (23). This implies that, at the double singular limit, an FSN1 becomes a CDH singularity [see Figs. 4(c) and 4(d)],
whereas an singularity is characterized by
According to Remarks II.2 and II.3, an singularity from the viewpoint converges to a CDH as and a DHB from the viewpoint converges to a CDH as . It is natural to expect that a CDH singularity point should serve as the interplay between the canard dynamics and the delayed Hopf bifurcation to produce MMOs, as seen in Ref. 26. However, this is not always the case. Specifically, we find that while the CDH on the upper fold surface [see Figs. 4(b) and 4(d)] supports MMOs with a high level of robustness due to the coexistence and interaction of two distinct MMO mechanisms, no MMO dynamics were observed near the lower CDH. It is worth highlighting that both the upper and lower CDH points in our system are of the same type as the CDH investigated in Ref. 26.
The CDH points in system (2) are singularities at the double singular limits. We prove in Appendix C that the CDH singularity in (2) is a novel type of folded saddle-node singularity as described in Ref. 26, with the center manifold of the CDH transverse to the fold 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 (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 , 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 . However, MMOs are not observed in the neighborhood of any lower CDH (see Sec. III for more discussions).
E. Singular orbit construction
To understand the dynamics of (2) using GSPT,18 we need to construct singular periodic orbits by concatenating solution segments of singular limit systems. According to GSPT, such a singular oscillation will generically perturb to a periodic solution of the full problem as we move away from the singular limit.
We now use our analysis from Subsections II A–II D to construct singular approximations of (2) in order to understand the full system trajectory . We let and denote trajectories of the fast layer problem (7) and the slow reduced problem (12) obtained from the singular limit. Let and denote trajectories of the slow layer problem (13) and the superslow reduced problem (15) obtained from the 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 , denote the fast, slow, and superslow flows at the double singular limits.
We start by showing the singular orbit construction for the -subsystem (Fig. 5), which is a relaxation oscillation that is independent of , , and . The singular orbit (Fig. 5, green trajectory) consists of a superslow excursion along the left branch of -nullcline ( , phase ), a slow jump at the lower fold of -nullcline up to its right branch ( , phase ), a superslow excursion through the active phase ( , phase ), and a slow jump back to its left branch ( , phase ). 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 -space, the -nullcline (red) corresponds to the superslow manifold and the singular orbit overlaps with the singular orbit from the double singular limits due to its independence on . For sufficiently small perturbation , perturbs to the full system trajectory shown by the black curve.
Figure 6(a) illustrates the singular orbit for in -space, together with the superslow manifold (red solid curve : attracting; red dashed curve : repelling). In other panels, the critical manifold (blue surface) and its folds (blue curves) are also plotted. is separated into three sheets by the folds , in which the upper ( ) and lower ( ) branches are stable and the middle branch ( ) is unstable.
Projection of singular periodic orbit (green curve) of system (2) for [(a)–(c)] and (d) to -space. Specifically, panels (a) and (c) show the singular orbits for , and , respectively. Panels (b) and (d) show the singular orbits at the double singular limit . Also shown are the superslow manifold (red curves), the critical manifold (blue surface), and folds of (blue curves). The solid (respectively, dashed) red curves are the branches of that are attracting (respectively, repelling) under the superslow flow. Yellow symbols represent transition points between slow and superslow flow for the oscillation as shown in Fig. 5. The red circle denotes the fast subsystem DHB . 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 and evolves along the lower branch of under the superslow reduced problem (15). After hitting the DHB where the stability of changes, the slow layer problem (13) takes over but with still evolving on a superslow timescale (see Fig. 5, phase ). Thus, the singular orbit during the rest of this phase is a continuum of relaxation oscillations. As the evolution speed of changes from superslow to slow at the yellow circle (beginning of phase ), a few more spikes occur before the slow flow travels to the yellow square on , at which phase begins. After that, the superslow reduced problem takes over until the singular orbit reaches the upper DHB at which becomes unstable. As such, a family of oscillations emerges as we observed during phase . As phase begins at the yellow triangle, several additional spikes occur before the slow flow 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 . Instead, we observe a continuum of spikes throughout phase . This is because the upper becomes unstable for all values as , which will be discussed further in Sec. III B 1. In (b), the singular orbit consists of (triple arrow) that are fast jumps from , (double arrow) which travels along stable branches of or the intersection of and the -nullcline, and (single arrow) that follows the stable branch of .
While the value at the slow/superslow transition (yellow circle or yellow triangle) is uniquely determined, the corresponding values can assume arbitrary positions on the upper or lower sheet of for that fixed . Thus, infinitely many singular orbit segments can be constructed during phases and , 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 during phases and .
Comparing (a) and (b) indicates that to obtain MMO dynamics in the perturbed system, cannot be too small. To explain the MMO dynamics for in Sec. III, we mainly make use of the singular orbit in panel (a) with . This orbit, aside from the upper superslow segment where the stability of differs between the two panels, can be viewed as an 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 jumps and (16), which describes the slow motion along .
For completeness, we also plot the singular orbit for but in Fig. 6(c). Instead of a continuum of spikes, we observe a finite number of spikes during phases and since . In this limit, the fast segment ( , triple arrow) is the same as , whereas the slow segments are perturbations of or from (b). The latter has also been illustrated in Fig. 5.
Finally, the singular orbit for at the double singular limits is shown in Fig. 6(d). There are two major differences between and : First, in contrast to the case where the upper DHB vanishes, it now converges to the upper CDH at the double singular limits in (d). Second, in panel (d) follows the upper stable branch of until reaching a saddle-node bifurcation at the yellow triangle where becomes unstable. This results in a unique construction of singular orbit segment during phase due to its unique starting position. Nonetheless, the construction of singular orbits during phase remains non-unique as discussed in the case. Later in Sec. IV, we show how this singular orbit perturbs to MMOs at for various perturbations.
III. ANALYSIS OF MMOs WHEN gsyn = 4.3 mS/cm2
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 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 MS and MSS
The MMO solution of (2) for from Fig. 1 is projected onto the -space in Fig. 7. The full system equilibrium (black circle) lies on and is unstable. The stability of the superslow manifold changes at the two DHBs (red circles). In particular, as decreases, the upper branch of 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 always lies above the upper branch of the superslow manifold so there is no upper CDH, whereas the lower fold intersects the lower branch of at a CDH singularity (blue diamond in the lower inset). This CDH is a folded focus and will become an FSN1 as as discussed in Sec. II D. Moreover, the nearby HB bifurcation (red circle in the lower inset) is close to this CDH.
Projection of an attracting MMO solution trajectory (black curve) of system (2) for from Fig. 1 to -space. Also shown are parts of the singular orbit from Fig. 6(a) (green curves), the critical manifold (blue surface), folds of (blue curves), and the superslow manifold (red curves). The upper inset shows the upper branch of is always below the upper fold of , 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 intersects the lower fold of . The lower HB bifurcation (red circle) is 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 , the upper superslow segment perturbs to , which displays two SAOs around the stable branch of . These SAOs soon transition to large-amplitude oscillations before crossing the DHB at the red circle to reach the unstable branch of . As phase begins at the green triangle, the slow jump down of brings the trajectory to a region of where there is a nearby stable that attracts the trajectory. From the green star, the trajectory follows on the superslow timescale to the lower fold of (see Fig. 7, lower inset), where it jumps up to on the fast timescale under (7), which corresponds to the onset of spikes in . The spikes persist until reaching the green circle, at which phase begins and the slow jump up of brings the trajectory to the green square, completing a full cycle.
We claim that the only mechanism underlying the MMOs at is the DHB mechanism. This is not surprising given that the upper switches to a continuum of big spikes when the upper DHB vanishes at the singular limit , as shown in Fig. 6(b). To understand why canard dynamics are not involved, we view the trajectory in -space (see Fig. 8). Unlike Fig. 7 where and curves of folded singularity ( ) overlap, the -projection captures the structure of the fold surfaces and lets us examine the position of the solution trajectory relative to . To illustrate how SAOs arise from the DHB mechanism, we look at the projection of the trajectory onto -plane, which includes the periodic orbits born at the upper Hopf bifurcation (see Fig. 9).
The solution from Fig. 7 projected onto ( )-space. Also shown are the projections of the fold surface (blue surface), superslow manifold (red curve), and folded singularities (green and yellow curves). The cyan star denotes 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.
Solutions of (2) for , and various values of and bifurcation diagrams of the slow layer problem (13), projected onto -plane. From (a) to (b) to (c), increases from to to its default value of . In addition to the trajectory (black), the singular orbit (green) and (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 exhibit the same values as the folds of the -nullcline (green circle and green triangle). The lower magenta circle (see the inset) represents the actual fold of (denoted as ) and is not a fold of the -nullcline. Other symbols have the same meanings as in Fig. 7. (d) Real part of the eigenvalues of the upper triangular block of (24), the Jacobian matrix of the slow layer problem (13) along the superslow manifold . The eigenvalues along are real when there are two branches of (blue curves) and complex when there is a single branch of (black curve).
1. Canard mechanism does not contribute to MMOs
In Fig. 8, the solution of (2) for is projected onto the space with two separate fold surfaces (blue) and two branches of folded singularities. In the lower folded singularity curve, there is an (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 . The upper folded singularity curve also has an (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 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 -plane, along with singular and perturbed orbits [see Figs. 9(a)–9(c)]. Both the upper and lower DHBs (red circle) are subcritical. As 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 as it coalesces with a third family of unstable periodic solutions born in the upper DHB.
Figures 9(a)–9(c) illustrate how the singular orbits from Fig. 6(a) perturb to MMOs for various perturbations . With small perturbations [Fig. 9(a), ], the trajectory closely follows the attracting side of (denoted as ), passes over the upper DHB (red circle) to the repelling side of (denoted as ), and undergoes a delay in which the trajectory traces before it jumps away. As the perturbation size increases, the delay is less substantial and the observed small oscillations around become larger and fewer [see Fig. 9(b)]. Panel (c) illustrates the projection of the full system trajectory under the default perturbation . 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 in panel (a), the delay is no longer present and the trajectory jumps away from before reaching the upper DHB. While this might initially suggest that the default perturbation is distant from the singular limit, we show below that this is not the case. Moreover, this specific value of aligns with the perturbation used in Ref. 35, where GSPT analysis has been successfully employed to elucidate the dynamics across various 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 approaches , 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 during phase , 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 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 is significantly shorter than that near .
![An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g013.jpg An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g013.jpg](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11087137/bin/CHAOEH-000034-053119_1-g013.jpg)
Effect of varying (or ) on time traces of solutions of (2) for and other parameters given in Table I. (a) Decreasing preserves MMOs, but the amplitude and the number of the small oscillations alternates between increasing and decreasing with . (b) Increasing leads to a total of five transitions between MMOs and non-MMOs. MMOs exist for [i.e., ] and are completely lost for [i.e., .
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 , the CDH singularity only exists on the lower fold surface of . 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).
2-parameter bifurcation curves for (13) projected onto ( )-space when . (a) The bifurcation curve of the upper HB (Fig. 7, upper red circle), with a horizontal asymptote at . (b) The bifurcation curves of the lower HB (red curve) and the lower fold of (magenta curve), which meet at a Bogdanov–Takens (BT) bifurcation. As (i.e., ), the lower HB converges to the lower CDH (blue).
Relation of (22), (23), and the CDH point on the lower fold for and various values of . (a) (blue star) converges to the (blue diamond) as (equivalently, ). (b) (cyan star) are located far away from the and converge to the leftmost cyan star at the intersection of the folded singularity curve and [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 are summarized by the two-parameter bifurcation diagrams of (13) projected onto -space (see Fig. 10). As decreases (or equivalently, as decreases), the upper Hopf moves to larger values and eventually vanishes for small enough [see Fig. 10(a)]. This explains why the upper singular orbit in Fig. 6(a) for switches to a continuum of spikes in Fig. 6(b) as . On the other hand, since there is a CDH on the lower , the lower Hopf will converge to that CDH as [see Fig. 10(b) and recall Remark II.3]. When increases, the lower Hopf and the lower fold of 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 [also see Figs. 9(a) and 9(c), lower magenta circle].
2. Effect of δ on FSN points
There is no CDH or on the upper fold, hence, we only examine the effect of on FSN singularities on the lower . Figure 11(a) shows as (or equivalently, ), the singularity converges to the lower CDH as demonstrated in our analysis [see (28)]. On the other hand, the 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 close to an and 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 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 and eventually vanish upon crossing the actual fold of [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 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 , i.e., close to a double zero eigenvalue at a BT bifurcation of the 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 along which the Jacobian matrix (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 , excluding the third eigenvalue given by . The eigenvalues are real when there are two branches of curves for and become complex when the curves coalesce and become a single branch. In panel (d), the eigenvalues on the stable lower branch of are initially real and negative. That is, the trajectory approaches the attracting 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 , 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 and vary by changing . Recall that increasing (respectively, decreasing) or slows down (respectively, speeds up) the fast variable , whereas increasing (respectively, decreasing) or speeds up (respectively, slows down) the superslow variable . Other (slow) variables are not affected.
Figure 3(a) summarizes the effects of on MMOs when . 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 . Specifically, we observe the following:
- For fixed at , slowing down the fast variable by increasing from its default value [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.
Effect of varying (or ) on time traces of solutions of (2) for and other parameters given in Table I. (a) Increasing from its default value (equivalently, increasing from ) preserves MMOs. The number of small oscillations increases with . (b) Decreasing from the default value leads to a transition from MMOs ( ) to non-MMOs ( ) to MMOs again ( ) to non-MMOs ( ).
- For fixed at , speeding up the fast variable via decreasing 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 at , slowing down the superslow variable by decreasing from the default value [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 decreases. Additionally, the amplitudes of the small oscillations also display a similar non-monotonic behavior [see Fig. 13(a)].
- For fixed at , speeding up the superslow variable by increasing 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 C1
Increasing slows down the fast variable and hence moves the three timescale (1F, 2S, 1SS) problem closer to (3S, 1SS) separation. As a result, the critical manifold and the folded singularities become less relevant with increased . Nonetheless, this does not affect the existence of MMOs, as they arise from the upper Hopf bifurcation. As is increased, the upper HB moves to smaller values [see Fig. 10(a)]. This change allows the trajectory to travel a longer distance along the stable part of during phase and generate more SAOs, as shown in Fig. 14. To simplify the presentation, we omit singular orbits and focus only on 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 An external file that holds a picture, illustration, etc.
Object name is CHAOEH-000034-053119_1-g014.jpg](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC11087137/bin/CHAOEH-000034-053119_1-g014.jpg)
Solutions of (2) for and various values of and bifurcation diagrams of the slow layer problem (13), projected to -space. From left to right, ( ), ( ), ( ). Increasing moves the upper HB bifurcation (red circle) to lower values and hence increases the number of small oscillations surrounding . Color codings and symbols are the same as in Figs. 9(a)–9(c).
Recall that the small oscillations occurring during phase switch to large spikes upon crossing the inner unstable periodic orbit branch that is born at the upper HB [see Fig. 9(c)]. Increasing causes the upper HB in -space (Fig. 14, red circle) to move to the left and become further away from the maximum of (green square), which remains unchanged. As a result, the trajectory with larger 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 [compare Figs. 9(c) and and14].14]. As is increased to , the trajectory passes over the HB to 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 at the green triangle.
2. Decreasing C1 leads to three MMOs/non-MMOs transitions
Contrary to the preservation of MMOs with increasing , MMOs appear to be sensitive to the decrease of . Specifically, speeding up the fast variable by decreasing leads to transitions from MMOs (e.g., at ) to non-MMOs (e.g., ), back to MMOs at , and then to non-MMOs for .
The initial decrease of [e.g., from to , 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 . Interestingly, a further decrease of to a range of , 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 and 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 grow in amplitude as the trajectory spirals away from . When is reduced to be below , the voltage 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), ].
Projections of the full system solutions and bifurcation diagrams of (13) for and (a) ( ), (b) ( ), (c) ( ), (d) ( ). The upper DHB (red circle) moves to larger values with decreased and eventually vanishes for . Color codings and symbols have the same meaning as in Figs. 9(a)–9(c).
To better understand why the SAOs for grow in amplitude [see Figs. 15(b) and 15(c)], we project the trajectory when onto space (see Fig. 16). The blue triangle denotes a saddle-focus of the subsystem near the green square. During phase 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 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 bifurcation48 and a singular Hopf bifurcation in two-timescale settings.1
The solution (black curve) from Fig. 15(c) projected onto space, together with (red curve). The blue triangle denotes a saddle-focus equilibrium on with maximum , which has one negative real eigenvalue and a pair of complex eigenvalues with positive real parts ( ). A stable manifold branch associated with the saddle-focus at the blue triangle is denoted by the magenta curve.
3. Decreasing ϕ2 preserves MMOs but causes non-monotonic effects on SAOs
Decreasing the perturbation parameter (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 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 decreases from to , 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 , the MMOs exhibit less SAOs with larger amplitudes than the SAOs at [see Fig. 13(a), the third and the fourth row]. The number of SAOs increases and their amplitudes decrease again as decreases from to . This alternating pattern of changes in the number and amplitude of SAOs repeats as continues to decrease [see Fig. 13(a)].
Our analysis reveals that as 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 , 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 at the beginning of phase , 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 during phase . 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 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 causes non-monotonic effects on SAOs.
4. Increasing ϕ2 leads to five MMOs/non-MMOs transitions
When increases, the superslow variable 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 large enough. Indeed, we observe a total of five transitions between MMOs and non-MMOs as increases, and eventually, MMOs are lost for [i.e., . The mechanism driving these MMOs/non-MMOs transitions over the increase of is similar to the mechanism underlying the non-monotonic effects on SAOs when is decreased. Specifically, if no additional big (full) spike before phase is lost with the increase of , there will be fewer SAOs with larger amplitudes or no MMOs as one would naturally anticipate [e.g., when increases from to , see Fig. 13(b), top two rows]. In contrast, if one full spike is lost during the process of increasing , changes to the SAOs will be reversed such that there will be more SAOs with smaller amplitude [e.g., when increases from to , see the top row in Figs. 13(a) and 13(b)]. Eventually, MMOs will be completely lost when for which the HB is no longer relevant.
For a more detailed discussion, we refer the readers to Appendix E.
IV. ANALYSIS OF MMOS WHEN gsyn = 4.4 mS/cm2
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 MS and MSS
The solution of (2) for in Fig. 1 is projected onto the -space, together with the singular orbit from Fig. 6(d) (green curve), critical manifold (blue surface), fold (blue curve), and the superslow manifold (red curves) (see Fig. 17). As the coupling strength increases from to , two new features regarding the upper emerge. First, the stability of the upper changes at a fold point (the green triangle) rather than a DHB. Specifically, as increases, the equilibrium along the upper 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 now intersects the upper branch of 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 close to a folded saddle-node singularity FSN1 (upper blue star) and close to a HB (upper red circle). This is further confirmed in Fig. 18, which shows that the upper and HB point converge to the same CDH point on the upper in the double singular limits .
Projection of a solution trajectory (black curve) of system (2) for from Fig. 1(b) to -space. Also shown are parts of the singular orbit from Fig. 6(d) (green curve), the critical manifold (blue surface), folds of (blue curves), and the superslow manifold (red curves). The solid (respectively, dashed) red curves represent the attracting (respectively, repelling) branches of . Yellow and green symbols mark the transitions between slow and superslow pieces of the oscillations as in Fig. 5. The blue diamonds denote CDH singularities—the intersection of and . 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 close to the folded saddle-node singularity (blue star) and 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.
The upper represents the interaction of canard and DHB mechanisms when . (Left) The upper singularity (blue star) converges to the upper (blue diamond) as (or equivalently ). (Right) As (or ) decreases, the upper HB (red curve) moves to larger values and approaches the upper (blue curve) in the singular limit .
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 , the singular orbit (green) traces 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 . 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 during phase . We omit explaining the solution dynamics of (2) for during other phases as they are similar to those observed when . In particular, the absence of SAOs near the lower CDH point is due to the same mechanism as discussed in Sec. III C.
Time traces of the solutions of (2) (left panels) and their projections onto the corresponding bifurcation diagrams in -space (right panels) for and (a) , (b) , and (c) ( ). Increasing moves the upper DHB (red circle) to the lower and manifests DHB-like features. Color coding and symbols of the bifurcation diagrams on the right panel are the same as in Fig. 17.
Time traces of the solutions of (2) (left panels) and their projections (black curves) onto -space (right panels) for and (a) (i.e., ), (b) (i.e., ), (c) (i.e., ), and (d) (i.e., ). Increasing preserves MMOs with more canard-like feature. Codings and symbols on the right panels are the same as in Fig. 20.
At default parameters, . The emergence of SAOs for 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 , as seen in examples considered in Refs. 26 and 27. To understand why folded node singularities play an important role in the occurrence of SAOs, we view the trajectory and folded singularities projected onto -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 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 , and their relations to the full trajectory.
Projections of the solution from Fig. 17 (black) and the singular 3D funnel volume corresponding to the curve of folded nodes onto - space (top) and - plane (bottom). The 3D funnel volume is bounded between the singular canard surface (magenta surface) and the fold surface . 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.
Projection of the solution trajectory and superslow manifold (red curve) from Fig. 17 onto ( )-plane. The solid (respectively, dashed) cyan curves denote stable (respectively, unstable) periodic orbit branches. Magenta circles represent saddle-node bifurcations of , in which the upper two have the same values as the folds of the -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 -space when . The upper folded singularity curve comprises folded singularities of two types: folded nodes (solid green) and folded saddles (dashed green). The blue star denotes , which occurs away from the upper CDH at the blue diamond. Different from the previous case ( ) where folded singularities did not contribute to the SAO dynamics, the current solution trajectory crosses the upper 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 (see Fig. 20). As discussed in Sec. II B, the funnel for a folded node point represents a two-dimensional trapping region on . The 2D funnels for all folded nodes on the upper fold together form a three-dimensional funnel volume. In the 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 -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 at a regular jump point and there are no SAOs. These observations suggest that the MMOs for 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 -plane. Starting at the green square, the trajectory exhibits SAOs as it follows the upper branch of toward the left. As the trajectory passes through the attracting region of , the oscillation amplitude decays in a typical DHB fashion. Moreover, the orbit experiences a delay along the repelling for an amount of time as it passes through the HB. It is worth noting that after the trajectory enters the unstable part of , 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 switches from a superslow to a slow timescale at the green triangle, as discussed earlier when .
Thus, for , the MMOs for 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 to make the canard mechanism less relevant. This is because increasing slows down the evolution speed of , which in turn drives the three-timescale system (2) closer to (3S, 1SS) splitting. We observed that MMOs persisted for as large as 80 ( ), at which the folded singularities were no longer relevant (see Fig. 22). Second, we increased by raising to drive the system closer to (1F, 3S) splitting, and observed that SAOs persisted for as large as 0.01 ( ), 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 to timescale variations, the interaction of canard and DHB mechanisms due to the existence of the upper CDH when makes MMOs much more robust, as discussed above and illustrated in Fig. 3(b). Specifically, MMOs persist over biologically relevant ranges of and (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 exhibit both canard- and DHB-like features, whereas by tuning , DHB-like features diminish and the canard mechanism dominates.
Below we summarize the effects of increasing the singular perturbations:
- For fixed , increasing (i.e., increasing ) enhances the DHB-like features of the SAOs (see Fig. 22).
- For fixed , increasing (i.e., increasing ) makes the canard mechanism dominant and causes DHB-like features to gradually vanish (see Fig. 23).
1. Increasing C1 makes DHB dominate
Increasing drives the three-timescale system (2) closer to (3S, 1SS) splitting. As a result, the critical manifold and folded singularities including the CDH point become less meaningful and eventually irrelevant for sufficiently large (e.g., ). 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 increases, while the influence of the canard mechanism becomes less significant.
Figure 22 shows the effect of increasing on voltage traces of the full system and the bifurcation diagrams of the fast subsystem (13). As increases from panel (A2) to panel (C2), the upper DHB (red circle) moves away from the upper CDH (blue diamond) to smaller values, whereas the CDH points and slow/superslow timescale transitions denoted by magenta and green symbols all remain unaffected by . As a result, the trajectory with larger begins small oscillations at a larger distance from the Hopf bifurcation point, similar to what we observed in the case of . Furthermore, we have noticed that trajectories for larger exhibit more pronounced DHB-like characteristics, including SAOs with decreasing amplitudes and a more extended travel distance along the unstable branch of . As discussed before, due to a switch of the timescale at the green triangle, there are no oscillations with growing amplitudes as one would expect in a typical DHB fashion.
2. Increasing ϕ2 makes the canard mechanism dominate
Increasing speeds up the superslow variable and hence drives the three-timescale system (2) closer to (1F, 3S) splitting. As a result, the superslow manifold and the DHB points become less relevant and eventually no longer meaningful for sufficiently large (e.g., ). Nonetheless, MMOs continue to persist due to the existence of the canard mechanism. Figure 23 reveals the presence of canard-like features for values across a wide range (from to ). 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 increases, small oscillations tend to pull away from and lose their DHB-like features (i.e., oscillations with decaying amplitude and a delay after passing the HB). Specifically, MMOs with [e.g., Figs. 19 and 23(a)] exhibit both canard and DHB characteristics. When [Fig. 23(b)], there are still some DHB-like features. Further increasing to or [Figs. 23(c) and 23(d)], the trajectories no longer closely follow and the amplitudes of small oscillations become almost constant, which reflects the absence of DHB-like features.
In summary, increasing makes the canard mechanism dominate and the SAOs exhibit fewer DHB-like features. Conversely, decreasing brings the solution and closer together and amplifies the DHB characteristics of the sustained MMOs (see Fig. 23).
V. DISCUSSION
Mixed-mode oscillations (MMOs) are commonly exhibited in dynamical systems that involve multiple timescales. These complex oscillatory dynamics have been observed in various areas of applications and are well studied in two-timescale settings.1,19–24,49,50 In contrast, progress on MMOs in three-timescale problems has been made only in the recent past (see, e.g., Refs. 30, 27, 26, 41, 11, 39, 38, and 40). In this work, we have contributed to the investigation of MMOs in three-timescale settings by considering a four-dimensional model system of coupled Morris–Lecar neurons that exhibit three distinct timescales. We have investigated two types of MMO solutions obtained with different synaptic strengths ( and ). 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 , only the 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 , both and 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 and the superslow manifold . The point where a fold of the critical manifold intersects the superslow manifold is referred to as the canard-delayed-Hopf (CDH) singularity, which naturally arises in problems that involve three different timescales. Reference 26 considered a common scenario in three-timescale systems where the CDH singularity exists and proved the existence of canard solutions near the CDH singularity for sufficiently small and . Moreover, small-amplitude oscillations (SAOs) constantly occur in the vicinity of a CDH singularity.26,27
In this work, we have reported several key findings that have not been previously observed in three-timescale systems. First, although we have identified the same type of CDH singularity as documented in Ref. 26 with its center manifold transverse to 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 did not support SAOs in either the case of or . 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 . 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 (see Fig. 2), increasing causes a shift of the CDH toward lower 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 spikes (data not shown). With higher values, the transition from spikes to plateau or SAOs can happen during phase instead of phase . 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 (see Sec. III). Ultimately, a sufficiently large value of will bring the upper CDH below the minimal of the trajectory, resulting in a complete loss of MMOs. A comprehensive investigation of these complex transitions induced by variations in 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 , 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 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 [see Figs. 15(b), 15(c), and and16].16]. Eventually, SAOs disappeared entirely when became so rapid that it failed to remain in proximity of to generate small oscillations.
TABLE II.
Effects of varying and δ on MMOs.
gsyn = 4.3 | gsyn = 4.4 | |
---|---|---|
Mechanisms | No upper CDH; Only DHB mechanism | Upper CDH exists; Canard and DHB mechanisms coexist |
Increasing C1 (or ) | MMOs are preserved with more DHB characteristics. | MMOs are preserved with more DHB characteristics. |
Decreasing C1 (or ) | Two MMOs/non-MMOs transitions are observed before a complete loss of MMOs. | MMOs are preserved and organized by both mechanisms. |
Increasing ϕ2 (or δ) | Four MMOs/non-MMOs transitions are observed before MMOs are entirely lost. | MMOs are preserved with more canard-like and less DHB features. |
Decreasing ϕ2 (or δ) | MMOs are preserved, but a non-monotonic effect on the small oscillations is observed. | MMOs are preserved with both DHB and canard characteristics. |
As one would expect, MMOs when 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 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 , an upper CDH occurs near the SAOs at . 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 and 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 , MMOs at 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 [Fig. 10(a)], decreasing at 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 singular limit. Hence, MMOs at 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 due to the presence of a nearby upper CDH, leads to significantly enhanced robustness against timescale variations compared with , 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 is decreased when , unlike what we observed in . 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 with reference timescale . This transforms (2) to the dimensionless system (4) given in the Introduction, namely,
where , , ,
and
From system (A1), we can see that the voltage evolves on a fast timescale , evolve on a slow timescale , whereas evolves on a superslow timescale .
APPENDIX B: FOLDED SADDLE-NODE (FSN) SINGULARITIES OF (18)
To derive the conditions for FSN singularities, we note that the eigenvalues of the folded singularities satisfy the following algebraic equation:
where and denote the trace and determinant of the Jacobian matrix of the desingularized system (18), respectively. As discussed before, along the folded singularity curve . This also directly follows from the following calculations:
where the entries in denote partial derivatives. The last equality in the above equations holds because [see (18)], where and . Hence, the remaining two nontrivial eigenvalues and satisfy
It follows that an is given by
Plugging in functions , , and from (18), we can rewrite the above condition as Eq. (21) in Sec. II B, namely,
where and .
Next, we prove there are two ways an can occur. In case 1, suppose . An can occur when , where . Note that is a folded singularity point (20), that is, . This implies , where . Thus, the first way that an can occur is described as the following:
which is Eq. (22) in Sec. II B. This implies that an point becomes a CDH when .
In case 2, suppose . An can occur when where . As a result, the second way that an can occur is defined as
which is Eq. (23) in Sec. II B. This implies that an is close to the intersection of and .
APPENDIX C: CDH AT DOUBLE SINGULAR LIMIT EXHIBITS TWO LINEARLY INDEPENDENT CRITICAL EIGENVECTORS
Recall the Jacobian matrix at an , 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 of the desingularized system at the double singular limit ( is given by
in which and . At a CDH, . It follows that . Thus, the Jacobian matrix at a CDH singularity becomes
which has a nullity of 2, implying the existence of two linearly independent critical eigenvectors associated with zero eigenvalues. Moreover, the center subspace at the CDH, given by the plane [see pink planes in Figs. 4(c) and 4(d)], is transverse to the fold surface of the critical manifold. This is because the normal vector of at the CDH, which is given by where is not parallel to the normal vector of the center subspace given by . Thus, the CDH singularity (an at the double singular limit) considered in system (2) is a folded saddle-node singularity, with the center manifold of the 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 (denoted by ) is always tangent to . Through direct computations, we obtain with
The dot product of and the unit normal vector of at the CDH is given by
It follows directly that is tangent to as expected.
APPENDIX D: DECREASING WHEN CAUSES NON-MONOTONIC EFFECTS ON SAOS
In this subsection, we explain the mechanism underlying the non-monotonic effects of decreasing 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 , the number of SAOs will decrease and their amplitudes will increase (e.g., when decreases from to ); (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 decreases from to ). Below, we examine the two opposite effects using the projections of the corresponding solutions and bifurcation diagrams onto space (see Fig. 24).
FIG. 24.
Projections of solutions of (2) with , , and different values onto -plane. (a) Projections of trajectories with (black) and (dark yellow). (b) Projections of trajectories with (black) and (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 (black) is smaller than SAOs with (dark yellow). There are also more black SAOs when [see Fig. 13(a)], which, however, is not obvious in the projection as the black SAOs are hardly visible. Compared with the black solution, the yellow one has a slower evolving rate of , which is slaved to , during phase . Thus, one extra full spike occurs for the yellow trajectory compared to the black trajectory during the jump at phase . As a result, the black trajectory, after its last big spike, approaches the black square at its maximum along . In contrast, the yellow trajectory, due to the extra full spike that occurs during the jump, approaches the maximum from the outside of the periodic orbit branch. Consequently, the yellow trajectory that is being pushed further away from exhibits fewer small oscillations with larger amplitudes compared with the black trajectory, which remains close to .
When is reduced from to , the solution trajectory changes to the green curve in Fig. 24(b). Slowing down of 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 , similar to the case when . 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 WHEN 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 in the case of (see Fig. 13 in Sec. III D), which is similar to the mechanism that causes the non-monotonic effects on SAOs when is decreased (see Appendix D).
During the increase of which speeds up , the big spikes produced during phase will gradually decrease. If one big (full) spike before phase is lost during this process, there will be more SAOs with smaller amplitudes [e.g., when increases from to in Fig. 25(a)] or the earlier lost MMOs will reappear [e.g., when increases from 0.0016 to 0.0018 in Fig. 25(c)]. As increases from to in Fig. 25(a), the number of full spikes during phases and decreases by one. Consequently, all three green full spikes occur within phase , while there is a large black spike during phase . As explained in the previous subsection, this causes the green trajectory to approach the maximum along 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 from to leads to a decrease in the number of big spikes that occur before phase , more SAOs are observed and the transition from non-MMOs to MMOs occurs [see Fig. 25(c)].
FIG. 25.
In contrast, increasing from to does not change the number of large spikes before phase [compare the green and dark yellow trajectories in Figs. 25(a) and 25(b)]. However, speeding up of 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 , 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.