# Relaxation Oscillations and Dynamical Properties in Two Time-Delay Slow-Fast Modified Leslie-Gower Models

Authors: Yufeng Wang; Youhua Qian; Bingwen Lin
Journal: Complexity
 (2020)

---

## Abstract

In this paper, we consider two kinds of time-delay slow-fast modified Leslie-Gower models. For the first system, we prove the existence and uniqueness of relaxation oscillation cycle through the geometric singular perturbation theory and entry-exit function. For the second system, we put forward a conjecture that the relaxation oscillation of the system is unique. Numerical simulation also verifies our results for the systems.


---

## Body

## 1. Introduction

In 2003, Alaoui and Okiye proposed the next modified Leslie–Gower model:(1)dxdT=r1−b1x−a1yx+k1x,dydT=r2−a2yx+k2y,where x is the density of a prey, y is the density of a predator, and the parameters r1,k1,a1,b1,r2,k2,anda2 are positive, which play biology roles.Consider the rescaling,(2)x=r1b1x¯,y=r1r2a2b1y¯,T=1r1t.Substitute (2) into (1), and we have x¯ and y¯ in terms of x and y, then we have(3)dxdt=x′=x1−x−axyx+e1,dydt=y′=εy1−yx+e2,where a=a1r2/a2r1,e1=b1k1/r1,e2=b1k2/r1, and ε=r2/r1. The parameter ε can be regard as a small parameter in some special models, such as prey hares and predator lynx.Consider the case whereε is small enough, system (3) could be treat as a slow-fast system which has the fast variable x and the slow one y. Applying the geometric singular perturbation theory, we could study the slow-fast systems [4, 5].The geometric singular perturbation theory is mathematical rigorous that can analyze dynamics of some slow-fast systems. The study of invariant manifolds in Fenichel’s theory [3] is the basis of geometric singular perturbation theory, and the theory guarantees that any compact normally hyperbolic submanifold So of the critical manifold Co could perturb a locally invariant manifold Sε (0 <  ε≪ 1), which is O(ε), close to So. In recent years, geometric singular perturbation theory contains a fairly broad class of geometric points used to study the slow-fast systems, see, for example, [6–12]. And, for the exchange lemma, see [13, 14].Recently, Wang and Zhang proved the existence and uniqueness of relaxation oscillation cycle of the slow-fast modified Leslie-Gower model [15]. Valery et al. studied the global dynamics in the Leslie-Gower model with the Allee Effect [16]. Du et al. considered two delays induced Hopf bifurcation and double Hopf bifurcation in a diffusive Leslie-Gower predator-prey system [17]. Karl et al. discussed the bifurcation of critical sets and relaxation oscillations in singular fast-slow systems [18]. Wang and Zhang considered the stability loss delay and smoothness of the return map in slow-fast systems [19]. Ambrosio et al. addressed the canard phenomenon in a slow-fast modified Leslie-Gower model [20]. Xia et al. discussed relaxation oscillations and the mechanism in a periodically excited vector field with pitchfork-Hopf bifurcation [21]. Atabaigi and Barati studied relaxation oscillations and canard explosion in a predator-prey system of Holling and Leslie types [22]. Ai and Sadhu considered the entry-exit theorem and relaxation oscillations in slow-fast planar systems [23].The rest of the paper is organized as follows: Section2 introduces the entry-exit function theory. Section 3 studies the first time-delay slow-fast modified Leslie-Gower model. Section 4 investigates another kind of time-delay slow-fast modified Leslie-Gower model.


## 2. Entry-Exit Function

We are going to consider a slow-fast vector field in the form of(4)dxdt=fx,y,εx,dydt=εgx,y,ε,where x,y∈ℝ2 are state-space variables and the parameter ε represents the ratio of time scales. We define the functions f and g satisfying(5)f0,y,0<0,fory>0,f0,y,0>0,fory<0,g0,y,0<0.For system (4), the y-axis consists of equilibria when ε=0 and Y+≔x,yx=0,y>0 attracting and Y−≔x,yx=0,y<0 repelling (Figure 1(a)). The Figure 1(b) is the case for ε>0, since the orbit of system (4) is attracted by the y-axis, the orbit starting at x0,y0 with y0>0 and x0>0 slightly moves toward the downside at the speed of ε. After the orbit passes positive x-axis, the y-axis begins repelling, then the orbit tends to move away from the negative y-axis at the point 0,p0y0, and then it intersects the line x=x0 at the point whose y-coordinate saying pεy0 satisfies limε⟶0pεy0=p0y0, where p0y0 is determined by(6)∫y0p0y0f0,y,0g0,y,0dy=0.Figure 1
(a) Forε=0, the orbits of system (4) is shown and the y-axis consists of equilibria. (b) For 0<ε≪1, a symbolic orbit of system (4) starts at x0,y0 and ends at x=x0,y0<0.
(a)(b)


## 3. System with Two Positive Time Delays

Consider a system with time delaysτ1 and τ2, which are sufficiently small constants:(7)dxdt=x′=x1−x−ayx+e1xt−τ1,dydt=y′=εy1−yt−τ2x+e2.Using the Taylor formula,(8)xt−τ1≈xt−x′tτ1,yt−τ2≈yt−y′tτ2.Substitute (8) into system (7), it yields(9)dxdt=x′=x1−x−axyx+e1x+e1x+e1−ayτ1≔xfx,y,ε,dydt=y′=εy1−yx+e2x+e2x+e2−εyτ2≔εgx,y,ε.


### 3.1. Equilibria

Lemma 1.
For system (9), we have the following result: system (9) and system (3) have the same equilibria. So, that system (9) has equilibria E10,0,E20,e2, and E31,0.
Moreover, ifae2−e1<0, system (9) has a unique positive equilibrium E∗x∗,y∗ with x∗=−a+e1−1+a+e1−12−4ae2−e1/2,y∗=x∗+e2.


### 3.2. Hopf Bifurcation

We consider the bifurcation of system (9) at the unique positive equilibrium E∗x∗,y∗ about parameter τ1. For convenience, we denote(10)Fx,y=x1−x−axyx+e1x+e1x+e1−ayτ1,Gx,y=εy1−yx+e2x+e2x+e2−εyτ2,so the linearized system of system (9) at the unique positive equilibrium E∗x∗,y∗ could be written as(11)x′y′=FxFyGxGyxy,with(12)Fx=−2x∗2−e1−1x∗x∗+e1−ay∗τ1,Fy=−ax∗x∗+e1−ay∗τ1,Gx=εx∗+e2x∗+e2−εy∗τ2,Gy=ε−21−ετ2.Then, the characteristic equation of system (11) is(13)λ2−Fx+Gyλ−GxFy=0.Lemma 2.
For equation (13), we have the following result. If τ1=1−ετ2/2−εae21+2−ε/1−ετ21−a−e1−2ax∗+e1/ae2, then equation (13) has a pair of pure imaginary roots.Proof.
The condition of equation (13) has a pair of pure imaginary roots:(14)Fx+Gy=0,GxFy<0.
Considerε,τ1,  andτ2 are sufficiently small constants. We can easily get Gx>0,Fy<0,so GxFy<0 holds.
According to the equationFx+Gy=0, we reach(15)−2x∗2−e1−1x∗x∗+e1−ay∗τ1+ε−21−ετ2=0.
Takingy∗=x∗+e2 into equation (15), then we get(16)2x∗2+e1−1x∗x∗+e1−ax∗+e2τ1=ε−21−ετ2,i.e.,(17)τ1=1−ετ22−εae21+2−ε1−ετ21−a−e1−2ax∗+e1ae2.
The proof is completed.
Letλτ=pτ±iqτ be the pair of pure imaginary roots of equation (13) that are satisfying pτ1=0, qτ1≠0. Then the following transversality condition holds.Lemma 3.
Ifa≠a+e1−12−4ae2−e1 holds, then λ'τ1≠0.Proof.
The condition of equation (13) has a pair of pure imaginary roots:(18)Redλdτ1=Fx+Gy2a=12adFxdτ1=x∗x∗+e1−ay∗τ122x∗+e1−1=x∗x∗+e1−ay∗τ12−a+a+e1−12−4ae2−e1≠0.
The proof is completed.
Denoteτ1′=1−ετ2/2−εae21+2−ε/1−ετ21−a−e1−2ax∗+e1/ae2. Combining Lemma 2 and Lemma 3, we conclude the following.Theorem 1.
For system (9), if τ1=τ1′ and a≠a+e1−12−4ae2−e1, then the equilibrium E∗x∗,y∗ is a Hopf bifurcation point and a limit cycle occurs.


### 3.3. Dynamics of Limit Systems

Consider the limit systems of system (9). Setting ε=0 in (9), we reach the fast subsystem:(19)dxdt=x′=x1−x−axyx+e1x+e1x+e1−ayτ1,dydt=y′=0.Consider the slow time scaleτ=εt and take the singular again, we get the slow subsystem:(20)0=x1−x−axyx+e1x+e1x+e1−ayτ1,dydτ=y˙=y1−yx+e2x+e2x+e2−εyτ2,which is a differential-algebraic equation on the critical set:(21)C0=x,yx=0ory=1−xx+e1a≔Fx.We restrict the parametersa,e1, ande2 to a subset(22)U=a,e1,e2a+e1−12−4ae2−e1<a<e1e2.Under this restriction,E∗x∗,y∗ and E20,e2 are located in the exclusion part of C0. Some studies have shown that the branch x,yy=Fx,x>0，y>0 of C0 has a unique generic fold point DxM,yM=1−e1/2,F1−e1/2. And, the two branches of C0 cross transversally at the transcritical point Ex,y=0,e1/a (Figure 2). Then, the critical set C0 is divided into four parts by points D and E:(23)C0a=x,y1−e12<x<1,y=Fx,C0+=x,yx=0,y>e1a,C0r=x,y0<x<1−e12,y=Fx,C0−=x,yx=0,0<y<e1a,where C0r and C0− are normally hyperbolic repelling and C0a and C0+ are normally hyperbolic attracting.Figure 2
Illustration for the relaxation oscillations of system (9). The critical manifold C0 contains repelling part C0r and attracting part C0a. The dashed black line is the predator isocline, the red dot is the fold point DxM,yM, the solid blue line is the slow-fast cycle γ0, the double and single arrows are fast and slow flows, and the solid red is the attracting relaxation oscillation γε.In fact, when the parameter is limited toU, E∗x∗,y∗ is on C0r and E20,e2 is on C0−.


### 3.4. Relaxation Oscillation

Lemma 4.
For system (9), there exists a unique y¯∗e2<y¯∗<e1/a such that,(24)∫y¯∗yMf0,y,0g0,y,0dy=0.Proof.
Considering the system (9), for y¯∈e2,e1/a,(25)Iy¯=∫y¯yMf0,y,0g0,y,0dy=∫y¯yM1−ay/e1y1−y/e2⋅e1e1−ayτ1dy.
So that(26)Iy¯⟶−∞asy¯⟶e2+.
Furthermore(27)Ie1a=∫e1/ayMf0,y,0g0,y,0dy=∫e1/ayM1−ay/e1y1−y/e2⋅e1e1−ayτ1dy>0.
Then, we combine (25) and (27) and conclude that there exist a unique y¯∗e2<y¯∗<e1/a such that,(28)∫y¯∗yMf0,y,0g0,y,0dy=0.
The proof is completed.
Let us definexr to be the x-coordinate of the intersection point of y=y¯∗ and C0a. Then, define a singular slow-fast cycle γ0, which contains two slow segments on C0a from xr,y¯∗ to xM,yM and on the positive y-axis from 0,yM to 0,y¯∗ and the two fast connections from xM,yM to 0,yM and 0,y¯∗ to xr,y¯∗, respectively. The next theorem explains the existence and uniqueness of the relaxation oscillation.Theorem 2.
For system (9), restrict the parameters to subset U and let V be a tubular neighborhood of γ0. Then, for each fixed ε>0 sufficiently small, system (9) has a unique limit cycle γε⊂V, which is strongly attracting. Moreover, the cycle γε is the unique limit cycle that converges to γ0 in the Hausdorff distance [3] as ε⟶0.Proof.
By Fenichel’s theory, the critical submanifoldC0a perturbs to a nearby slow manifold Cεa, which is Oε near C0a. By Theorem 2.1 of [9] on the analysis of a jump point, the slow manifold Cεa can be continued and passes the fold point DxM,yM, and then it jumps to another attracting branch C0+.
Letxl be a sufficiently small positive number. Define two vertical sections Δin and Δout as shown in Figure 2:(29)Δin≔xl,y∈Vy∈Iin,Δout≔xl,y∈Vy∈Iout,where Iin and Iout are closed intervals centered at yM and yout=p0yM, respectively. Through the flow of system (4), we define the transition map Π:Δin⟶Δin, which is a composition of the next two maps:(30)Π1:Δin⟶Δout,Π2:Δout⟶Δin.
Then,Π:Δin⟶Δin is given by the composition Π:=Π2∘Π1.
Now, we analyze the properties of these two mapsΠ1 and Π2.(a)
Analysis ofΠ1. Use the same proof as Lemma 4, for each y0∈Iin, we can define p0y0, with e2<p0y0<e1/a through the following formula:(31)∫p0y0y0f0,y,0g0,y,0dy=0.(b)
Analysis ofΠ2. Consider two orbits γε1 and γε2 starting on Δout, from Fenichel’s theory, γε1 and γε2 will be attracted to Cεa at the exponential rate Oe−1/ε. By Theorem 2.1 of [9], γε1 and γε2 pass by the generic fold point DxM,yM contracting exponentially toward each other. Then, they fly to Δin.
Therefore, according to the results of (a) and (b), we obtain that the transition mapΠ:Δin⟶Δin is a contraction at the exponential rate Oe−1/ε. Then, it follows from the contraction mapping theorem that Π has a unique fixed point in Δin, which must be stable. This fixed point provides a unique relaxation oscillation cycle γε⊂V of system (9) passing Δin for each 0<ε≪1.
According to Fenichel’s theory and Theorem2.1 of [9], we obtain that the relaxation oscillation cycle γε converges to the slow-fast cycle γ0 as ε⟶0 in the Hausdorff distance.
So, we can get that the relaxation oscillation cycleγε is a unique limit cycle of system (9) located in V for each 0<ε≪1.
The theorem is proved.


### 3.5. Numerical Simulation

First of all, we provide an example to illustrate Theorem2.Select the parameter valuesa,e1,e2,ε,τ1,τ2=1,0.2,0.022,0.01,0.2,0.2, and it is easy to state that the system (9) has a unique positive equilibrium point E∗0.33,0.36. Numerical simulation shows that it has a unique limit cycle. Figure 3 shows the phase portrait and the time series of the relaxation oscillation cycle.(a)
The orbit is shown in red in phase space(b)
Time series with the black dashed curve (prey) and red full curve (predator)Figure 3
Numerical simulation of system (9).
(a)(b)In addition, we compare the phase portrait and time series of system (9) with system (7). Figure 4 presents the phase portrait of system (9) and system (7). We can see that the two curves are very close to each other. Figure 5 shows the time series of system (9) and system (7).(a)
Time series for the relaxation oscillation of the predator with the red dashed curve (system (7)) and the blue full curve (system (9))(b)
Time series for the relaxation oscillation of the prey with the blue dashed curve (system (7)) and the red full curve (system (9)). Time series for the relaxation oscillation of prey with the blue dashed curve (system (4)) and the red full curve (system (6)).Figure 4
Phase space: the dashed blue curve shows the orbit of system (9) and the red full curve shows the orbit of system (7).Figure 5
Time series of system (9) and system (7).
(a)(b)


## 3.1. Equilibria

Lemma 1.
For system (9), we have the following result: system (9) and system (3) have the same equilibria. So, that system (9) has equilibria E10,0,E20,e2, and E31,0.
Moreover, ifae2−e1<0, system (9) has a unique positive equilibrium E∗x∗,y∗ with x∗=−a+e1−1+a+e1−12−4ae2−e1/2,y∗=x∗+e2.


## 3.2. Hopf Bifurcation

We consider the bifurcation of system (9) at the unique positive equilibrium E∗x∗,y∗ about parameter τ1. For convenience, we denote(10)Fx,y=x1−x−axyx+e1x+e1x+e1−ayτ1,Gx,y=εy1−yx+e2x+e2x+e2−εyτ2,so the linearized system of system (9) at the unique positive equilibrium E∗x∗,y∗ could be written as(11)x′y′=FxFyGxGyxy,with(12)Fx=−2x∗2−e1−1x∗x∗+e1−ay∗τ1,Fy=−ax∗x∗+e1−ay∗τ1,Gx=εx∗+e2x∗+e2−εy∗τ2,Gy=ε−21−ετ2.Then, the characteristic equation of system (11) is(13)λ2−Fx+Gyλ−GxFy=0.Lemma 2.
For equation (13), we have the following result. If τ1=1−ετ2/2−εae21+2−ε/1−ετ21−a−e1−2ax∗+e1/ae2, then equation (13) has a pair of pure imaginary roots.Proof.
The condition of equation (13) has a pair of pure imaginary roots:(14)Fx+Gy=0,GxFy<0.
Considerε,τ1,  andτ2 are sufficiently small constants. We can easily get Gx>0,Fy<0,so GxFy<0 holds.
According to the equationFx+Gy=0, we reach(15)−2x∗2−e1−1x∗x∗+e1−ay∗τ1+ε−21−ετ2=0.
Takingy∗=x∗+e2 into equation (15), then we get(16)2x∗2+e1−1x∗x∗+e1−ax∗+e2τ1=ε−21−ετ2,i.e.,(17)τ1=1−ετ22−εae21+2−ε1−ετ21−a−e1−2ax∗+e1ae2.
The proof is completed.
Letλτ=pτ±iqτ be the pair of pure imaginary roots of equation (13) that are satisfying pτ1=0, qτ1≠0. Then the following transversality condition holds.Lemma 3.
Ifa≠a+e1−12−4ae2−e1 holds, then λ'τ1≠0.Proof.
The condition of equation (13) has a pair of pure imaginary roots:(18)Redλdτ1=Fx+Gy2a=12adFxdτ1=x∗x∗+e1−ay∗τ122x∗+e1−1=x∗x∗+e1−ay∗τ12−a+a+e1−12−4ae2−e1≠0.
The proof is completed.
Denoteτ1′=1−ετ2/2−εae21+2−ε/1−ετ21−a−e1−2ax∗+e1/ae2. Combining Lemma 2 and Lemma 3, we conclude the following.Theorem 1.
For system (9), if τ1=τ1′ and a≠a+e1−12−4ae2−e1, then the equilibrium E∗x∗,y∗ is a Hopf bifurcation point and a limit cycle occurs.


## 3.3. Dynamics of Limit Systems

Consider the limit systems of system (9). Setting ε=0 in (9), we reach the fast subsystem:(19)dxdt=x′=x1−x−axyx+e1x+e1x+e1−ayτ1,dydt=y′=0.Consider the slow time scaleτ=εt and take the singular again, we get the slow subsystem:(20)0=x1−x−axyx+e1x+e1x+e1−ayτ1,dydτ=y˙=y1−yx+e2x+e2x+e2−εyτ2,which is a differential-algebraic equation on the critical set:(21)C0=x,yx=0ory=1−xx+e1a≔Fx.We restrict the parametersa,e1, ande2 to a subset(22)U=a,e1,e2a+e1−12−4ae2−e1<a<e1e2.Under this restriction,E∗x∗,y∗ and E20,e2 are located in the exclusion part of C0. Some studies have shown that the branch x,yy=Fx,x>0，y>0 of C0 has a unique generic fold point DxM,yM=1−e1/2,F1−e1/2. And, the two branches of C0 cross transversally at the transcritical point Ex,y=0,e1/a (Figure 2). Then, the critical set C0 is divided into four parts by points D and E:(23)C0a=x,y1−e12<x<1,y=Fx,C0+=x,yx=0,y>e1a,C0r=x,y0<x<1−e12,y=Fx,C0−=x,yx=0,0<y<e1a,where C0r and C0− are normally hyperbolic repelling and C0a and C0+ are normally hyperbolic attracting.Figure 2
Illustration for the relaxation oscillations of system (9). The critical manifold C0 contains repelling part C0r and attracting part C0a. The dashed black line is the predator isocline, the red dot is the fold point DxM,yM, the solid blue line is the slow-fast cycle γ0, the double and single arrows are fast and slow flows, and the solid red is the attracting relaxation oscillation γε.In fact, when the parameter is limited toU, E∗x∗,y∗ is on C0r and E20,e2 is on C0−.


## 3.4. Relaxation Oscillation

Lemma 4.
For system (9), there exists a unique y¯∗e2<y¯∗<e1/a such that,(24)∫y¯∗yMf0,y,0g0,y,0dy=0.Proof.
Considering the system (9), for y¯∈e2,e1/a,(25)Iy¯=∫y¯yMf0,y,0g0,y,0dy=∫y¯yM1−ay/e1y1−y/e2⋅e1e1−ayτ1dy.
So that(26)Iy¯⟶−∞asy¯⟶e2+.
Furthermore(27)Ie1a=∫e1/ayMf0,y,0g0,y,0dy=∫e1/ayM1−ay/e1y1−y/e2⋅e1e1−ayτ1dy>0.
Then, we combine (25) and (27) and conclude that there exist a unique y¯∗e2<y¯∗<e1/a such that,(28)∫y¯∗yMf0,y,0g0,y,0dy=0.
The proof is completed.
Let us definexr to be the x-coordinate of the intersection point of y=y¯∗ and C0a. Then, define a singular slow-fast cycle γ0, which contains two slow segments on C0a from xr,y¯∗ to xM,yM and on the positive y-axis from 0,yM to 0,y¯∗ and the two fast connections from xM,yM to 0,yM and 0,y¯∗ to xr,y¯∗, respectively. The next theorem explains the existence and uniqueness of the relaxation oscillation.Theorem 2.
For system (9), restrict the parameters to subset U and let V be a tubular neighborhood of γ0. Then, for each fixed ε>0 sufficiently small, system (9) has a unique limit cycle γε⊂V, which is strongly attracting. Moreover, the cycle γε is the unique limit cycle that converges to γ0 in the Hausdorff distance [3] as ε⟶0.Proof.
By Fenichel’s theory, the critical submanifoldC0a perturbs to a nearby slow manifold Cεa, which is Oε near C0a. By Theorem 2.1 of [9] on the analysis of a jump point, the slow manifold Cεa can be continued and passes the fold point DxM,yM, and then it jumps to another attracting branch C0+.
Letxl be a sufficiently small positive number. Define two vertical sections Δin and Δout as shown in Figure 2:(29)Δin≔xl,y∈Vy∈Iin,Δout≔xl,y∈Vy∈Iout,where Iin and Iout are closed intervals centered at yM and yout=p0yM, respectively. Through the flow of system (4), we define the transition map Π:Δin⟶Δin, which is a composition of the next two maps:(30)Π1:Δin⟶Δout,Π2:Δout⟶Δin.
Then,Π:Δin⟶Δin is given by the composition Π:=Π2∘Π1.
Now, we analyze the properties of these two mapsΠ1 and Π2.(a)
Analysis ofΠ1. Use the same proof as Lemma 4, for each y0∈Iin, we can define p0y0, with e2<p0y0<e1/a through the following formula:(31)∫p0y0y0f0,y,0g0,y,0dy=0.(b)
Analysis ofΠ2. Consider two orbits γε1 and γε2 starting on Δout, from Fenichel’s theory, γε1 and γε2 will be attracted to Cεa at the exponential rate Oe−1/ε. By Theorem 2.1 of [9], γε1 and γε2 pass by the generic fold point DxM,yM contracting exponentially toward each other. Then, they fly to Δin.
Therefore, according to the results of (a) and (b), we obtain that the transition mapΠ:Δin⟶Δin is a contraction at the exponential rate Oe−1/ε. Then, it follows from the contraction mapping theorem that Π has a unique fixed point in Δin, which must be stable. This fixed point provides a unique relaxation oscillation cycle γε⊂V of system (9) passing Δin for each 0<ε≪1.
According to Fenichel’s theory and Theorem2.1 of [9], we obtain that the relaxation oscillation cycle γε converges to the slow-fast cycle γ0 as ε⟶0 in the Hausdorff distance.
So, we can get that the relaxation oscillation cycleγε is a unique limit cycle of system (9) located in V for each 0<ε≪1.
The theorem is proved.


## 3.5. Numerical Simulation

First of all, we provide an example to illustrate Theorem2.Select the parameter valuesa,e1,e2,ε,τ1,τ2=1,0.2,0.022,0.01,0.2,0.2, and it is easy to state that the system (9) has a unique positive equilibrium point E∗0.33,0.36. Numerical simulation shows that it has a unique limit cycle. Figure 3 shows the phase portrait and the time series of the relaxation oscillation cycle.(a)
The orbit is shown in red in phase space(b)
Time series with the black dashed curve (prey) and red full curve (predator)Figure 3
Numerical simulation of system (9).
(a)(b)In addition, we compare the phase portrait and time series of system (9) with system (7). Figure 4 presents the phase portrait of system (9) and system (7). We can see that the two curves are very close to each other. Figure 5 shows the time series of system (9) and system (7).(a)
Time series for the relaxation oscillation of the predator with the red dashed curve (system (7)) and the blue full curve (system (9))(b)
Time series for the relaxation oscillation of the prey with the blue dashed curve (system (7)) and the red full curve (system (9)). Time series for the relaxation oscillation of prey with the blue dashed curve (system (4)) and the red full curve (system (6)).Figure 4
Phase space: the dashed blue curve shows the orbit of system (9) and the red full curve shows the orbit of system (7).Figure 5
Time series of system (9) and system (7).
(a)(b)


## 4. System with Two Time Delays

In this section, we consider time delaysτ1 and τ2 in system (3), then we have(32)dxdt=x′=x1−x−axx+e1yt−τ1,dydt=y′=εy1−yxt−τ2+e2,where τ1 stands for delayed maturation of the predator and τ2 stands for the time needed to digest the prey.


### 4.1. Equilibria and Hopf Bifurcation

System (32) has the same equilibria with system (9). So, system (32) has a unique positive equilibrium E∗x∗,y∗withx∗=−a+e1−1+a+e1−12−4ae2−e1/2,y∗=x∗+e2.In order to get the conditions for Hopf bifurcation, we reach the linearized system of system (32) at the equilibrium E∗x∗,y∗:(33)x′y′=Axy+Bxt−τ1yt−τ1+Cxt−τ2yt−τ2,where(34)A=1−2x∗−ae1x+e12yt−τ100ε1−2y∗xt−τ2+e2,B=0−ax∗x∗+e100,C=00εy∗2xt−τ2+e220.Substitutingy∗=x∗+e2 into (34) yields(35)A=1−2x∗−ae1x+e12yt−τ100ε1−2x∗+e2xt−τ2+e2,B=0−ax∗x∗+e100,C=00εx∗+e22xt−τ2+e220,i.e.,(36)A=x∗x∗+e12x∗+e1−100−ε,B=0−ax∗x∗+e100,C=00ε0.Denoteτ=τ1+τ2, the associated characteristic equation of system (33) can be obtained as follows:(37)λ2−x∗x∗+e12x∗+e1−1+ελ+εax∗x∗+e1e−λτ+x∗x∗+e12x∗+e1−1=0.Lemma 5.
For equation (37), if a+e1−12−4ae2−e1≥2a, then equation (37) ha no pairs of pure imaginary roots. If a+e1−12−4ae2−e1<2a, then equation (37) has a pair of pure imaginary roots.Proof.
Letλ=±ωi be the pair of pure imaginary roots of equation (37), then we take λ=ωi into equation (37) with e−λτ=e−ωiτ=cosωτ−isinωτ, it yields(38)−ω2−x∗x∗+e12x∗+e1−1+εωi,+εax∗x∗+e1cosωτ−isinωτ+x∗x∗+e12x∗+e1−1=0.
Separating the real and imaginary parts of equation (38), then we get(39)−ω2+εx∗x∗+e12x∗+e1−1+aεx∗x∗+e1cosωτ=0,ωx∗x∗+e12x∗+e1−1+ε−aεx∗x∗+e1sinωτ=0.
So that(40)cosωτ=ω2−εx∗/x∗+e12x∗+e1−1aεx∗/x∗+e1≜Q0,sinωτ=ωx∗/x∗+e12x∗+e1−1+εaεx∗/x∗+e1.
According tosin2ωτ+cos2ωτ=1, we have(41)ω4+Dω2+E=0,with(42)D=x∗x∗+e12x∗+e1−1+ε2−2εx∗x∗+e12x∗+e1−1=x∗x∗+e12x∗+e1−12+ε2>0,E=εx∗x∗+e12x∗+e1−12−aεx∗x∗+e12=εx∗x∗+e122x∗+e1−12−a2,=εx∗x∗+e122x∗+e1−1+a2x∗+e1−1−a.
Since(43)2x∗+e1−1+a=a+e1−12−4ae2−e1>0,2x∗+e1−1−a=a+e1−12−4ae2−e1−2a.
We can conclude that ifa+e1−12−4ae2−e1≥2a, thatE≥0, which implies that Equation (41) has no positive roots, so equation (37) has no pairs of pure imaginary roots. Besides, if a+e1−12−4ae2−e1<2a, then E<0 and equation (41) has a positive root:(44)ω=−D+D2−4E2,which implies equation (37) has a pair of pure imaginary roots when(45)τ0j=arccosQ0+2jπω,j=0,1,2,…
The proof is completed.
Denoteλτ=pτ±iqτ be the pair of pure imaginary roots of equation (37) that are satisfying pτ0j=0, qτ0j=ω.Lemma 6.
(transversality condition) Ifa+e1−12−4ae2−e1<2a holds, then p′τ0j>0, j=0,1,2,….Proof.
Differentiating both sides of equation (37) with respect to τ, we reach(46)2λdλdτ+x∗x∗+e12x∗+e1−1+εdλdτ−λεax∗x∗+e1e−λτ=0,then(47)dλdτ−1=2λ+x∗/x∗+e12x∗+e1−1+ελεax∗/x∗+e1e−λτ.
So(48)Redλdττ=τ0j−1=Re2λ+x∗/x∗+e12x∗+e1−1+ελεax∗/x∗+e1e−λττ=τ0j=2ωcosωτ+ax∗/x∗+e12x∗+e1−1+εsinωτωεax∗/x∗+e1,=2ω2+Dεax∗/x∗+e12>0.
To sum up, we getp′τ0j>0.
The proof is completed.Theorem 3.
For system (32), if τ1+τ2=τ0j, j=0,1,2,…, system (32) undergoes a Hopf bifurcation and τ0j are the Hopf bifurcation values.Proof.
Ifτ1+τ2=τ0j, j=0,1,2,… then the equation (37) has a pair of purely imaginary roots which satisfy the transversality condition. Thus system (32) undergoes a Hopf bifurcation and τ0j are the Hopf bifurcation values.
The theorem is proved.Remark 1.
If0<τ1+τ2<τ00, then E∗x∗,y∗ is locally asymptotically stable and if τ1+τ2>τ00, then E∗x∗,y∗ is unstable.


### 4.2. Analysis of Limit Systems

In order to get rid of the delayτ1 and τ2, we denote x=x1, xt−τ2=x2, y=y1yt−τ1=y2. Then, we get(49)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1yt−τ1−τ2,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2xt−τ1−τ2+e2,i.e.,(50)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1y1t−τ,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2x1t−τ+e2.Whenτ=0, system (50) goes to(51)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1y1,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2x1+e2.Consider the case whereε is small enough, system (51) could be viewed as a slow-fast system which has the fast variable X=x1x2 and the slow variable Y=y1y2. We have(52)fX,Y,ε=f1,f2=x11−x1−ax1x1+e1y2,x21−x2−ax2x2+e1y1.Then, we arrive at the critical set:(53)C0=x,y∈ℝ2×ℝ2fX,Y,0=0=C01∪C02∪C03∪C04.with(54)C01=x1=0，x2=0,C02=x1=0，y1=1a1−x2x2+e1,C03=y2=1a1−x1x1+e1,x2=0,C04=y2=1a1−x1x1+e1,y1=1a1−x2x2+e1.Since(55)DXf=1−2x1−ae1y2x1+e12001−2x2−ae1y1x2+e12,Then, the eigenvalues of matrix (55) are(56)μ1=1−2x1−ae1y2x1+e12,μ2=1−2x2−ae1y1x2+e12.Thus, we reach the attractive part and the repulsive part ofC01 are, respectively,(57)C01a=x1=0,x2=0,y1>e1a,y2>e1a,C01r=x1=0,x2=0,y1<e1a,y2<e1a.The attractive part and the repulsive part ofC02 are, respectively,(58)C02a=x1=0,x2>1−e12,y1=1a1−x2x2+e1,y2>e1a,C02r=x1=0,0<x2<1−e12,y1=1a1−x2x2+e1,y2<e1aThe attractive part and the repulsive part ofC03 are, respectively,(59)C03a=x1>1−e12,x2=0,y1>e1a,y2=1a1−x1x1+e1,C03r=0<x1<1−e12,x2=0,y1<e1a,y2=1a1−x1x1+e1.and the attractive part and the repulsive part of C04 are, respectively,(60)C04a=x1>1−e12,x2>1−e12,y1=1a1−x2x2+e1,y2=1a1−x1x1+e1,C04r=0<x1<1−e12,0<x2<1−e12,y1=1a1−x2x2+e1,y2=1a1−x1x1+e1.


### 4.3. Conjecture of the System

In this section, we combine the result of Subsection4.1 and Subsection 4.2 and propose a conjecture.Conjecture 1.
Ifx∗>1−e1/2 and x∗/x∗+e12x∗+e1−1+ε=0, then the equilibrium E∗x∗,y∗is a Hopf bifurcation point of system (51) and a periodic solution is bifurcated from E∗x∗,y∗. Moreover, the relaxation oscillation cycle of the system (51) is unique.


### 4.4. Numerical Simulation

We show some numerical simulations in this part to verify our theoretical results.For system (32), we choose the parameter values a,e1,e2,ε=0.805,0.2,0.022,0.01, then we get x∗,y∗=0.422,0.444. A series of calculations gives us that τ00≈11.50. Thus, we choose τ1,τ2=3,1, and then the equilibrium E∗x∗,y∗ of system (32) is asymptotically stable (Figure 6). Besides, we choose τ1,τ2=3,13 and get the periodic solutions of system (32) (Figure 7).Figure 6
Chooseτ1,τ2=3,1 satisfying τ1+τ2=4<τ00, then the equilibrium E∗x∗,y∗ of system (32) is asymptotically stable. Time series with the black dashed curve (prey) and red full curve (predator).Figure 7
Chooseτ1,τ2=3,13 satisfying τ1+τ2=16>τ00, and then reach the periodic solutions of system (32) (time series with the black dashed curve (prey) and red full curve (predator)).For system (51), we choose the parameter values a,e1,e2,ε=1,0.2,0.022,0.01. Figure 8 shows the phase portrait and the time series of the relaxation oscillation cycle.Figure 8
Numerical simulation of system (51). (a) The orbit is shown in red in the phase space. (b) Time series with the black dashed curve (prey) and red full curve (predator).
(a)(b)


## 4.1. Equilibria and Hopf Bifurcation

System (32) has the same equilibria with system (9). So, system (32) has a unique positive equilibrium E∗x∗,y∗withx∗=−a+e1−1+a+e1−12−4ae2−e1/2,y∗=x∗+e2.In order to get the conditions for Hopf bifurcation, we reach the linearized system of system (32) at the equilibrium E∗x∗,y∗:(33)x′y′=Axy+Bxt−τ1yt−τ1+Cxt−τ2yt−τ2,where(34)A=1−2x∗−ae1x+e12yt−τ100ε1−2y∗xt−τ2+e2,B=0−ax∗x∗+e100,C=00εy∗2xt−τ2+e220.Substitutingy∗=x∗+e2 into (34) yields(35)A=1−2x∗−ae1x+e12yt−τ100ε1−2x∗+e2xt−τ2+e2,B=0−ax∗x∗+e100,C=00εx∗+e22xt−τ2+e220,i.e.,(36)A=x∗x∗+e12x∗+e1−100−ε,B=0−ax∗x∗+e100,C=00ε0.Denoteτ=τ1+τ2, the associated characteristic equation of system (33) can be obtained as follows:(37)λ2−x∗x∗+e12x∗+e1−1+ελ+εax∗x∗+e1e−λτ+x∗x∗+e12x∗+e1−1=0.Lemma 5.
For equation (37), if a+e1−12−4ae2−e1≥2a, then equation (37) ha no pairs of pure imaginary roots. If a+e1−12−4ae2−e1<2a, then equation (37) has a pair of pure imaginary roots.Proof.
Letλ=±ωi be the pair of pure imaginary roots of equation (37), then we take λ=ωi into equation (37) with e−λτ=e−ωiτ=cosωτ−isinωτ, it yields(38)−ω2−x∗x∗+e12x∗+e1−1+εωi,+εax∗x∗+e1cosωτ−isinωτ+x∗x∗+e12x∗+e1−1=0.
Separating the real and imaginary parts of equation (38), then we get(39)−ω2+εx∗x∗+e12x∗+e1−1+aεx∗x∗+e1cosωτ=0,ωx∗x∗+e12x∗+e1−1+ε−aεx∗x∗+e1sinωτ=0.
So that(40)cosωτ=ω2−εx∗/x∗+e12x∗+e1−1aεx∗/x∗+e1≜Q0,sinωτ=ωx∗/x∗+e12x∗+e1−1+εaεx∗/x∗+e1.
According tosin2ωτ+cos2ωτ=1, we have(41)ω4+Dω2+E=0,with(42)D=x∗x∗+e12x∗+e1−1+ε2−2εx∗x∗+e12x∗+e1−1=x∗x∗+e12x∗+e1−12+ε2>0,E=εx∗x∗+e12x∗+e1−12−aεx∗x∗+e12=εx∗x∗+e122x∗+e1−12−a2,=εx∗x∗+e122x∗+e1−1+a2x∗+e1−1−a.
Since(43)2x∗+e1−1+a=a+e1−12−4ae2−e1>0,2x∗+e1−1−a=a+e1−12−4ae2−e1−2a.
We can conclude that ifa+e1−12−4ae2−e1≥2a, thatE≥0, which implies that Equation (41) has no positive roots, so equation (37) has no pairs of pure imaginary roots. Besides, if a+e1−12−4ae2−e1<2a, then E<0 and equation (41) has a positive root:(44)ω=−D+D2−4E2,which implies equation (37) has a pair of pure imaginary roots when(45)τ0j=arccosQ0+2jπω,j=0,1,2,…
The proof is completed.
Denoteλτ=pτ±iqτ be the pair of pure imaginary roots of equation (37) that are satisfying pτ0j=0, qτ0j=ω.Lemma 6.
(transversality condition) Ifa+e1−12−4ae2−e1<2a holds, then p′τ0j>0, j=0,1,2,….Proof.
Differentiating both sides of equation (37) with respect to τ, we reach(46)2λdλdτ+x∗x∗+e12x∗+e1−1+εdλdτ−λεax∗x∗+e1e−λτ=0,then(47)dλdτ−1=2λ+x∗/x∗+e12x∗+e1−1+ελεax∗/x∗+e1e−λτ.
So(48)Redλdττ=τ0j−1=Re2λ+x∗/x∗+e12x∗+e1−1+ελεax∗/x∗+e1e−λττ=τ0j=2ωcosωτ+ax∗/x∗+e12x∗+e1−1+εsinωτωεax∗/x∗+e1,=2ω2+Dεax∗/x∗+e12>0.
To sum up, we getp′τ0j>0.
The proof is completed.Theorem 3.
For system (32), if τ1+τ2=τ0j, j=0,1,2,…, system (32) undergoes a Hopf bifurcation and τ0j are the Hopf bifurcation values.Proof.
Ifτ1+τ2=τ0j, j=0,1,2,… then the equation (37) has a pair of purely imaginary roots which satisfy the transversality condition. Thus system (32) undergoes a Hopf bifurcation and τ0j are the Hopf bifurcation values.
The theorem is proved.Remark 1.
If0<τ1+τ2<τ00, then E∗x∗,y∗ is locally asymptotically stable and if τ1+τ2>τ00, then E∗x∗,y∗ is unstable.


## 4.2. Analysis of Limit Systems

In order to get rid of the delayτ1 and τ2, we denote x=x1, xt−τ2=x2, y=y1yt−τ1=y2. Then, we get(49)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1yt−τ1−τ2,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2xt−τ1−τ2+e2,i.e.,(50)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1y1t−τ,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2x1t−τ+e2.Whenτ=0, system (50) goes to(51)dx1dt=x11−x1−ax1x1+e1y2,dx2dt=x21−x2−ax2x2+e1y1,dy1dt=εy11−y1x2+e2,dy2dt=εy21−y2x1+e2.Consider the case whereε is small enough, system (51) could be viewed as a slow-fast system which has the fast variable X=x1x2 and the slow variable Y=y1y2. We have(52)fX,Y,ε=f1,f2=x11−x1−ax1x1+e1y2,x21−x2−ax2x2+e1y1.Then, we arrive at the critical set:(53)C0=x,y∈ℝ2×ℝ2fX,Y,0=0=C01∪C02∪C03∪C04.with(54)C01=x1=0，x2=0,C02=x1=0，y1=1a1−x2x2+e1,C03=y2=1a1−x1x1+e1,x2=0,C04=y2=1a1−x1x1+e1,y1=1a1−x2x2+e1.Since(55)DXf=1−2x1−ae1y2x1+e12001−2x2−ae1y1x2+e12,Then, the eigenvalues of matrix (55) are(56)μ1=1−2x1−ae1y2x1+e12,μ2=1−2x2−ae1y1x2+e12.Thus, we reach the attractive part and the repulsive part ofC01 are, respectively,(57)C01a=x1=0,x2=0,y1>e1a,y2>e1a,C01r=x1=0,x2=0,y1<e1a,y2<e1a.The attractive part and the repulsive part ofC02 are, respectively,(58)C02a=x1=0,x2>1−e12,y1=1a1−x2x2+e1,y2>e1a,C02r=x1=0,0<x2<1−e12,y1=1a1−x2x2+e1,y2<e1aThe attractive part and the repulsive part ofC03 are, respectively,(59)C03a=x1>1−e12,x2=0,y1>e1a,y2=1a1−x1x1+e1,C03r=0<x1<1−e12,x2=0,y1<e1a,y2=1a1−x1x1+e1.and the attractive part and the repulsive part of C04 are, respectively,(60)C04a=x1>1−e12,x2>1−e12,y1=1a1−x2x2+e1,y2=1a1−x1x1+e1,C04r=0<x1<1−e12,0<x2<1−e12,y1=1a1−x2x2+e1,y2=1a1−x1x1+e1.


## 4.3. Conjecture of the System

In this section, we combine the result of Subsection4.1 and Subsection 4.2 and propose a conjecture.Conjecture 1.
Ifx∗>1−e1/2 and x∗/x∗+e12x∗+e1−1+ε=0, then the equilibrium E∗x∗,y∗is a Hopf bifurcation point of system (51) and a periodic solution is bifurcated from E∗x∗,y∗. Moreover, the relaxation oscillation cycle of the system (51) is unique.


## 4.4. Numerical Simulation

We show some numerical simulations in this part to verify our theoretical results.For system (32), we choose the parameter values a,e1,e2,ε=0.805,0.2,0.022,0.01, then we get x∗,y∗=0.422,0.444. A series of calculations gives us that τ00≈11.50. Thus, we choose τ1,τ2=3,1, and then the equilibrium E∗x∗,y∗ of system (32) is asymptotically stable (Figure 6). Besides, we choose τ1,τ2=3,13 and get the periodic solutions of system (32) (Figure 7).Figure 6
Chooseτ1,τ2=3,1 satisfying τ1+τ2=4<τ00, then the equilibrium E∗x∗,y∗ of system (32) is asymptotically stable. Time series with the black dashed curve (prey) and red full curve (predator).Figure 7
Chooseτ1,τ2=3,13 satisfying τ1+τ2=16>τ00, and then reach the periodic solutions of system (32) (time series with the black dashed curve (prey) and red full curve (predator)).For system (51), we choose the parameter values a,e1,e2,ε=1,0.2,0.022,0.01. Figure 8 shows the phase portrait and the time series of the relaxation oscillation cycle.Figure 8
Numerical simulation of system (51). (a) The orbit is shown in red in the phase space. (b) Time series with the black dashed curve (prey) and red full curve (predator).
(a)(b)


## 5. Conclusion

In this paper, we mainly study two kinds of time-delay slow-fast modified Leslie-Gower models. For the first system, through the geometric singular perturbation theory and entry-exit function, we prove the existence and uniqueness of relaxation oscillation cycle. For another system, whenτ=τ1+τ2 crosses some critical values τ0j, the Hopf bifurcation occurs. Meanwhile, we put forward a conjecture that the relaxation oscillation cycle of the system is unique when τ=0. Numerical simulation verified our theoretical results and indicated that our method is effective. Our results show that delay affects the stability of the positive equilibrium and produces more complex dynamics than the model without delay.[1, 2]



---

*Source: 1351397-2020-09-25.xml*