Canard Solutions in a Predator-Prey Model

Guojian Lin

School of Mathematics, Renmin University of China, Beijing, China.

DOI: 10.4236/jamp.2022.105116                

Abstract

The canard explosion phenomenon in a predator-prey model with Michaelis-Menten functional response is analyzed in this paper by employing the geometric singular perturbation theory. First, some turning points, such as, fold point, transcritical point, pitchfork point, canard point, are identified; then Hopf bifurcation, relaxation oscillation, together with the canard transition from Hopf bifurcation to relaxation oscillation are discussed.

Keywords

Canard Explosion, Relaxation Oscillation, Predator-Prey Model, Geometric Singular Perturbation Theory

1. Introduction

To formulate the population dynamics, in 1963, MacArthur and Rosenzweig proposed the following predator-prey model with Michaelis-Menten functional response [1]

$\begin{array}{l}\stackrel{\dot{}}{u}=\gamma u\left(1-\frac{u}{K}\right)-\frac{m}{\mu }\left(\frac{vu}{a+u}\right),\\ \stackrel{\dot{}}{v}=v\left(\frac{mu}{a+u}-d\right),\end{array}$ (1.1)

where the function $v\left(t\right)$ is the population of the predator at time t, the function $u\left(t\right)$ is the population of the prey at time t, m is the maximum growth (birth) rate of the predator, d is the death rate of the predator, $\mu$ is the yield factor of the predator feeding on the prey and a is the half-saturation constant of the predator, which is the prey density at which the functional response of the predator is half maximal. The parameters $\gamma$ and K are the intrinsic rate of increase and the carrying capacity for the prey population, respectively. The parameters $\gamma \mathrm{,}K\mathrm{,}m\mathrm{,}\mu \mathrm{,}a$ and d are positive constants. Additionally, it is assumed that $u\left(t\right)\ge \mathrm{0,}v\left(t\right)\ge 0$ for biological meanings.

Using the following rescaling

$\epsilon =\frac{1}{\gamma },\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\beta =\frac{a}{K},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}x=\frac{u}{K},\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}y=\frac{v}{\gamma \mu K},$ (1.2)

it follows that system (1.1) becomes

$\begin{array}{l}\epsilon \stackrel{\dot{}}{x}=x\left(1-x-\frac{my}{\beta +x}\right),\\ \stackrel{\dot{}}{y}=y\left(\frac{mx}{\beta +x}-d\right).\end{array}$ (1.3)

In this paper, it is assumed that $\gamma$ is sufficiently large; then $\epsilon$ is a small parameter; therefore Equation (1.3) is a standard singularly perturbed system. From biological meanings, this assumption implies the prey population in the model (1.3) grows much faster than the predator population.

By switching to the fast time scale

$\tau =\frac{t}{\epsilon },$ (1.4)

one obtains the following equivalent system

$\begin{array}{l}{x}^{\prime }=x\left(1-x-\frac{my}{\beta +x}\right),\\ y^{\prime }=\epsilon y\left(\frac{mx}{\beta +x}-d\right).\end{array}$ (1.5)

If $m>d$, let $\rho =\frac{\beta d}{m-d}$, then $\rho >0$ and system (1.5) are rewritten as the follows

$\begin{array}{l}{x}^{\prime }=x\left(1-x-\frac{mdy}{\rho \left(m-d\right)+dx}\right),\\ y^{\prime }=\epsilon y\left(\frac{mdx}{\rho \left(m-d\right)+dx}-d\right).\end{array}$ (1.6)

Accordingly, the results in [2] [3] [4] are reformulated in the following form:

1) If $m\le d$, then the equilibrium $\left(\mathrm{1,0}\right)$ of system (1.5) is asymptotically stable.

2) If $\rho \ge 1$, then the equilibrium $\left(\mathrm{1,0}\right)$ of system (1.6) is asymptotically stable.

3) If $\frac{d}{m+d}\le \rho <1$, then the equilibrium $\left(\rho \mathrm{,}\frac{\rho \left(1-\rho \right)}{d}\right)$ of system (1.6) is asymptotically stable.

4) If $0<\rho <\frac{d}{m+d}$, then the equilibrium $\left(\rho \mathrm{,}\frac{\rho \left(1-\rho \right)}{d}\right)$ of system (1.6) is unstable and system (1.6) possesses a unique large-amplitude periodic solution.

In this paper, it is always assumed that $m>d$ and $0<\rho <\max \left\{\frac{d}{m-d},1\right\}$.

Additionally, it is also assumed $x\ge \mathrm{0,}y\ge 0$ in system (1.6) for biological meanings.

By using the geometric singular perturbation theory, it will be proved in this paper that the canard explosion phenomenon happens in system (1.6) as the parameter

$\rho$ decreases through $\rho =\frac{d}{m+d}$. This canard explosion phenomenon can

explain the reason why the sudden transition from a small-amplitude periodic

solution, which bifurcates from the equilibrium $\left(\rho \mathrm{,}\frac{\rho \left(1-\rho \right)}{d}\right)$ via the supercritical Hopf bifurcation at $\rho =\frac{d}{m+d}$, to a large-amplitude relaxation oscillation which emerges at $\rho <\frac{d}{m+d}$.

There are a great deal of articles [5] - [10], which are related to study the dynamics of predator-prey systems, such as bifurcations, stability, and so on.

Canard solutions were first analyzed by Benoit, Callot, Diener and Diener [11] using non-standard analysis in van der Pol equations. A canard solution is a solution of a singularly perturbed system which follows an attracting slow manifold, passes close to a non-hyperbolic point of the critical manifold, and then follows, rather surprisingly, a repelling slow manifold for a considerable amount of time before being repelled. The existence of a canard solution can lead to canard explosion, that is, a transition from a small limit cycle to a relaxation oscillation through a sequence of canard cycles upon variation of a parameter. Afterward, Eckhaus [12] studied the existence of canard solutions for van der Pol equation by employing the method of matched asymptotic expansion. A breakthrough in geometric explanation of canard cycles and canard explosion came with the work of Dumortier and Roussarie [13], who analyzed these phenomena in van der Pol’s equation by means of blow up technique and foliation of center manifolds in detail. From the work of Dumortier and Roussarie, it became apparent that blow up technique was the right tool for analyzing non-hyperbolic points of the slow manifold in a singularly perturbed system. Motivated by their work, Krupa and Szmolyan extended the standard normally hyperbolic geometric singular perturbation [14] [15] to non-hyperbolic points [16] [17] [18] by employing the blow up technique. Recently, canard solutions of a singularly perturbed system are extensively studied [19] - [24]. An introduction to basic knowledge on the geometric singular perturbation theory can be also founded in [25].

The paper is organized as follows. Section 2, Section 3 and Section 5 identify the fold point, the transcritical point, the Pitchfork Point and the canard point of system (1.6) respectively; Section 4 discusses Hopf bifurcations and relaxation oscillations of system (1.6); Section 6 analyzes the canard explosion phenomenon of system (1.6). Finally, some concluding remarks are given in Section 7.

2. Fold Point

Let $u=x-\frac{d-\rho \left(m-d\right)}{2d},v=y-\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}$, system (1.6) becomes

$\begin{array}{l}{u}^{\prime }=\left(u+\frac{d-\rho \left(m-d\right)}{2d}\right)\left(\frac{-2d{u}^2-2dmv}{2du+d+\rho \left(m-d\right)}\right),\\ v^{\prime }=\epsilon \left(v+\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)\left(\frac{2dmu+md-\rho m\left(m-d\right)}{2du+d+\rho \left(m-d\right)}-d\right).\end{array}$ (2.1)

Let $x=-u,y=-v$, then system (2.1) reduces to

$\begin{array}{l}{x}^{\prime }=\left(x-\frac{d-\rho \left(m-d\right)}{2d}\right)\left(\frac{2dmy-2d{x}^2}{d+\rho \left(m-d\right)-2dx}\right),\\ y^{\prime }=\epsilon \left(\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}-y\right)\left(\frac{2mdx-md+\rho m\left(m-d\right)}{d+\rho \left(m-d\right)-2dx}+d\right).\end{array}$ (2.2)

where $x\le \frac{d-\rho \left(m-d\right)}{2d}$ and $y\le \frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}$.

Let

$\begin{array}{l}f\left(x,y,\epsilon \right)=\left(x-\frac{d-\rho \left(m-d\right)}{2d}\right)\left(\frac{2dmy-2d{x}^2}{d+\rho \left(m-d\right)-2dx}\right),\\ g\left(x,y,\epsilon \right)=\left(\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}-y\right)\left(\frac{2mdx-md+\rho m\left(m-d\right)}{d+\rho \left(m-d\right)-2dx}+d\right).\end{array}$

then system (2.2) can be rewritten as

$\begin{array}{l}{x}^{\prime }=f\left(x,y,\epsilon \right),\\ y^{\prime }=\epsilon g\left(x,y,\epsilon \right).\end{array}$ (2.3)

Setting $\epsilon =0$ in system (2.3) results in the layer problem

$\begin{array}{l}{x}^{\prime }=f\left(x,y,0\right),\\ y^{\prime }=0.\end{array}$ (2.4)

In term of the time rescaling $\tau =\frac{t}{\epsilon }$, system (2.3) becomes

$\begin{array}{l}\epsilon \stackrel{\dot{}}{x}=f\left(x,y,\epsilon \right),\\ \stackrel{\dot{}}{y}=g\left(x,y,\epsilon \right).\end{array}$ (2.5)

Setting $\epsilon =0$ in system (2.5) results in the reduced problem

$\begin{array}{l}0=f\left(x,y,0\right),\\ \stackrel{\dot{}}{y}=g\left(x,y,0\right).\end{array}$ (2.6)

Let

$S=\left\{\left(x,y\right):\left(x-\frac{d-\rho \left(m-d\right)}{2d}\right)\left(my-x^2\right)=0\right\}$

be the slow manifold, which consists of two parts $S_1$ and $S_2$, where

$S_1=\left\{\left(x,y\right)\in S:y=\frac{1}{m}{x}^2\right\}$

and

$S_2=\left\{\left(x,y\right)\in S:x=\frac{d-\rho \left(m-d\right)}{2d}\right\}.$

Let $S_1^a=\left\{\left(x,y\right)\in S_1:x<0\right\}$ and $S_1^r=\left\{\left(x,y\right)\in S_1:x>0\right\}$.

Let

$S_2^a=\left\{\left(x,y\right)\in S_2:y<\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right\}$

and

$S_2^r=\left\{\left(x,y\right)\in S_2:y>\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right\}.$

Assume that $0<\rho <\frac{d}{m+d}$, then it can verified that

$\begin{array}{l}f\left(0,0,0\right)=0,g\left(0,0,0\right)=\frac{\left(m-d\right)\left[d+\rho \left(m-d\right)\right]\left[\rho \left(m+d\right)-d\right]}{4m{d}^2}<0,\\ \frac{\partial f}{\partial x}\left(0,0,0\right)=0,\frac{\partial f}{\partial y}\left(0,0,0\right)=-\frac{m\left[d-\rho \left(m-d\right)\right]}{d+\rho \left(m-d\right)}>0,\\ \frac{{\partial }^{2}f}{\partial x^2}\left(0,0,0\right)=\frac{2\left[d-\rho \left(m-d\right)\right]}{d+\rho \left(m-d\right)}>0.\end{array}$

Therefore, by the definition of a fold point [16], (0, 0) is a fold point of system (2.3).

Under the assumption that $0<\rho <\frac{d}{m+d}$, it can be verified that the branch

$S_1^a$ is attracting and the branch $S_1^r$ is repelling for the layer problem. The origin (0, 0) is nonhyperbolic, weakly attracting from the left and weakly repelling to the right. Moreover, the reduced flow on $S_a^1$ and $S_r^1$ is directed towards the fold point (0, 0), see Figure 1 for the dynamics of the layer problem and the reduced problem.

Figure 1. Slow-fast dynamics of system (2.3) for $\epsilon =0$ in the case that $0<\rho <\frac{d}{m+d}$.

The standard normally hyperbolic geometric singular perturbation [14] implies that outside an arbitrarily small neighborhood of (0, 0), the manifolds $S_1^a$ and $S_1^r$ perturb smoothly to locally invariant manifolds $S_1^{a\mathrm{,}\epsilon }$ and $S_1^{r\mathrm{,}\epsilon }$, which are simply solutions to system (2.3).

Let

${\Delta }^{out}=\left\{\left(\delta \mathrm{,}y\right)\mathrm{,}y\in J\right\}$

be a section transverse to the fast fiber, where $J\in ℝ$ is a suitable interval and $\delta >0$ is a suitable constant, see Figure 1.

By theorem 2.1 in [16], it follows that

Proposition 2.1. There exists ${\epsilon }_0>0$ such that for $\epsilon \in \left(\mathrm{0,}{\epsilon }_0\right)$ , the manifold $S_1^{a\mathrm{,}\epsilon }$ passes through ${\Delta }^{out}$ at a point $\left(\delta \mathrm{,}h\left(\epsilon \right)\right)$ with $h\left(0\right)=0$ .

3. Transcritical Point and Pitchfork Point

Let $u=x-\frac{d-\rho \left(m-d\right)}{2d},v=y-\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}$, then system (2.2) becomes

$\begin{array}{l}{u}^{\prime }=\frac{mduv-d{u}^3-u^2\left[d-\rho \left(m-d\right)\right]}{\rho \left(m-d\right)-du},\\ v^{\prime }=\epsilon \left(\frac{\rho \left(m-d\right)}{md}-v\right)\left(\frac{mdu}{\rho \left(m-d\right)-du}+d\right).\end{array}$ (3.1)

Let $x=-u,y=v$, then system (3.1) becomes

$\begin{array}{l}{x}^{\prime }=\frac{mdxy-d{x}^3+x^2\left[d-\rho \left(m-d\right)\right]}{\rho \left(m-d\right)+dx},\\ y^{\prime }=\epsilon \left(\frac{\rho \left(m-d\right)}{md}-y\right)\left(\frac{-mdx}{\rho \left(m-d\right)+dx}+d\right).\end{array}$ (3.2)

Let

$\begin{array}{l}H\left(x,y,\epsilon \right)=\frac{mdxy-d{x}^3+x^2\left[d-\rho \left(m-d\right)\right]}{\rho \left(m-d\right)+dx},\\ I\left(x,y,\epsilon \right)=\left(\frac{\rho \left(m-d\right)}{md}-y\right)\left(\frac{-mdx}{\rho \left(m-d\right)+dx}+d\right).\end{array}$

then it follows that system (3.2) can be rewritten as

$\begin{array}{l}{x}^{\prime }=H\left(x,y,\epsilon \right),\\ y^{\prime }=\epsilon I\left(x,y,\epsilon \right).\end{array}$ (3.3)

Under the assumption that $0<\rho <\frac{d}{m-d}$, it can be calculated that

$\begin{array}{l}H\left(0,0,0\right)=0,I\left(0,0,0\right)=\frac{\rho \left(m-d\right)}{m}>0,\frac{\partial H}{\partial x}\left(0,0,0\right)=0,\\ \frac{\partial H}{\partial y}\left(0,0,0\right)=0,\frac{{\partial }^{2}H}{\partial x^2}\left(0,0,0\right)=\frac{2\left[d-\rho \left(m-d\right)\right]}{\rho \left(m-d\right)}>0,\\ \frac{{\partial }^{2}H}{\partial y^2}\left(0,0,0\right)=0,\frac{{\partial }^{2}H}{\partial x\partial y}\left(0,0,0\right)=\frac{md}{\rho \left(m-d\right)}>0.\end{array}$

It follows that

$\begin{array}{c}\left|\begin{array}{cc}\frac{{\partial }^{2}H}{\partial x^2}\left(0,0,0\right)& \frac{{\partial }^{2}H}{\partial x\partial y}\left(0,0,0\right)\\ \frac{{\partial }^{2}H}{\partial x\partial y}\left(0,0,0\right)& \frac{{\partial }^{2}H}{\partial y^2}\left(0,0,0\right)\end{array}\right|=\left|\begin{array}{cc}\frac{2\left[d-\rho \left(m-d\right)\right]}{\rho \left(m-d\right)}& \frac{md}{\rho \left(m-d\right)}\\ \frac{md}{\rho \left(m-d\right)}& 0\end{array}\right|\\ =-\frac{m^2{d}^2}{{\rho }^2{\left(m-d\right)}^2}<0.\end{array}$

By the definition of a transcritical point in [17], it can be seen that the point (0, 0) is a transcritical point of system (3.3), which implies that the point

$\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ is a transcritical point of system (2.2).

Remark 3.1. If $\rho =\frac{d}{m-d}$ , then $\frac{{\partial }^{2}H}{\partial x^2}\left(0,0,0\right)=0$ . Furthermore, it is can be calculated that $\frac{{\partial }^{3}H}{\partial x^3}\left(0,0,0\right)=-\frac{6d^2}{{\rho }^2{\left(m-d\right)}^2}<0$ . By the definition of a pitchfork

point in [17], it can be seen that the point (0, 0) is a pitchfork point of system

(3.3), which implies that the point $\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ is a pitchfork point of system (2.2).

If $\rho >\frac{d}{m-d}$, then it can be seen that the point $\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ is also a transcritical point of system (2.2).

However,in the cases that $\rho =\frac{d}{m-d}$ and $\rho >\frac{d}{m-d}$ ,canard explosion phenomena do not happen in (2.2).

Under the assumption that $0<\rho <\frac{d}{m-d}$, it can be verified that the branch

$S_2^a$ is attracting and the branch $S_2^r$ is repelling for the layer problem. The point (1, 1) of system (2.2) is nonhyperbolic, weakly repelling from the left and weakly attracting to the right, see Figure 1.

The standard normally hyperbolic geometric singular perturbation [14] implies that outside an arbitrarily small neighborhood of (1, 1) of system (2.2), the manifolds $S_2^a$ and $S_2^r$ perturb smoothly to locally invariant manifolds $S_2^{a\mathrm{,}\epsilon }$ and $S_2^{r\mathrm{,}\epsilon }$. In the following, it will be analyzed that how does $S_2^{a\mathrm{,}\epsilon }$ pass through

a neighbourhood of the transcritical point $\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ of system (2.2).

Let $u=x+\frac{md}{2\left[d-\rho \left(m-d\right)\right]}y,v=y$, then system (3.2) becomes

$\begin{array}{l}{u}^{\prime }=\frac{d-\rho \left(m-d\right)}{\rho \left(m-d\right)}{u}^2-\frac{m^2{d}^2}{4\rho \left(m-d\right)\left[d-\rho \left(m-d\right)\right]}{v}^2\\ +\frac{d\rho \left(m-d\right)}{2\left[d-\rho \left(m-d\right)\right]}\epsilon +h_1\left(u,v,\epsilon \right),\\ v^{\prime }=\epsilon \left(\frac{\rho \left(m-d\right)}{m}+h_2\left(u,v,\epsilon \right)\right),\end{array}$ (3.4)

where $h_1\left(u,v,\epsilon \right)=O\left(u^3,u^{2}v,u{v}^2,v^3,\epsilon u,\epsilon v,{\epsilon }^2\right)$,

$h_2\left(u\mathrm{,}v\mathrm{,}\epsilon \right)=O\left(u\mathrm{,}v\mathrm{,}\epsilon \right)$.

Letting

$\tilde{x}=\frac{d-\rho \left(m-d\right)}{\rho \left(m-d\right)}u,\tilde{y}=\frac{md}{2\rho \left(m-d\right)}v,\tilde{\epsilon }=\frac{d}{2}\epsilon ,$ (3.5)

and substituting (3.5) into Equation (3.4), then by directly calculating and dropping the tildes, Equation (3.4) becomes

$\begin{array}{l}{x}^{\prime }=x^2-y^2+\sigma \epsilon +h_1\left(x,y,\epsilon \right),\\ y^{\prime }=\epsilon \left[1+h_2\left(x,y,\epsilon \right)\right]\end{array}$ (3.6)

with $\sigma =1$, $h_1\left(x\mathrm{,}y\mathrm{,}\epsilon \right)=O\left(x^3\mathrm{,}{x}^{2}y\mathrm{,}x{y}^2\mathrm{,}{y}^3\mathrm{,}\epsilon x\mathrm{,}\epsilon y\mathrm{,}{\epsilon }^2\right)$ and $h_2\left(x\mathrm{,}y\mathrm{,}\epsilon \right)=O\left(x\mathrm{,}y\mathrm{,}\epsilon \right)$.

Therefore, by using a result in [17], it follows that

Proposition 3.2. There exists ${\epsilon }_0>0$ and a function ${\sigma }_c\left(\sqrt{\epsilon }\right)$ with ${\sigma }_c\left(0\right)=1$ such that for $\sigma ={\sigma }_c\left(\sqrt{\epsilon }\right)$ , the slow manifold $S_2^{a\mathrm{,}\epsilon }$ extend to $S_2^{r\mathrm{,}\epsilon }$ for sufficiently small $\epsilon >0$ .

4. Relaxation Oscillation and Hopf Bifurcation

Based on the local dynamics nearby the fold point (0, 0) and the transcritical

point $\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ of system (2.2), the following result can be obtained.

Theorem 4.1. Assume that $0<\rho <\frac{d}{m+d}$ , then for sufficiently small $\epsilon >0$ , system (2.2) has a stable large-amplitude limit cycle ${\Gamma }_{\epsilon }$ .

Proof. Let ${\Delta }^{in}$ be a section of the flow defined as a small horizontal interval

intersecting $S_1^a$ at a point between $\left(-\frac{d+\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ and (0,

0), see Figure 1. Consider tracking a trajectory starting in ${\Delta }^{in}$ for $0<\epsilon \ll 1$. Initially this trajectory will be attracted to $S_1^{a\mathrm{,}\epsilon }$ and then pass beyond the fold point (0, 0) until it reach the section ${\Delta }^{out}$. As this trajectory arrives in the vicinity of $S_2^a$, it will be attracted to $S_2^{a\mathrm{,}\epsilon }$ and then pass beyond the transcritical point (1, 1). Rather surprising, this trajectory will follow $S_2^{r\mathrm{,}\epsilon }$ for a considerable amount of time until it is repelled by $S_2^{r\mathrm{,}\epsilon }$. Therefore this trajectory will come close to $S_1^a$ and it will follow $S_1^{a\mathrm{,}\epsilon }$ until it reaches ${\Delta }^{in}$. Let $\pi \mathrm{<mo>\; :\; </mo>}{\Delta }^{in}\to {\Delta }^{in}$ be the return map. By the geometric singular perturbation theory, it follow that for $0<\epsilon \ll 1$, $\pi$ is a contraction map. By the implicit function theorem, there exists a unique and attracting fixed point of $\pi$ in ${\Delta }^{in}$. This fixed point gives rise to a stable large-amplitude limit cycle ${\Gamma }_{\epsilon }$.

Theorem 4.2. There exist ${\rho }_0=\frac{d}{m+d}$ such that system (2.2) has a unique and stable small-amplitude limit cycle bifurcating from the equilibrium (0, 0) via the supercritical Hopf bifurcation for $\rho <{\rho }_0=\frac{d}{m+d}$ .

Proof. System (2.2) has an equilibrium $E=\left(\frac{d-\rho \left(m-d\right)}{2d},\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$. The linearization of system (2.2) at E has the following form

$\left|\begin{array}{cc}\frac{d-\rho \left(m+d\right)}{m}& -d\\ \frac{\left(1-\rho \right)\left(m-d\right)\epsilon }{m}& 0\end{array}\right|$

which has eigenvalues

$\alpha \left(\rho \right)\pm i\beta \left(\rho \right)=\frac{d-\rho \left(m+d\right)}{2m}\pm \frac{i\sqrt{4md\left(1-\rho \right)\left(m-d\right)\epsilon -{\left[d-\rho \left(m+d\right)\right]}^2}}{2m}.$

If $\rho ={\rho }_0=\frac{d}{m+d}$, then $\alpha \left({\rho }_0\right)=0$, $\beta \left({\rho }_0\right)=\sqrt{\frac{d\left(m-d\right)\epsilon }{m+d}}>0$, ${\alpha }^{\prime }\left({\rho }_0\right)=-\frac{m+d}{2m}\ne 0$.

Lengthy calculations show that the first Liapunov coefficient

$l_1\left(0\right)=-\frac{d{\left(m+d\right)}^2}{4m\left(m-d\right)\epsilon }<0.$

Therefore, the results in theorem 4.2 are justified.

5. Canard Point

In this section, it is assumed that $\rho ={\rho }_0=\frac{d}{m+d}$.

Let

$\begin{array}{l}f\left(x,y,\rho ,\epsilon \right)=\left(x-\frac{d-\rho \left(m-d\right)}{2d}\right)\left(\frac{2dmy-2d{x}^2}{d+\rho \left(m-d\right)-2dx}\right),\\ g\left(x,y,\rho ,\epsilon \right)=\left(\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}-y\right)\left(\frac{2mdx-md+\rho m\left(m-d\right)}{d+\rho \left(m-d\right)-2dx}+d\right),\end{array}$

then system (2.2) can be rewritten as the following form.

$\begin{array}{l}{x}^{\prime }=f\left(x,y,\rho ,\epsilon \right),\\ y^{\prime }=\epsilon g\left(x,y,\rho ,\epsilon \right).\end{array}$ (5.1)

It can verified that

$\begin{array}{l}f\left(0,0,{\rho }_0,0\right)=0,g\left(0,0,{\rho }_0,0\right)=0,\frac{\partial f}{\partial x}\left(0,0,{\rho }_0,0\right)=0,\frac{\partial f}{\partial y}\left(0,0,{\rho }_0,0\right)=-d<0,\\ \frac{\partial g}{\partial x}\left(0,0,{\rho }_0,0\right)=\frac{m-d}{m+d}>0,\frac{\partial g}{\partial \rho }\left(0,0,{\rho }_0,0\right)=\frac{m-d}{2d}>0,\frac{{\partial }^{2}f}{\partial x^2}\left(0,0,{\rho }_0,0\right)=\frac{2d}{m}>0.\end{array}$

By the definition of a canard point in [16], it can be seen that the point $\left(x,y\right)=\left(0,0\right)$ is a canard point of system (5.1).

The reduced dynamics on $S_1$ is governed by the equation

$\begin{array}{l}\stackrel{\dot{}}{x}=\left(\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4d^2}-x^2\right)\frac{d\left(m-d\right)}{d+\rho \left(m-d\right)-2dx}.\hfill \end{array}$ (5.2)

It follows that the right-hand side of system (5.2) is a smooth function at the origin. Let $x_0\left(t\right)$ denote a maximal solution of system (5.2) with the property $x_0\left(0\right)=0$. It follows that $x_0\left(t\right)$ exists and passes through the origin, see Figure 2 for the dynamics of the layer problem and the reduced problem.

By theorem 3.1 in an article [16], it follows that

Proposition 5.1. There exists ${\epsilon }_0>0$ and a smooth function ${\rho }_0+{\rho }_c\left(\sqrt{\epsilon }\right)$ defined on $\left[\mathrm{0,}{\epsilon }_0\right]$ such that for $\epsilon \in \left(\mathrm{0,}{\epsilon }_0\right]$ , a solution starting in $S_1^{a\mathrm{,}\epsilon }$ connects to $S_1^{r\mathrm{,}\epsilon }$ if and only if $\rho ={\rho }_0+{\rho }_c\left(\sqrt{\epsilon }\right)$ with ${\rho }_c\left(0\right)=0$ .

6. Canard Explosion

Let

$\phi \left(x,y,\epsilon \right)=\frac{2dmy-2d{x}^2}{d+\rho \left(m-d\right)-2dx}.$

Theorem 6.1. For any $y_0\in \left[\mathrm{0,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)\subset S_2^a$ , there exists a unique $\tau >0$ such that

Figure 2. Slow-fast dynamics of system (2.2) for $\epsilon =0$ in the case that $\rho ={\rho }_0=\frac{d}{m+d}$.

${\displaystyle {\int }_0^{\tau }}\phi \left(\frac{d-\rho \left(m-d\right)}{2d},y\cdot s,0\right)\text{d}s=0,$

where $y\cdot s$ denotes the solution of system (2.5) on $S_2$ at $\epsilon =0$ .

Proof. For the limiting slow dynamics on $S_2$, system (2.5) is reduced to

$\stackrel{\dot{}}{y}=-dy+\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4md}.$

Therefore, for any $y_0\in S_2^a$, it can be calculated that

$y\left(s\right)=\left(y_0-\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right){\text{e}}^{-ds}+\frac{{\left[d+\rho \left(m-d\right)\right]}^2}{4m{d}^2}.$

It follows that

$\begin{array}{l}{\displaystyle {\int }_0^{\tau }}\phi \left(\frac{d-\rho \left(m-d\right)}{2d},y\cdot s,0\right)\text{d}s={\displaystyle {\int }_0^{\tau }}\left(1+\frac{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}{4d\rho \left(m-d\right)}{\text{e}}^{-ds}\right)\text{d}s\\ =-\frac{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}{4\rho \left(m-d\right)}\left({\text{e}}^{-d\tau }-1-\frac{4d^2\rho \left(m-d\right)}{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}\tau \right)=0,\end{array}$

which is equivalent to

${\text{e}}^{-d\tau }-1-\frac{4d^2\rho \left(m-d\right)}{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}\tau =0.$

Let

$F\left(\tau \right)={\text{e}}^{-d\tau }-1-\frac{4d^2\rho \left(m-d\right)}{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}\tau .$

Then for any $y_0\in \left[\mathrm{0,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)\subset S_2^a$, it follows that

$F\left(0\right)=0,F^{\prime }\left(\tau \right)=-d{\text{e}}^{-d\tau }-\frac{4d^2\rho \left(m-d\right)}{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}.$

Therefore, there exists a unique ${\tau }_0=-\frac{1}{d}\ln \left(-\frac{4d\rho \left(m-d\right)}{4m{d}^2{y}_0-{\left[d+\rho \left(m-d\right)\right]}^2}\right)>0$ such that

$\text{For}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\tau >{\tau }_0,F^{\prime }\left(\tau \right)>0;\tau <{\tau }_0,F^{\prime }\left(\tau \right)<0.$

It follows that there exists a unique $\tau >0$ such that $F\left(\tau \right)=0$. Thus theorem 6.1 is proved.

Define a map $P\mathrm{<mo>\; :\; </mo>}{S}_2^a\to S_2^r$ by $P\left(y_0\right)=y_1=y\cdot \tau$. where $y\cdot \tau$ denotes the solution on $S_2$ at $\epsilon =0$ and $\tau$ is determined by theorem 6.1.

Remark 6.2. At $y_0=\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}$ , define $P\left(y_0\right)=\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}$ .

Define singular canard cycles

$\Gamma \left(s\right)=\left\{\left(x,\frac{1}{m}{x}^2\right):x\in \left[-\sqrt{ms},\sqrt{ms}\right]\right\}\cup \left\{\left(x,s\right):x\in \left[-\sqrt{ms},\sqrt{ms}\right]\right\},$

for $s\in \left[\mathrm{0,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right]$,

$\begin{array}{c}\Gamma \left(s\right)=\left\{\left(x,\frac{1}{m}{x}^2\right):x\in \left[-\sqrt{m\cdot P\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)},\sqrt{m\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)}\right]\right\}\\ \cup \left\{\left(x,\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right):x\in \left[\sqrt{m\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)},\frac{d-\rho \left(m-d\right)}{2d}\right]\right\}\\ \cup \left\{\left(\frac{d-\rho \left(m-d\right)}{2d},y\right):y\in \left[\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s,P\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)\right]\right\}\\ \cup \left\{\left(x,P\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)\right):x\in \left[-\sqrt{m\cdot P\left(\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}-s\right)},\frac{d-\rho \left(m-d\right)}{2d}\right]\right\},\end{array}$

for $s\in \left[\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}\right]$. See Figure 3 for an illustration.

Remark 6.3. As $\epsilon \to 0$ , the large-amplitude limit cycle ${\Gamma }_{\epsilon }$ in theorem 4.1 converges to $\Gamma \left(s\right)$ at $s=\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}$ in the Hausdorff distance.

Figure 3. Singular canard cycles $\Gamma \left(s\right)$. Left: $\Gamma \left(s\right)$ for $s\in \left[\mathrm{0,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right]$. Right: $\Gamma \left(s\right)$ for $s\in \left[\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}\right]$.

Figure 4. Blow-up singular cycles for $\rho ={\rho }_0=\frac{d}{m+d}$.

By blowing up the transcritical point $\left(\frac{d-\rho \left(m-d\right)}{2d}\mathrm{,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{4m{d}^2}\right)$ of

system (2.2) and the canard point (0, 0) of system (2.2), see Figure 4. The following results can be obtained by theorem 3.3 in [18].

Theorem 6.4. Fix ${\epsilon }_0$ sufficiently small. Then for $\epsilon \in \left(\mathrm{0,}{\epsilon }_0\right)$ , system (2.2) possesses a family of periodic orbits

$s\to \left({\rho }_0+\rho \left(s\mathrm{,}\sqrt{\epsilon }\right)\mathrm{,}\Gamma \left(s\mathrm{,}\sqrt{\epsilon }\right)\right)\mathrm{,}s\in \left(\mathrm{0,}\frac{{\left[d-\rho \left(m-d\right)\right]}^2}{2m{d}^2}\right)\mathrm{,}$

whichis smooth in $\left(s\mathrm{,}\sqrt{\epsilon }\right)$ ,and such that:

1) As $\epsilon \to 0$ ,the family $\Gamma \left(s\mathrm{,}\sqrt{\epsilon }\right)$ converges uniformly in Hausdorff distance to $\Gamma \left(s\right)$ .

2) Any periodic orbit passing sufficiently close to the slow manifold is a member of the family $\Gamma \left(s\mathrm{,}\sqrt{\epsilon }\right)$ or a relaxation oscillation.

3) All canard cycles are stable and the function $\rho \left(s\mathrm{,}\sqrt{\epsilon }\right)$ ismonotonic in s.

7. Conclusions

As shown in this paper, canard explosion phenomenon in the predator-prey model with Michaelis-Menten functional response happens due to the interactions between the local dynamics nearby turning points, such as, fold point, transcritical point, canard point, and the global return mechanism induced by the slow manifold in system (1.6). Additionally, canard explosion phenomenon in two-dimensional singularly perturbed autonomous dynamical system is a codimension one bifurcating phenomenon, in which the parameter $\rho$ is chosen as

a bifurcating parameter, as the parameter $\rho$ decreases through $\rho =\frac{d}{m+d}$, the sudden transition from a small-amplitude periodic solution, which bifurcates from the equilibrium $\left(\rho \mathrm{,}\frac{\rho \left(1-\rho \right)}{d}\right)$ via the supercritical Hopf bifurcation at $\rho =\frac{d}{m+d}$, to a large-amplitude relaxation oscillation which emerges at $\rho <\frac{d}{m+d}$, takes place by canard explosions.

However, the global return mechanism in the predator-prey model is slightly different from that in van der Pol’ equations analyzed by Krupa and Szmolyan [18]; the latter is S shape; the former is not S shape.

Additionally, canard explosion phenomenon in two dimensional singularly perturbed autonomous dynamical system is a codimension one bifurcating phenomenon. In this paper, the parameter $\rho$ is selected as a bifurcating parameter, and canard explosion phenomenon in system (1.6) is demonstrated. Actually, the parameter $\beta$ in system (1.5) can be also chosen as a bifurcating parameter, and it can be shown that canard explosion phenomenon happens in

system (1.5) as the parameter $\beta$ decreases through $\beta =\frac{m-d}{m+d}$.

Acknowledgements

This work was supported by the NNSFC 11971477.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Rosenzweig, M.L. and MacArthur, R.H. (1963) Graphical Representation and Stability Conditions of Predator-Prey Interactions. The American Naturalist, 97, 209-223. https://doi.org/10.1086/282272
[2]	Cheng, K.S. (1981) Uniqueness of a Limit Cycle for a Predator-Prey System. SIAM Journal on Mathematical Analysis, 12, 541-548. https://doi.org/10.1137/0512047
[3]	Hsu, S.B., Hubbell, S.P. and Waltman, P. (1978) Competing Predators. SIAM Journal on Applied Mathematics, 35, 617-625. https://doi.org/10.1137/0135051
[4]	Hsu, S.B., Hubbell, S.P. and Waltman, P. (1978) A Contribution to the Theory of Competing Predators. Ecological Monographs, 48, 337-349. https://doi.org/10.2307/2937235
[5]	Wolkowicz, G.S.K. (1988) Bifurcation Analysis of a Predator-Prey System Involving Group Defence. SIAM Journal on Applied Mathematics, 48, 592-606. https://doi.org/10.1137/0148033
[6]	Ruan, S. and Xiao, D. (2001) Global Analysis in a Predator-Prey System with Nonmonotonic Functional Response. SIAM Journal on Applied Mathematics, 61, 1445-1472. https://doi.org/10.1137/S0036139999361896
[7]	Zhu, H., Campbell, S.A. and Wolkowicz, G.S.K. (2002) Bifurcation Analysis of a Predator-Prey System with Nonmonotonic Functional Response. SIAM Journal on Applied Mathematics, 63, 636-682. https://doi.org/10.1137/S0036139901397285
[8]	Cantrell, R.S., Cosner, C. and Ruan, S. (2004) Intraspecific Interference and Consumer-Resource Dynamics. Discrete and Continuous Dynamical Systems B, 4, 527-546. https://doi.org/10.3934/dcdsb.2004.4.527
[9]	Abbakar, K., Yang, Y., Mohamed, A., Xia, S., Abkar, M. and Hassan, O. (2021) Stability Analysis and Hopf Bifurcation for ODE System of Predator-Prey Model with Mutual Interference. Applied Mathematics, 12, 793-802. https://doi.org/10.4236/am.2021.129053
[10]	Lu, M. and Huang, J. (2021) Global Analysis in Bazykins Model with Holling II Functional Response and Predator Competition. Journal of Differential Equations, 280, 99-138. https://doi.org/10.1016/j.jde.2021.01.025
[11]	Benoit, E., Callot, J.L., Diener, F. and Diener, M. (1981) Chasse au canards. Collectanea Mathematica, 31-32, 37-119.
[12]	Eckhaus, W. (1983) Relaxation Oscillations Including a Standard Chase on French Ducks. In: Verhulst, F., Ed., Asymptotic Analysis II, Springer, Berlin, 449-497. https://doi.org/10.1007/BFb0062381
[13]	Dumortier, F. and Roussarie, R. (1996) Canard Cycles and Center Manifolds. Memoirs of the American Mathematical Society, No. 577. https://doi.org/10.1090/memo/0577
[14]	Fenichel, N. (1979) Geometric Singular Perturbation Theory for Ordinary Differential Equations. Journal of Differential Equations, 31, 53-98. https://doi.org/10.1016/0022-0396(79)90152-9
[15]	Jones, C.K.R.T. (1995) Geometric Singular Perturbation Theory. In: Johnson, R., Ed., Dynamical Systems, Springer, Berlin, 44-118. https://doi.org/10.1007/BFb0095239
[16]	Krupa, M. and Szmolyan, P. (2001) Extending Geometric Singular Perturbation Theory to Nonhyperbolic Points—Fold and Canard Points in Two Dimensions. SIAM Journal on Mathematical Analysis, 33, 286-314. https://doi.org/10.1137/S0036141099360919
[17]	Krupa, M. and Szmolyan, P. (2001) Extending Slow Manifolds near Transcritical and Pitchfork Singularities. Nonlinearity, 14, 1473-1491. https://doi.org/10.1088/0951-7715/14/6/304
[18]	Krupa, M. and Szmolyan, P. (2001) Relaxation Oscillation and Canard Explosion. Journal of Differential Equations, 174, 312-368. https://doi.org/10.1006/jdeq.2000.3929
[19]	Brons, M. and Kaasen, R. (2010) Canards and Mixed-Mode Oscillations in a Forest Pest Model. Theoretical Population Biology, 77, 238-242. https://doi.org/10.1016/j.tpb.2010.02.003
[20]	Krupa, M. and Wechselberger, M. (2010) Local Analysis near a Folded Saddle-Node Singularity. Journal of Differential Equations, 248, 2841-2888. https://doi.org/10.1016/j.jde.2010.02.006
[21]	Li, C. and Zhu, H. (2013) Canard Cycles for Predator-Prey Systems with Holling Types of Functional Response. Journal of Differential Equations, 254, 879-910. https://doi.org/10.1016/j.jde.2012.10.003
[22]	Vo, T. and Wechselberger, M. (2015) Canards of Folded Saddle-Node Type I. SIAM Journal on Mathematical Analysis, 47, 3235-3283. https://doi.org/10.1137/140965818
[23]	Jardon-Kojakhmetov, H., Broer, H.W. and Roussarie, R. (2016) Analysis of a Slow-Fast System near a Cusp Singularity. Journal of Differential Equations, 260, 3785-3843. https://doi.org/10.1016/j.jde.2015.10.045
[24]	Mitry, J. and Wechselberger, M. (2017) Folded Saddles and Faux Canards. SIAM Journal on Applied Dynamical Systems, 16, 546-596. https://doi.org/10.1137/15M1045065
[25]	Kuehn, C. (2015) Multiple Time Scale Dynamics. Volume 191, Springer. https://doi.org/10.1007/978-3-319-12316-5