Complex oscillatory patterns near singular Hopf bifurcation in a two time-scale ecosystem

TL;DR

This paper analyzes a two-time-scale ecological predator-prey model exhibiting complex oscillatory behaviors, including mixed mode oscillations and chaos, near singular Hopf bifurcations, with implications for understanding ecosystem regime shifts.

Contribution

It introduces a detailed bifurcation analysis of a predator-prey model with explicit predator interference, revealing novel oscillatory phenomena near singular Hopf bifurcations.

Findings

Presence of mixed mode oscillations and relaxation oscillations.

Long transient chaotic dynamics near bifurcation points.

Oscillations driven by variations in predator intraspecific competition.

Abstract

We consider an ecological model consisting of two species of predators competing for their common prey with explicit interference competition. With a proper rescaling, the model is portrayed as a singularly perturbed system with one-fast (prey dynamics) and two-slow variables (dynamics of the predators). The model exhibits variety of rich and interesting dynamics, including, but not limited to mixed mode oscillations (MMOs), featuring concatenation of small and large amplitude oscillations, relaxation oscillations and bistability between a semi-trivial equilibrium state and a coexistence oscillatory state. Existence of co-dimenison two bifurcations such as fold-Hopf and generalized Hopf bifurcations make the system further intriguing. More interestingly, in a neighborhood of {\emph{singular Hopf}} bifurcation, long lasting transient dynamics in form of chaotic MMOs or relaxation oscillations are observed as the system approaches the periodic attractor born out of supercritical Hopf bifurcation or a semi-trivial equilibrium state respectively. The transient dynamics could persist for hundreds or thousands of generations before the ecosystem experiences a regime shift. The time series of population cycles with different types of irregular oscillations arising in this model stem from a biological realistic feature, namely, by the variation in the intraspecific competition amongst the predators. To explain these oscillations, we use bifurcation analysis and methods from {\emph{geometric singular perturbation theory}}.

Equations

$$\displaystyle\left\{\begin{array}[]{ll}\frac{dX}{dT}&=rX\left(1-\frac{X}{K}\right)-\frac{p_{1}XY}{H_{1}+X}-\frac{p_{2}XZ}{H_{2}+X}\\ \frac{dY}{dT}&=\frac{b_{1}p_{1}XY}{H_{1}+X}-d_{1}Y-a_{12}YZ\\ \frac{dZ}{dT}&=\frac{b_{2}p_{2}XZ}{H_{2}+X}-d_{2}Z-a_{21}YZ-mZ^{2}\end{array}\right.$$

$$X(0)=\tilde{X}\ge 0,Y(0)=\tilde{Y}\ge 0,Z(0)=\tilde{Z}\ge 0,$$

$$t$$

$${\beta }_2$$

$$Y_0=\frac{rK}{p_1},Z_0=\frac{rK}{p_2},$$

$$\displaystyle\left\{\begin{array}[]{ll}{x^{\prime}}&=x\left(1-x-\frac{y}{\beta_{1}+x}-\frac{z}{\beta_{2}+x}\right)\\ {y^{\prime}}&=\zeta_{1}y\left(\frac{x}{\beta_{1}+x}-c-\alpha_{12}z\right)\\ {z^{\prime}}&=\zeta_{2}z\left(\frac{x}{\beta_{2}+x}-d-\alpha_{21}y-hz\right),\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}{x^{\prime}}&=x\left(1-x-\frac{y}{\beta_{1}+x}-\frac{z}{\beta_{2}+x}\right):=x\phi(x,y,z)\\ {y^{\prime}}&=\zeta y\left(\frac{x}{\beta_{1}+x}-c-\alpha_{12}z\right):=\zeta y\chi(x,z)\\ {z^{\prime}}&=\zeta z\left(\frac{x}{\beta_{2}+x}-d-\alpha_{21}y-hz\right):=\zeta z\psi(x,y,z),\end{array}\right.$$

$$\displaystyle\left\{\begin{array}[]{ll}\zeta\dot{x}&=x\phi(x,y,z)\\ \dot{y}&=y\chi(x,z)\\ \dot{z}&=z\psi(x,y,z),\end{array}\right.$$

$$\zeta =0.01,{\beta }_2=0.35,c=0.4,d=0.21,{\alpha }_{12}=0.5,{\alpha }_{21}=0.1,$$

$$\zeta =0.01,{\beta }_2=0.6,c=0.4,d=0.21,{\alpha }_{12}=0.5,{\alpha }_{21}=0.1,$$

$$\left\{\begin{array}[]{r c l}{x^{\prime}}&=&x\phi(x,y,z)\\ {y^{\prime}}&=&0\\ {z^{\prime}}&=&0.\\ \end{array}\right.$$

$$\left\{\begin{array}[]{r c l}0&=&x\phi(x,y,z)\\ \dot{y}&=&y\chi(x,z)\\ \dot{z}&=&z\psi(x,y,z).\\ \end{array}\right.$$

$$\mathcal{M}=\left\{(x,y,z):x=0\text{leavevmode}\mathrm{or}\text{leavevmode}\varphi (x,y,z)=0\right\}:=T\cup S,$$

$$\displaystyle\left\{\begin{array}[]{ll}x&=0\\ \dot{y}&=y\chi(0,z)\\ \dot{z}&=z\psi(0,y,z).\end{array}\right.$$

$$\displaystyle\begin{pmatrix}-\phi_{x}\dot{x}\\ \dot{y}\end{pmatrix}=\begin{pmatrix}\phi_{y}y\chi+\phi_{z}z\psi\\ y\chi\end{pmatrix}\bigg{|}_{z=\theta(x,y)}.$$

$$\displaystyle\begin{pmatrix}\dot{x}\\ \dot{y}\end{pmatrix}=\begin{pmatrix}\phi_{y}y\chi+\phi_{z}z\psi\\ -\phi_{x}y\chi\end{pmatrix}\bigg{|}_{z=\theta(x,y)},$$

$${\int }_0^{{\tau }_0}\varphi (0,y(s),z(s))ds=0,$$

Full text

Key Words. Mixed-mode oscillations, singular Hopf bifurcation, long transients, basins of attraction, fold-Hopf bifurcation, chaotic dynamics.

1. Introduction

Variability in population abundance is one of the most ubiquitous phenomena in natural populations. These dynamics typically involve rapid transitions from one oscillatory state to another, featuring temporal patterns of small fluctuations interspersed with large amplitude oscillations [27, 37], as observed in population cycles of Pacific sardine and northern anchovy [3], forest insects such as larch bud moth, gypsy moth and cankerworms [2, 23], agricultural pests such as desert locusts [11, 53] and so forth. Using singular perturbation theory, regular cycles of population outbreaks and collapses are modeled by relaxation oscillation cycles, commonly known as boom and bust cycles [30, 36, 41]. However, since events of population outbreaks are not regular, vary in intensities and generally not periodic in time, they are perhaps not best represented by relaxation oscillations. A more realistic representation of these cycles can be mixed-mode oscillations (MMOs), which are complex oscillatory patterns consisting of one or more small amplitude oscillations (SAOs) followed by large amplitude oscillations (LAOs) which are of relaxation type [12]. MMOs have been widely studied in various physical systems, including, but not limited to computational neuroscience models, autocatalytic and electrochemical reactions (see [5, 8, 22, 24, 29, 31, 42] and the references therein for examples of mathematical models exhibiting MMOs). In ecology, there have been relatively few studies on such complex oscillatory patterns [6, 27, 37, 44, 45, 46].

Traditionally, theoretical ecological studies have focussed on asymptotic analysis to study population dynamics. However, there has been an increasing recognition of importance of understanding the properties of transient dynamics in complex ecosystems [19, 34]. One of the major issues in ecology is predicting sudden changes in population abundance, commonly referred to as regime shifts [49, 50], as they often result in population collapse and extinctions. Long term transient dynamics on ecological timescales can offer novel perspectives in underlying mechanisms behind regime shifts in an apparently constant environment. The present work unifies the two phenomena that appear in natural populations, namely MMOs featuring oscillations on different timescales, and their occurrence as long transient dynamics over ecological timescales.

In this paper, we study the dynamics between three interacting species, namely two predators competing for their common prey with explicit Lotka-Volterra type interference competition between the predators. Taking the ratios between birth rates of the predators to the prey extremely small, separation of time scales is introduced into the system as a singular parameter. This then brings the model in the framework of singularly perturbed system of equations with one-fast and two-slow variables, whereby the prey exhibits fast dynamics and the predators have slow dynamics. Similar assumptions have been made in many ecological models, including but not limited to forest pest models [6], tri-trophic food chain models [9, 10], age-structured predator-prey models with dormancy of predators [27] and the well-known Rosenzweig-MacArthur predator-prey model [38]. Though the model studied in this paper is a generic ecological model, it can give rise to a host of new and interesting patterns of behavior. These dynamics seem to represent natural population cycles. A primary goal of this paper is to study various kinds of complex oscillatory dynamics that emerges in this model and use tools from geometric singular perturbation theory and bifurcation analysis to explain them.

Treating the intraspecific competition amongst one of the species of predators as the main control parameter and the predation efficiency of the other species of predators as the secondary parameter ($\beta_{1}$), the local and global bifurcation structures in the system are explored. We find different dynamical regimes on which interesting ecological phenomena occur, and compare the bifurcation structures of the model by considering different values of the predation efficiency of the other predator (denoted by $\beta_{2}$). Bifurcations such as saddle-node, Hopf, period-doubling, torus, fold-Hopf and generalized Hopf occur in the system and contribute to unfolding different interesting phenomena exhibited by the model.

One of the interesting dynamics observed in this model is the existence of long term transients in form of chaotic MMOs consisting of very long epochs of SAOs near the singular-Hopf point, also referred to as folded saddle-node of type II (FSN II) bifurcation point [25, 26]. The transient MMO dynamics persists for thousands of generations and eventually asymptotes to the small amplitude periodic attractor born out of a supercritical Hopf bifurcation. The stable periodic attractor possesses a large basin of attraction and interferes with the global return mechanism that is responsible for the relaxation dynamics of the MMOs. The SAOs associated with the chaotic transients are organized by a slow passage near special points in phase space such as folded node singularities [26, 52] and are further influenced by the local vector field around a saddle-focus equilibrium which lies in a close vicinity to the folded node singularity. The existence of long lasting transient dynamics reflects that the dynamics on an ecological timescale can be totally different than asymptotic dynamics. The occurrence of a supercritical Hopf bifurcation, which is regarded as a safe bifurcation in context of ecological regime shifts [49, 50] can be misleading on an ecological timescale.

In another parameter regime, where the coexistence equilibrium undergoes a subcritical Hopf bifurcation, transient dynamics in form of relaxation oscillations are observed before the system approaches a boundary equilibrium state, where one of the predators becomes extinct. In a nearby neighborhood of the subcritical Hopf bifurcation, a torus bifurcation occurs, giving birth to a family of stable periodic orbits. In this regime, the system exhibits a bistable behavior, namely depending on the initial densities, the model could either result into coexistence of species in an oscillatory state or lead to a competitive exclusion of one of the species of predators.

Long term transients in form of MMOs or relaxation oscillations shed light on the underlying inherent uncertainties associated with population dynamics. These can be ecologically alarming, as abrupt large population fluctuations may lead to outbreaks of harmful pests such as desert locusts or forest insects or bring a species perilously close to extinction under a seemingly constant environment. In fact, such long lasting transient dynamics can provide an alternate explanation of a sudden transition from one state to another in an ecosystem which may show an apparently stable dynamics or long statis [19]. In this model, the local dynamics near the equilibrium can be used to identify early warning signs of a large fluctuation which has been explored in a companion paper [47].

In multiple timescale systems with one fast and two or more slow variables, the two main mechanisms attributed to the generation of MMOs are the generalized canard phenomenon [7, 52] and a singular Hopf bifurcation [12, 15]. In [15], the differences in the characteristics of SAOs in an MMO orbit due to singular Hopf bifurcation and those created due to folded node singularity is discussed: namely, the small oscillations associated with singular Hopf start with low amplitude and then grow before entering into relaxation regime, whereas the SAOs associated with folded nodes decrease and then increase in amplitudes. However in this model, MMOs with both characteristics are seen near FSN II bifurcation, where a folded node lies in the vicinity of the equilibrium. Similar dynamics have been previously observed in neuroscience models [8, 42]. More interestingly in an another parameter regime, MMOs with another unique characteristic is observed: namely, the durations of the episodes of small oscillations between two large oscillations vary by orders of magnitude, coupled with large variations in the amplitudes of the small oscillations. A major focus on this paper is to study different MMO patterns that arise in this model near FSN II bifurcation. We use the presence of different timescales in the system to analyze the mechanism for the complex oscillations. The analyses performed at the singular limit are used to explain the behaviors seen away from the limit.

The paper is organized as follows. In Section 2, we introduce and properly scale the model. The assumptions and physical significance of each parameter are also discussed. Section 3 focusses on examples of different MMO patterns exhibited by the model either as asymptotic dynamics or as long term transients. In parameter regimes, where MMOs occur as prolonged transients, “transient basins of attraction” are computed and compared. Relaxation oscillations as long term transients in a different parameter regime are also studied in Section 3. Section 4 focusses on bifurcation analysis of the model, where interesting dynamics near the folded saddle-node singularity of type II are explored. A geometric singular perturbation approach and relevant background review is given in Section 5, which provides us a framework of local mechanisms that are responsible for MMOs. The numerical results are explained using this framework. Finally, we discuss our results and summarize our conclusions in the last section of the paper.

2. The Model

The model studied in this paper is an extension of the model considered in [43] with additional direct interactions between the predators. The model reads as follows:

[TABLE]

under the initial conditions

[TABLE]

where $X$ represents the population density of the prey and $Y$, $Z$ represent the densities of the two species of predators. The parameters $r$ and $K$ represent the intrinsic growth rate and the carrying capacity of the prey, $p_{1}$ is the maximum per-capita predation rate of $Y$, $H_{1}$ is the semi-saturation constant which represents the prey density at which $Y$ reaches half of its maximum predation rate ($p_{1}/2$), $b_{1}$ and $d_{1}$ are the birth-to-consumption ratio and per-capita natural death rate of $Y$ respectively. The parameter $a_{12}$ is the rate of adverse effect of $Z$ on $Y$, which we will refer to as the interspecific competition rate. The other parameters $p_{2},b_{2},d_{2},H_{2},a_{21}$ are defined analogously for $Z$. We assume that $Z$ is territorial and/or is a social predator, and therefore experiences more competition for space and resources in comparison to $Y$ (which is assumed to be non-territorial or is a solitary predator). The term $mZ^{2}$ in the $Z$ equation accounts for this intraspecific competition, which may include lethal fighting and cannibalism, and measures the density dependent mortality rate in the class of $Z$. Examples of species governed by (4) may include small mammals such as rabbits/rodents preyed upon by raptors such as eagles ($Y$) and larger mammals such as wolves ($Z$) or planktonic and small nektonic species such as fish and squids ($X$) preyed upon by sharks ($Y$) and dolphins ($Z$).

The density dependent mortality term in system (4) is the defining term that distinguishes the possibility of existence of a positive equilibrium state from competitive exclusion of species. In case of two competitors competing for a common species, the stronger competitor can outcompete the other; a phenomenon commonly referred to as the principle of competitive exclusion [18]. Such dynamics are known to occur in Rosenzweig-MacArthur and Armstrong and McGhee predator-prey models [21, 32, 36]. In these models, the only coexistence state that the systems admit are in the form of large oscillation cycles [30]. However, a density dependent mortality term can prevent competitive exclusion and ensure long-term survival of the species besides guaranteeing the existence of its competitor [43]. Thus, the term $mZ^{2}$ is ecologically significant and we will treat the coefficient of intraspecific competition as the key parameter in this work.

With the following change of variables and parameters:

[TABLE]

where

[TABLE]

system $(\ref{maineq})$ takes the following dimensionless form:

[TABLE]

where the primes denote differentiation with respect to the time variable $t$. Similar scaling variables were considered in a three-trophic food chain model by Deng [9]. The quantity $Y_{0}$ will be interpreted as the maximum predation capacity of $Y$ [9]. More precisely, at this capacity, $Y$ will consume $rK$ number of prey (which measures the reproduction capacity of the prey population) at the maximum predation rate $p_{1}$. Similarly, $Z_{0}$ will be considered as the maximum predation capacity of $Z$. We will assume the following conditions on the parameters:

(A) The maximum per capita birth rate of the prey is much higher than the per capita birth rates of the predators, i.e. $b_{1}p_{1}<<r$ and $b_{2}p_{2}<<r$, thus yielding $0<\zeta_{1},\zeta_{2}<<1$. This is usually observed in many ecosystems [40]. Examples include insects and their avian predators, rodents and their aerial or ground-based predators and so forth. For simplicity, we will assume in our model that $\zeta_{1}=\zeta_{2}=\zeta$ (say). Similar dynamics are obtained when $\zeta_{1}\neq\zeta_{2}$, but have the same order.

(B) The parameters $c$ and $d$ satisfy the inequality $0<c,\ d<1$, which implies that the growth rates of the predators are greater than their death rates. This is a default assumption otherwise the predators would die out faster than they could reproduce even at their maximum reproduction rate.

(C) The parameters $\beta_{1}$ and $\beta_{2}$ are dimensionless semi-saturation constants measured against the prey’s carrying capacity. We will assume that both predators are efficient, and hence they will reach the half of their maximum predation rates before the prey population reaches its carrying capacity, thus yielding $0<\beta_{1},\beta_{2}<1$.

(D) The parameters $\alpha_{12}$ and $\alpha_{21}$ are dimensionless interspecific competition coefficients measuring the interference effect of $Z$ on $Y$ and of $Y$ on $Z$ respectively. We will assume that the maximum interference effect of one species on the other does not exceed the intrinsic growth rate of the other species, yielding $0<\alpha_{12},\alpha_{21}<1$.

Under the assumptions (A)-(D), system $(\ref{nondim1})$ transforms to a singular perturbed system of equations with two time scales, where the prey exhibits fast dynamics and the predators exhibit slow dynamics.

We rewrite system $(\ref{nondim1})$ as

[TABLE]

where $\phi=0$, $\chi=0$, and $\psi=0$ are the nontrivial $x$, $y$, and $z$-nullclines respectively. On rescaling $t$ by $\zeta$ and letting $s=\zeta t$, system $(\ref{nondim2})$ can be reformulated as

[TABLE]

where the overdot denotes differentiation with respect to the variable $s$. The variables $t$ and $s$ are referred to as the fast and slow time variables respectively. The parameter $\zeta$ can be regarded as the separation of time scales.

Throughout the paper, we will fix the parameter values to

[TABLE]

or to

[TABLE]

and treat the intraspecific competition coefficient $h$ as the primary control parameter, and the predation efficiency $\beta_{1}$ as the secondary parameter. The parameter values chosen here are for illustrative purposes to demonstrate the complex oscillatory patterns that could occur in a complex ecosystem as shown in figure 1.

System (17) has $8$ free parameters which make the qualitative analysis very challenging. In the absence of interference competition, i.e. $\alpha_{12}=\alpha_{21}=0$, there exists a unique coexistence equilibrium in a suitable parameter regime. MMO orbits and chaotic dynamics exist as a control parameter is varied. This case has been studied in details in [45]. In this paper, we investigate the dynamics in a neighborhood of the singular Hopf point for $\alpha_{12},\alpha_{21}\neq 0$. To this end, we take all the parameter values to be nonzero such that the intersection of the non-trivial nullclines $\phi=0,\chi=0$ and $\psi=0$ produces equilibria that lie in the positive octant. These equilibria will be referred to as coexistence or non-trivial equilibria and will be denoted by $E_{i}^{*}$.

3. Complex oscillatory dynamics in system (17)

3.1. MMOs as asymptotic dynamics

Fixing $\beta_{1}=0.25$ and the other parameter values as in (18), we consider the bifurcations that arise by varying the intraspecific competition $h$. Using AUTO, a one-parameter bifurcation diagram was calculated as shown in figure (2), where $h$ is the continuation parameter and $L_{2}$ norm refers to the standard Euclidean norm. We note from the bifurcation diagram in figure 2 that the coexistence equilibrium point $E_{1}^{*}$ exists for $h>0.7425$. In the singular limit of system (17), FSN II bifurcation occurs at $h\approx 0.7785$ (also see figure 17), where $E_{1}^{*}$ is still an attractor of the full system. A supercritical Hopf bifurcation occurs at $h\approx 0.7803$, giving birth to a family of stable periodic orbits $\Gamma_{h}$. The equilibrium $E_{1}^{*}$ is now a saddle-focus with one negative and two complex (with positive real parts) eigenvalues. The periodic orbits $\Gamma_{h}$ grow in size and undergo a cascade of period-doubling bifurcations for $h\in(0.797,0.798)$, resulting into small amplitude chaotic invariant sets.

Past the regime of period doubling bifurcations of $\Gamma_{h}$, chaotic dynamics featuring small amplitude and large amplitude oscillations are observed. On further increasing $h$, a sequence of period doubling and saddle-node bifurcations of limit cycles produce periodic MMO orbits predominantly with signatures $1^{s}$ where $s\in\mathbb{N}$ can be very large, and other complicated signatures, which are usually a mix of $1^{i}$ and $1^{i+1}$ as shown in figure 4. The SAOs in such an MMO orbit are organized by the rotational properties of the weak canard as well as by the unstable manifold of $E_{1}^{*}$. We will discuss the local mechanisms generating these SAOs in Section 5. Relaxation oscillations are observed for $h>0.98$. A detailed view of the bifurcation structure of the MMOs is given in figure 3.

System (17) exhibits MMOs for different combinations of parameter values. Figure 5 represents an MMO orbit corresponding to $\beta_{1}=0.1977$ and $h=0.633$ and other parameter values as in (18). These parameter values are chosen from a regime near a co-dimension 2 bifurcation point (see figure 10). Note that the SAOs associated with the MMO orbit in figure 4 have much smaller amplitudes than the ones seen in figure 5. The small oscillations start with low amplitude, decrease and then increase before entering into the relaxation regime. In contrast, the small oscillations of the MMO orbit in figure 5 start with moderate amplitude and then grow slowly with increasing amplitude. MMOs with former characteristics have been attributed to canards at folded nodes, while the latter characteristics have been attributed to singular Hopf points in [15]. However, the presence of a folded node near the singular Hopf point in this system (see figure 27) may allow MMOs with only increasing-in-amplitude oscillations but not decreasing-in-amplitude, and hence we see both types of features in this model. Similar characteristics have been observed in other models in the literature such as in [42].

3.2. MMOs as prolonged transients

Considering $\beta_{1}=0.25$ and other parameter values as in (18), we noted that a supercritical Hopf bifurcation of the coexistence equilibrium occurs at $h\approx 0.7803$, giving birth to a stable limit cycle $\Gamma_{h}$.

Interestingly, in the parameter regime where $\Gamma_{h}$ is stable, i.e. for $0.7803<h<0.797$, long term transients in form of chaotic MMOs are observed before the system approaches its asymptotic state. The duration of the transient depends sensitively on the initial values of the state variables, and may last for a significantly long amount of time for certain initial conditions. Also, with the same initial condition, chaotic transients with a correspondingly distinct duration are observed as the control parameter is varied as shown in figure 6.

Figure 7 demonstrates the transition from chaotic MMOs to the asymptotic dynamics, i.e. the Hopf limit cycle, for a fixed set of parameter values. Such complex oscillatory patterns in population densities persist for thousands of generation before the system undergoes a sudden change to a completely different state. A “transient basin” can be thus defined as the set of initial conditions which produce chaotic transients whose duration is within a given interval $[0,N]$ [51]. Numerical approximations of transient basins for a fixed value of $h$ are plotted in figure 8 by choosing equally spaced points on a grid of size $[0.08,0.115]\times[0.325,0.36]$ on the plane $x=0.3428$ and integrating trajectories starting at each point on this grid. Points on the grid which lead to transients of duration less than $N$ are plotted in blue, and are referred to as “basins of short transients”. On the other hand, points in red will be referred to as “basins of long transients”. The transient basins are densely intermixed as can be seen in figure 8.

3.3. Relaxation oscillations as long term transients

Another interesting feature observed in system (17) is the existence of relaxation oscillations as prolonged transients before the system approaches the boundary equilibrium state $E_{xz}$. This occurs when one of the predators, namely $z$ is assumed to be more efficient than $y$. To this end, we let $\beta_{1}=0.1923$ and consider parameter values as in (19).

A subcritical Hopf bifurcation of the coexistence equilibrium occurs at $h=0.2646$, followed by a supercritical torus bifurcation at $h=0.2648$. For $h>0.2648$, the semi-trivial equilibrium state $E_{xz}$ is locally asymptotically stable and coexists with the periodic orbit born out of torus bifurcation. The periodic orbit has a relatively smaller basin of attraction and the system can either approach the periodic orbit or $E_{xz}$. Transient dynamics in form of relaxation oscillations are observed as the system approaches $E_{xz}$ as shown in figure 9. These types of transient dynamics which may apparently seem sustainable on ecologically relevant timescales could result into species extinction, a property that has been studied in many ecosystems [19].

4. Bifurcation analysis

In this section, we explore the bifurcation structure of system (17) and find other interesting dynamics that are generated by this generic model. We begin by treating $h$ as the primary bifurcation parameter and $\beta_{1}$, the predation efficiency of $y$, as the secondary parameter, keeping the other parameter values fixed at (18) or (19). Besides studying the effect of $h$ and $\beta_{1}$ on the system, we will also explore the role of $\beta_{2}$. The two different values of $\beta_{2}$ considered in (18) and (19) will elucidate the influence of predation efficiency of $z$ on the system. We will compute two-parameter bifurcation diagrams in $(h,\beta_{1})$ and then compare the bifurcation structures of system (17) for the two sets of parameters to see how the structure changes with the predation efficiency of $z$.

4.1. A two-parameter bifurcation of system (17) with similar predation efficiencies.

Keeping all the parameter values fixed as in (18), we consider a two-parameter bifurcation in $(h,\beta_{1})$ parameter space. The bifurcation diagram is shown in figure 10. The equilibrium point $E_{xz}$, lying on the $xz$-plane, undergoes a supercritical Hopf bifurcation along the vertical Hopf line $H$ (shown in blue), and the coexistence (non-trivial) equilibrium undergoes a supercritical Hopf bifurcation along the Hopf curve, also denoted by $H$. The FSN II curve and the Hopf curve lie very close to each other. A saddle-node bifurcation of the non-trivial equilibria of (17) occurs along the SN curve (shown in red). The equilibria lying on the SN curve may or may not be biologically feasible (as they could have one or more negative components). In between the SN curve and the branch curve (BR), system (17) could have exactly two or no nontrivial equilibria. A unique nontrivial equilibrium of (17) exists to the right of BR. In the region lying in between BR and the Hopf curve, the nontrivial equilibrium is locally asymptotically stable. To the right of the Hopf curve, oscillatory dynamics such as MMOs and relaxation oscillations appear. The PD curve, corresponding to the period-doubling bifurcation of the periodic orbit born out of the Hopf bifurcation demarcates MMOs from relaxation oscillations in the parameter space. The branch curve (BR) meets with the Hopf curve at the fold-Hopf bifurcation point (FH) which occurs at $(\beta_{1},h)=(0.197019,0.60803)$. Near the fold-Hopf bifurcation, which is a co-dimension two bifurcation [16, 28], other bifurcations may occur including saddle-node bifurcations of periodic orbits and bifurcation of Shil’nikov homoclinic orbit to a saddle-focus. The latter type of bifurcation can serve as an organizing center for the MMOs near the fold-Hopf point, and will be a subject for future investigation.

4.1.1. Dynamics near the fold-Hopf point

In a vicinity of the FH point, above the Hopf curve $H$ and below the branch curve BR in figure 10, system (17) possesses two coexistence equilibria $E_{1}^{*}$ and $E_{2}^{*}$. The equilibrium $E_{1}^{*}$ is a stable focus-node and $E_{2}^{*}$ is a saddle-focus with two-dimensional stable manifold and one-dimensional unstable manifold. As the control parameter $h$ is increased, $E_{1}^{*}$ undergoes a supercritical Hopf bifurcation and thereafter becomes a saddle-focus with one-dimensional stable manifold, $W^{s}(E_{1}^{*})$, and two-dimensional unstable manifold, $W^{u}(E_{1}^{*})$. The equilibrium $E_{2}^{*}$ no longer lies in the first octant and the system asymptotes to a periodic orbit that emerges from the Hopf bifurcation of $E_{1}^{*}$ (see figure 11). On further increasing $h$, MMOs (both aperiodic and periodic orbits) appear (see figure 5). The SAOs associated with the MMOs are canard induced. Such dynamics will be discussed in the next section. The amplitude of the SAOs depend sensitively on $\beta_{1}$, which will be also explored in this paper.

4.1.2. Horizontal bifurcations

Treating the two predators to be equally efficient, i.e. by letting $\beta_{1}=0.35$ and the other parameter values as in (18), a similar bifurcation diagram as in figure 2 is obtained. A supercritical Hopf bifurcation of the coexistence equilibrium occurs at $h\approx 1.5$. MMO orbits with signature $1^{s}$, where $s\leq 5$ are more commonly seen. An MMO orbit of signature $1^{5}$ is shown in figure 1.

Figure 12 illustrates the sequence of MMOs that occurs in this case. The bifurcation diagram demonstrates how the number of small oscillations between two large oscillations vary with the intraspecific competition. Since the return time between successive population outbreaks (to a large extent) can be determined by the number of the SAOs in an MMO orbit, and therefore by the epoch of the SAOs, we note that the frequency of extreme events of outbreaks increases with the intraspecific competition amongst the predators. This reflects an important biological feature.

4.2. A two-parameter bifurcation of system (17) with $z$ considerably more efficient than $y$.

In this subsection, we study the dynamics generated by system (17) by considering a higher value of the predation efficiency of $z$. Keeping the other parameter values fixed as defined by (19), we again consider a two-parameter bifurcation of system (17) in $(h,\beta_{1})$ parameter space. The bifurcation diagram is shown in figure 13. Comparing it with figure 10, we note that the two bifurcation diagrams are similar in many aspects, but the latter contains additional complex bifurcation points such as the generalized Hopf bifurcation.

The equilibrium $E_{xz}$ on the $xz$ plane is locally asymptotically stable in regions 1 and 2i in figure 13 and undergoes a supercritical Hopf bifurcation along the vertical line denoted by $H$. In a small parameter regime between the BR and SN curves, denoted by regions 2i and 2ii, the system has two coexistence equilibria $E_{1}^{*}$ and $E_{2}^{*}$. In region 2i, the equilibrium $E_{2}^{*}$ undergoes a subcritical Hopf bifurcation as it crosses the Hopf curve $H$, while in region 2ii, $E_{1}^{*}$ undergoes a supercritical Hopf bifurcation as it crosses the Hopf curve. A co-dimension 2 generalized Hopf bifurcation, denoted by GH, occurs where the subcritical and the supercritical Hopf curves meet. In a vicinity of the subcritical Hopf curve, a torus bifurcation curve, TR, emanates from GH and gives rise to periodic oscillations. The TR curve terminates at the fold-Hopf (FH) bifurcation. For parameter values near the TR curve in region 2i, long lasting transient dynamics with a peculiar approach to $E_{xz}$ are seen. The temporal development of the species are shown in figure 14.

Below the TR curve and above the yellow curve in figure 13, the equilibrium $E_{xz}$ is locally asymptotically stable and the system exhibits bistable dynamics with periodic orbits born out of torus bifurcation as shown in figure 15. The coexistence equilibrium $E_{2}^{*}$ in this regime is saddle-focus, with one dimensional unstable manifold, $W^{u}(E_{2}^{*})$, and two-dimensional stable manifold $W^{s}(E_{2}^{*})$. The approach of a trajectory to the periodic orbit that starts in a neighborhood of $E_{2}^{*}$ is shown in figure 15(A). Note the difference in the approaches to the periodic orbits between figures 11(A) and 15(A).

In the region bounded by the PD curve in figure 13, aperiodic and periodic MMO orbits of different signatures are observed. Interesting new MMO patterns, not of the standard types seen in systems with two slow variables and one fast variable, featuring significant variations (by orders of magnitude) in the number and amplitude of SAOs between successive spikes (see figure 18) are observed. These complicated dynamics demonstrate the uncertainty in the return time between population outbreaks, a realistic phenomena seen in many natural populations. These dynamics will be briefly discussed in the next section.

5. THE GEOMETRIC SINGULAR PERTURBATION THEORY APPROACH

In this section, we will take a geometrical singular perturbation approach to analyze the complex dynamics exhibited by system $(\ref{nondim3})$ (or equivalently system $(\ref{nondim2})$). The foundation of such geometric approach to analyze systems with a clear separation in time scales was given by Fenichel [14]. We will also review the appropriate canard theory to explain the mechanism behind the MMO patterns.

As $\zeta\to 0$ the trajectories of (13) during fast epochs approach to the solutions of the “layer equations” given by

[TABLE]

On the other hand, during slow epochs trajectories of (17) converge to the solutions of the “reduced problem” given by

[TABLE]

The subsystems (20) and (21) will be used to study the dynamics of the full system (13) or (17) (see [45] for details).

5.1. The reduced and layer problems

The algebraic equation in (21) defines the critical set

[TABLE]

where $T=\{(0,y,z):y,z\geq 0\}$ and $S=\{(x,y,z)\in{\mathbb{R}^{3}}^{+}:\phi(x,y,z)=0\}$. The critical set $\mathcal{M}$, consisting of two pieces $T$ and $S$, is the nullsurface of the layer system (20). The two pieces meet along the line $\mathcal{TC}=\{(0,y,z):y/\beta_{1}+z/\beta_{2}=1\}$ which divides the plane $T$ into two normally hyperbolic sheets $T^{a}=\{(0,y,z):\phi(0,y,z)<0\}$ and $T^{r}=\{(0,y,z):\phi(0,y,z)>0\}$. The surface $S$ is also divided into two normally hyperbolic sheets $S^{a}=S\cap\{\phi_{x}(x,y,z)<0\}$ and $S^{r}=S\cap\{\phi_{x}(x,y,z)>0\}$ by the curve $\mathcal{F}=S\cap\{\phi_{x}(x,y,z)=0\}$. A linearization of (20) at $\mathcal{M}\setminus\mathcal{F}$ and $\mathcal{M}\setminus\mathcal{TC}$ yield that the sheets $T^{a}$ and $S^{a}$ are attracting while $T^{r}$ and $S^{r}$ are repelling (see figure 16(A)). Saddle-node and transcritical bifurcations of equilibria of the fast subsystem (20) occur along the fold curve $\mathcal{F}$ and the transcritical curve $\mathcal{TC}$ respectively. By Fenichel’s theory [14], the normally hyperbolic segments of the critical manifold $\mathcal{M}$ perturb to locally invariant attracting and repelling slow manifolds ${T_{\zeta}^{a}}\cup{S_{\zeta}^{a}}$ and $T_{\zeta}^{r}\cup S_{\zeta}^{r}$ respectively for $\zeta>0$, and the slow flow restricted to these manifolds is an $O(\zeta)$ perturbation of the reduced flow on $\mathcal{M}$. However, the theory breaks down in neighborhoods of $\mathcal{F}$ and $\mathcal{TC}$ and interesting dynamics such as relaxation oscillations and MMOs occur.

Note that the reduced flow is restricted to the plane $T$ or to the surface $S$. On the plane $T$, the reduced dynamics solves the system

[TABLE]

Since $\dot{y}<-cy$ and $\dot{z}<-dz$ on $T$ with $(0,0,0)$ being the global attractor of (25), the reduced flow descends along this plane and approaches $(0,0,0)$. As the reduced flow descends, it crosses $\mathcal{TC}$ from $T^{a}$ to $T^{r}$ with finite speed, giving rise to singular canards. The plane $T$ is invariant for all $\zeta>0$, hence the connection between stable and unstable persists for $\zeta>0$, i.e. canards persist for the full system.

Since $\phi_{z}\neq 0$, by the implicit function theorem, the surface $S$ can be locally written as a graph of $z=\theta(x,y)$, i.e. $\phi(x,y,\theta(x,y))=0$. Hence we can project the dynamics of (21) onto the $(x,y)$ coordinate chart. Differentiating $\phi(x,y,z)=0$ implicitly with respect to time gives us the relationship $\phi_{x}\dot{x}+\phi_{y}\dot{y}+\phi_{z}\dot{z}=0$. Thus, the reduced flow (21) restricted to $S$, where $S$ is considered as the graph of $z=\theta(x,y)$ reads as

[TABLE]

System (26) has singularities when $\phi_{x}=0$ and its solutions blow-up in finite time at $\mathcal{F}$. Hence standard existence and uniqueness results do not hold. Points on $\mathcal{F}$ for which $\phi_{y}y\chi+\phi_{z}z\psi\neq 0$ are called “jump points” and they satisfy the “normal switching condition” [7]. At these points, a solution of (17) exits into relaxation after reaching $\mathcal{F}$ giving rise to relaxation dynamics or more commonly, referred to as boom and bust cycles in ecology [45] . On the other hand, points of $\mathcal{F}$ where the normal switching condition is violated can give rise to canards as discussed below.

To analyze the solutions where the normal switching condition fails, we rescale the time $s$ by a phase-space-dependent time transformation factor $-\phi_{x}$, i.e. $ds=-\phi_{x}dt_{s}$ [12]. This removes the finite-time blow up of solutions and system (26) transforms to the desingularized system

[TABLE]

where the overdot denotes $t_{s}$ derivatives. System (27) is topologically equivalent to system (26) on $S^{a}$. However, the phase-space-dependent time transformation reverses the orientation of the orbits on $S^{r}$, therefore, the flow of (26) on $S^{r}$ is obtained by reversing the direction of orbits of (27). Hence the reduced flow is either directed towards the fold or away from it (see figure 16(B)).

The set of equilibrium points of (27) that do not lie on the fold curve $\mathcal{F}$ (i.e. for which $\phi_{x}\neq 0$) are ordinary singularities. Such points lie on the curves $\{\chi=0\}\cap\{z=0\}\cap\{\phi_{x}\neq 0\}$, $\{\psi=0\}\cap\{y=0\}\cap\{\phi_{x}\neq 0\}$ and $\{\chi=0\}\cap\{\psi=0\}\cap\{\phi_{x}\neq 0\}$ (see figure 17). On the other hand, equilibrium points of (27) that lie on $\mathcal{F}$ are folded singularities or canard points [12, 52]. The set of folded singularities form isolated points of $\mathcal{F}$. The trajectories of the reduced system (21) may pass through the canard points and can thus cross from $S^{a}$ to $S^{r}$ with finite speed, giving rise to singular canards. The classification of a folded singularity as a folded node or a folded saddle or a folded focus or a degenerate folded node is based on the linearization of the folded singularity when considered as an equilibrium of (27) [52]. More precisely, if $\lambda_{1,2}$ are eigenvalues of the linearization of (27) at a folded singularity $p_{f}$, then $p_{f}$ is a folded focus if $\lambda_{1,2}\in\mathbb{C}$. If $\lambda_{1,2}\in\mathbb{R}$ such that $\lambda_{1}\lambda_{2}<0$, then $p_{f}$ is a folded saddle, while if $\lambda_{1},\lambda_{2}<0$, then $p_{f}$ is a folded node. Further degeneracies may occur if one of the eigenvalues pass through [math] and can give rise to folded saddle node (FSN) bifurcation of types I and II [12].

Figure 17 shows the curves of ordinary and folded singularities of system (27) as a function of the control parameter $h$, where the other parameter values are held constant. There exists two curves of ordinary singularities, namely the equilibria curves with $y$ absent and the curve of positive equilibria (coexistence equilibrium state). The latter intersects with the curve of folded singularities at a transcritical bifurcation, also referred to as the FSN II bifurcation. Here the equilibrium $E^{*}$ of the full system (13) and a folded singularity of the desingularized system (27) merge together and then split again, interchanging their type and stability. The equilibrium $E^{*}$ crosses the fold curve $\mathcal{F}$ at this bifurcation. The folded singularity switches to a folded node from a folded saddle, and $E^{*}$ switches from ordinary node to ordinary saddle.

The folded node possesses strong and weak eigenvalues $\lambda_{s}$ and $\lambda_{w}$, respectively. The singular strong canard $\gamma_{s}$ is the unique trajectory to (27) tangent to the strong eigendirection $v_{s}$, at $p$, while the singular weak canard $\gamma_{w}$ is tangential to the weak eigendirection, $v_{w}$. The fold curve $\mathcal{F}$ and the strong singular canard $\gamma_{s}$ form a trapping region (singular funnel) on $S^{a}$ as shown in figure 16(B), such that all solutions in the funnel converge to the folded node $p$. In fact, the folded node allows a sector of singular canards to flow from $S^{a}$ to $S^{r}$ and can give rise to local oscillations for $0<\zeta<<1$ [12].

5.2. Singular orbits, MMOs and relaxation oscillations

The dynamics of the layer problem (20) and the reduced flow (21) can be combined to construct singular orbits. Such orbits are used to study oscillatory dynamics such as relaxation oscillations or MMOs of the full system (17) [24]. Typically, a singular periodic orbit for relaxation oscillations in slow-fast systems with “S-shaped” critical manifolds with two folds is constructed by concatenating the reduced flow, occurring along the attractive branches of the critical manifold, with fast fibers, along which transitions between the branches occur when jump points on the fold curves are reached. However, the existence of singular canards in $T$ makes the construction of the singular orbit more subtle for system (17). In this case, the reduced flow on $S$ occurs along $S^{a}$ until it reaches a jump point on $\mathcal{F}$, where it gets connected to $T^{a}$ by a fast fiber of (20). On $T$, as discussed before, the reduced flow crosses $\mathcal{TC}$ as it descends along $T^{a}$ and stays on $T^{r}$ for a while, until it reaches a point $(0,y_{\tau_{0}},z_{\tau_{0}})\in T^{r}$, where a fast orbit concatenates with it. This phenomenon of delay is referred to as the Pontryagin’s delay of stability loss [1, 39, 40, 48]. The delay map $\mathcal{P}^{0}:T^{a}\to T^{r}$ is defined by $\mathcal{P}^{0}(y_{0},z_{0})=(y(\tau_{0}(y_{0},z_{0})),z(\tau_{0}(y_{0},z_{0})))$, where the delay $\tau_{0}$ is expressed by the integral

[TABLE]

where $(y(s),z(s))$ solves (25) with initial value $(y_{0},z_{0})\in T^{a}$. The existence of the Pontryagin’s delay of stability loss point on $T^{r}$ is incorporated in construction of the slow piece of the singular orbit [1, 48]. Singular relaxation cycles are then constructed as continuous concatenation of layer and reduced flow orbit segments and are singular representation of relaxation oscillation cycles. In a similar way, singular orbits which pass through the singular funnel and filtered through folded nodes are singular representation of MMO orbits [7].

In the parameter regime, where folded node singularities exist for the desingularized system (27), complex oscillatory dynamics such as mixed-mode oscillations (MMOs) are observed. The SAOs associated with the MMOs in this model occur near the fold $\mathcal{F}$. System (17) exhibits MMOs for different combinations of the parameter values as seen in figures 1, 4 and 5. When $\zeta>0$, a singular Hopf bifurcation occurs at a distance $O(\zeta)$ in parameter space from FSN II bifurcation (see figure 10). Moreover, the desingularized system (27) possesses a folded node singularity (see figure 17). As discussed earlier, according to the Fenichel’s theory, the attracting and repelling sheets of the critical set $\mathcal{M}$, perturb smoothly to locally invariant manifolds $S^{a}_{\zeta}$, $S^{r}_{\zeta}$, $T^{a}_{\zeta}$ and $T^{r}_{\zeta}$ away from the fold curve $\mathcal{F}$ and the transcritical curve $\mathcal{TC}$ for sufficiently small $\zeta>0$. Extending these perturbed slow manifolds by the flow into the vicinity of a folded node singularity leads to locally twisted intersection of the attracting and the repelling sheets [12], generating the small amplitude oscillations in an MMO orbit. A trajectory originating in $S^{a}_{\zeta}$ to one side of the strong canard is trapped by $S^{r}_{\zeta}$ and returns towards $S^{a}_{\zeta}$. This leads to small rotations and the primary weak canard serves as the axis of rotation. A trajectory that originates in $S^{a}_{\zeta}$ to the other side of the strong canard is prevented by $S^{r}_{\zeta}$ from returning to $S^{a}_{\zeta}$ and follows a fast direction as it jumps to $T^{a}_{\zeta}$. This region of rotation is the funnel. Figure 16(B) gives a view of the funnel in the singular limit as $\zeta\to 0$. The inset in figure 4(A) gives a zoomed view of the local dynamics near the folded node singularity for $\zeta>0$. As a trajectory gets trapped in the funnel, it initially experiences rotations imparted by the primary weak canard. This behavior persists during its passage through the funnel until it exits the vicinity of the folded node. Such small amplitude oscillations are referred to as canard-induced oscillations [7, 8]. As the trajectory exits the funnel, the saddle-focus equilibrium $E_{1}^{*}$ which lies in $S^{r}$ plays a role in generating additional oscillations. The local vector field around the equilibrium $E_{1}^{*}$ and more precisely, the unstable manifold $W^{u}(E_{1}^{*})$ now organizes the SAOs. These oscillations are referred to as singular Hopf induced SAOs. For MMOs to persist, a trajectory that eventually exits the funnel must return within the region of attraction inside the funnel. In other words, there must exist a global return mechanism. A singular periodic orbit that starts with a fast fiber segment at the folded node and returns within the singular funnel of the folded node formed from the concatenation of orbits of the reduced and layer problems is a candidate for the global return map. For parameter values considered in figure 4, such a singular orbit exists.

It is also worth mentioning that the SAOs associated with the MMO orbit in figure 1 appear quite different from figure 4. The parameter regime chosen for both figures is just past the cascade of period doubling bifurcation of the small amplitude limit cycle $\Gamma_{h}$. However the structure of the slow manifold $S^{\zeta}$ changes with $\beta_{1}$, and the amplitude of SAOs are governed by the extent to which $S_{a}^{\zeta}$ and $S_{r}^{\zeta}$ twist while guiding the path of a trajectory through the folded node regime. A tight twisting of $S_{a}^{\zeta}$ and $S_{r}^{\zeta}$ forces a decrease in amplitude of the SAOs as observed in figure 4.

More interestingly, for parameter values within the PD curve in figure 13, chaotic MMOs with two prominent epochs of SAOs are observed as shown in figure 18. The number and amplitude of the small oscillations of the MMOs could vary significantly from one LAO cycle to the next and have a combined characteristics of the SAOs seen in figure 1, figure 4 and figure 7(C). The exact mechanism behind these MMO dynamics is left for future study.

5.3. Transient MMO dynamics revisited

In Section 3, we showed that prolonged transients in form of MMOs are observed in an immediate neighborhood of FSN II bifurcation. These transients lasted for thousands of generations before the system reached its asymptotic state, namely a periodic orbit born out of a supercritical Hopf bifurcation. We noted that the SAOs in the transient MMO dynamics have long epochs, where some of the oscillations are too small to be detectable (see figure 7(C)). These MMOs are aperiodic and their chaotic nature can be illustrated by considering a suitable Poincaré map near the folded node singularity. During its passage through the funnel, these orbits visit different rotational sectors defined by the secondary canards [12], generating MMOs of mixed signatures $1^{i_{1}}1^{i_{2}}\ldots 1^{i_{k}}$, where $i_{i}\in\mathbb{N}$ is very large. To this end, we consider the Poincaré section $\Sigma_{h}$, defined by a plane transverse to the critical manifold $S$ through the folded node $p_{n}$ containing the weak eigenvector.

A zoomed view of the SAOs of an MMO trajectory along with $\Sigma_{h}$ is included in figure 19. During the epoch of SAOs, the $y$-coordinate increases while the $z$-coordinate decreases as the trajectory moves toward to $E_{1}^{*}$, until the trajectory gets repelled by $W^{u}(E_{1}^{*})$ and reaches a jump point on $\mathcal{F}$. Starting with initial condition $(0.01,0.01,0.12)$, the $(x,z)$ coordinates of the intersections of this trajectory in the direction of decreasing $x$ with the plane $\Sigma_{h}$ are recorded. The corresponding Poincaré map is shown in figure 20. The Poincaré map indicates that the trajectory repeatedly enters the funnel and stays in a close neighborhood of the folded node singularity $p_{n}=(0.3383,0.0923,0.3474)$, while twisting along the primary weak canard and getting pushed towards the equilibrium $E_{1}^{*}=(0.3299,0.1004,0.3378)$. Moreover, a zoomed view of the map shows that the trajectory while moving towards $E_{1}^{*}$, intersects with the Poincaré section at innumerous points forming a dense array of segments as shown in figure 20(B). A detailed analysis of the underlying mechanism responsible for transition from the transient state to the asymptotic state is done in the companion paper [47].

6. Discussion and future outlook

In this paper, we used singular perturbation theory, theory of canards and numerical simulations to study complex oscillatory dynamics that arise in a three species predator-prey model in a two-trophic ecosystem. Taking into account that each species operate on different timescales, separation of timescales is introduced in the model, which gives rise to a slow-fast system. The presence of intraspecific and direct interference competition terms and the nonlinear functional responses of the predators contribute to the complexities in the model and make the dynamics more realistic. Two-parameter bifurcation diagrams of the system reveal several local and global bifurcations such as period-doubling, torus, saddle-node bifurcation for equilibria, fold-Hopf and generalized Hopf bifurcations, reflecting the rich structure possessed by the system.

The model studied in this paper is an extension of the model considered in [45], where interesting transient dynamics were observed. One such example was a sudden transition from small amplitude oscillations to large amplitude relaxation oscillations. Such transients can offer an alternative explanation to regime shifts [19] and appear in many ecosystems [33]. In this paper, long lasting transients in form of chaotic MMOs are observed near FSN II bifurcation, which eventually asymptote to a stable limit cycle born out of a supercritical Hopf bifurcation. Though a supercritical Hopf bifurcation is considered as a“safe bifurcation” in the asymptotic limit, the dynamics exhibited by this model shows that this may not be the case on an ecological timescale. The presence of a folded node singularity in vicinity of an unstable equilibrium of saddle-focus type gives rise to complexities. The folded node allows certain trajectories on the attracting slow manifold to cross into the repelling slow manifold by creating a funnel. The passage through the funnel induces SAOs in MMOs. Moreover, the proximity of the stable manifold of the saddle-focus equilibrium to the primary weak canard brings a trajectory close to a neighborhood of the equilibrium. As a result, trajectories follow the unstable manifold of the equilibrium leading to additional small rotations before jumping to the other attracting branch of the slow manifold leading to a large excursion in phase space. Thus, finding early warning signs that could help predict when a population will drastically increase is very important. This requires a detailed understanding of the local mechanism that generates the MMOs and finding the precise locations of the “jump points” on the fold curve that allow switching from one attracting branch of the slow manifold to the other. This has been considered in a companion paper [47].

The predation efficiencies, $\beta_{1},\beta_{2}$, defined as the ratios of the semi-saturation constants to the carrying capacity of the prey, also played important roles in governing the dynamics. Additional co-dimension 2 bifurcations such as the generalized Hopf bifurcation occurs when one of the predators is considered to be more efficient than the other. Furthermore, complicated MMO dynamics with different characteristics timescales and epochs of SAOs between two LAOs are observed as the intraspecific competition strength is varied when $\beta_{2}$ is assumed to be significantly higher than $\beta_{1}$. These types of oscillations are unique in the framework of MMOs and represent the inherent uncertainties associated with population dynamics and warrants further analysis. The structure of the slow manifold varies with the predation efficiencies and thereby influence the resulting dynamics. In fact, the occurrence (or absence) of MMOs, the epochs of the SAOs in an MMO orbit, and hence formation of canards depend on the twisting properties of the attracting and repelling branches of the slow manifold. Holding the intraspecific competition rate and $\beta_{2}$ constant, we also noted that as the predation efficiency of $y$ increases, the competition between the predators go up and eventually leads to competitive exclusion of $y$. The emergence of this boundary equilibrium through torus-like oscillations is remarkable. Such form of oscillations can be related to permanence or persistence in ecological communities and perhaps bear some connection with bursting phenomena, most commonly seen in neuroscience models [20].

The transient basins, defined as the set of initial conditions which produce chaotic MMO transients whose duration is within a given interval, form complicated structures in phase space, and it is plausible that a small perturbation (such as an environmental fluctuation) can very easily lead to a transition to another basin. Ecologically, this may be helpful to control population outbreaks, as an external intervention can shift the dynamics from an excited state to a more stable state.

The Koper model [22], a reduced form of the Hodgkin-Huxley model [12] and a reduced neuronal competition model [8] are amongst the few examples of slow-fast systems where a detailed analysis has been carried out near the FSN II points. In [8], the authors studied the interaction of canard and singular Hopf mechanism for generation of SAOs near the FSN II bifurcation. In their work, they classified the SAOs in an MMO orbit as “canard-induced” or induced by the local vector field around a saddle-focus equilibrium by numerically computing the “way-in/way-out” function [26] which describes the maximal delay expected for generic solutions passing through FSN singularity. Similar analysis can be possibly carried out in this model in the regime far from the fold-Hopf bifurcation, where the SAOs in an MMO trajectory are distinctly detectable. Further analysis has been done on the Koper model to study the existence of a homoclinic orbit [17, 35]. The detection of trajectories with extraordinary large number of small oscillations between large amplitude oscillations and a chaotic return map suggests that the trajectories possibly lie in a small neighborhood of a Shil’nikov homoclinic orbit in this model. Explicitly detecting a homoclinic orbit to the saddle-focus equilibrium remains for future study. The unfolding of the fold-Hopf bifurcation and numerical investigation of global bifurcations responsible for such chaotic dynamics are also subjects of future study.

7. Acknowledgement

I am grateful to Andrew Morozov for providing constructive feedback on the manuscript. I also thank Chris Cosner and Sebastian Wieczorek for helpful conversations at the workshop in “New Mathematical Methods for Complex Systems in Ecology” at Banff International Research Station in July 2019. I am thankful to Ryan Lawson, a former student of Georgia College for helping with a numerical code, Saikat Chakraborty Thakur for several helpful discussions and Georgia College Faculty Scholarship Support Program for supporting this research.

Remark 7.1**.**

Most of the numerical simulations in this paper were done in MATLAB. We used the predefined routine ODE$45$ with relative and absolute error tolerances $10^{-11}$ and $10^{-12}$ respectively.