
In this paper, we investigate the stability and bifurcation of a Leslie-Gower predator-prey model with a fear effect and nonlinear harvesting. We discuss the existence and stability of equilibria, and show that the unique equilibrium is a cusp of codimension three. Moreover, we show that saddle-node bifurcation and Bogdanov-Takens bifurcation can occur. Also, the system undergoes a degenerate Hopf bifurcation and has two limit cycles (i.e., the inner one is stable and the outer is unstable), which implies the bistable phenomenon. We conclude that the large amount of fear and prey harvesting are detrimental to the survival of the prey and predator.
Citation: Hongqiuxue Wu, Zhong Li, Mengxin He. Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting[J]. Mathematical Biosciences and Engineering, 2023, 20(10): 18592-18629. doi: 10.3934/mbe.2023825
[1] | Rongjie Yu, Hengguo Yu, Chuanjun Dai, Zengling Ma, Qi Wang, Min Zhao . Bifurcation analysis of Leslie-Gower predator-prey system with harvesting and fear effect. Mathematical Biosciences and Engineering, 2023, 20(10): 18267-18300. doi: 10.3934/mbe.2023812 |
[2] | Zhenliang Zhu, Yuming Chen, Zhong Li, Fengde Chen . Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey. Mathematical Biosciences and Engineering, 2023, 20(6): 10977-10999. doi: 10.3934/mbe.2023486 |
[3] | Mengyun Xing, Mengxin He, Zhong Li . Dynamics of a modified Leslie-Gower predator-prey model with double Allee effects. Mathematical Biosciences and Engineering, 2024, 21(1): 792-831. doi: 10.3934/mbe.2024034 |
[4] | Manoj K. Singh, Brajesh K. Singh, Poonam, Carlo Cattani . Under nonlinear prey-harvesting, effect of strong Allee effect on the dynamics of a modified Leslie-Gower predator-prey model. Mathematical Biosciences and Engineering, 2023, 20(6): 9625-9644. doi: 10.3934/mbe.2023422 |
[5] | Saheb Pal, Nikhil Pal, Sudip Samanta, Joydev Chattopadhyay . Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model. Mathematical Biosciences and Engineering, 2019, 16(5): 5146-5179. doi: 10.3934/mbe.2019258 |
[6] | Shuo Yao, Jingen Yang, Sanling Yuan . Bifurcation analysis in a modified Leslie-Gower predator-prey model with fear effect and multiple delays. Mathematical Biosciences and Engineering, 2024, 21(4): 5658-5685. doi: 10.3934/mbe.2024249 |
[7] | Sangeeta Kumari, Sidharth Menon, Abhirami K . Dynamical system of quokka population depicting Fennecaphobia by Vulpes vulpes. Mathematical Biosciences and Engineering, 2025, 22(6): 1342-1363. doi: 10.3934/mbe.2025050 |
[8] | Christian Cortés García . Bifurcations in a discontinuous Leslie-Gower model with harvesting and alternative food for predators and constant prey refuge at low density. Mathematical Biosciences and Engineering, 2022, 19(12): 14029-14055. doi: 10.3934/mbe.2022653 |
[9] | Yuhong Huo, Gourav Mandal, Lakshmi Narayan Guin, Santabrata Chakravarty, Renji Han . Allee effect-driven complexity in a spatiotemporal predator-prey system with fear factor. Mathematical Biosciences and Engineering, 2023, 20(10): 18820-18860. doi: 10.3934/mbe.2023834 |
[10] | Peng Feng . On a diffusive predator-prey model with nonlinear harvesting. Mathematical Biosciences and Engineering, 2014, 11(4): 807-821. doi: 10.3934/mbe.2014.11.807 |
In this paper, we investigate the stability and bifurcation of a Leslie-Gower predator-prey model with a fear effect and nonlinear harvesting. We discuss the existence and stability of equilibria, and show that the unique equilibrium is a cusp of codimension three. Moreover, we show that saddle-node bifurcation and Bogdanov-Takens bifurcation can occur. Also, the system undergoes a degenerate Hopf bifurcation and has two limit cycles (i.e., the inner one is stable and the outer is unstable), which implies the bistable phenomenon. We conclude that the large amount of fear and prey harvesting are detrimental to the survival of the prey and predator.
Fishing is a method used in the industry to acquire fish products from natural or artificial bodies of water. With the development of fisheries, fishing has become more common. We see that harvesting for economic gain is a relatively regular occurrence in nature and that it significantly affects both the ecological balance and system dynamics. It is crucial to develop biological resources at their maximum sustainable yield while preserving the survival of all interacting populations, both ecologically and economically. However, if a species is overharvested, it can lead to ecological problems, as some people may prioritize profit over protecting the environment. Thus, the authors of [1,2] built mathematical models to analyze these problems, whose dynamical behaviors have attracted the interest of many scholars. There are three forms of harvesting: 1) constant harvesting, h(x)=h; 2) linear harvesting, h(x)=qEx; 3) nonlinear harvesting, h(x)=qExm1E+m2x, which is also called Michaelis-Menten-type harvesting.
Leslie and Gower studied the predator-prey relationship between two species and developed the famous Leslie-Gower predator-prey model [3], which has been widely discussed [4,5,6]. For example, the Leslie-Gower predator-prey model with the Allee effect and a generalist predator was studied in [7], where the authors found that the system exhibits a multi-stability phenomenon and undergoes various bifurcations. Huang et al. [8] studied the Leslie-Gower-type predator-prey model with constant-yield harvesting, and they found that the dynamical behavior of the model is very sensitive to the constant yield harvest of the predator.
Gupta et al. [9] considered the following Leslie-Gower predator-prey model with Michaelis-Menten-type prey-harvesting:
˙x=rx(1−xK)−αxy−qExm1E+m2x,˙y=sy(1−ynx),if(x,y)≠(0,0),˙y=0,if(x,y)=(0,0), | (1.1) |
where x and y are the prey and predator population densities, respectively. They studied the stability and bifurcation (saddle-node bifurcation and Hopf bifurcation) of system (1.1). Also, the existence of bionomic equilibria and optimal singular control were investigated. Based on system (1.1), Gupta and Chandra [10] introduced the Holling type Ⅱ functional response and obtained the bistable situation. The model exhibits several local bifurcations (saddle-node, Hopf, homoclinic and Bogdanov-Takens) which are ecologically important. Considering the group defense and nonlinear harvesting in prey, Kumar and Kharbanda [11] obtained that the density of the predator increases as the harvest rate of the predator decreases. Caraballo Garrido et al. [12] investigated the predator prey model with nonlinear harvesting with both constant and distributed delay by varying parameters. Some scholars [13,14,15,16,17] have combined other functional responses and harvesting to obtain more complex dynamical behaviors.
As with direct killing, indirect killing also has a significant impact on the dynamic behaviors of the system. The studies mentioned above, however, solely take into account the predator's direct killing. Predation danger may drive the prey to engage in anti-predation behaviors, such as habitat modifications or foraging, which may lower the prey's birth rate. Hence, Wang et al. [18] incorporated the fear effect into the reproduction of prey animals and obtained the following prey-predator model:
˙x=r0xf(k,y)−dx−ax2−pxy,˙y=cpxy−my, | (1.2) |
where f(k,y)=11+k0y accounts for the cost of anti-predator defense due to fear. They studied a model with a linear functional response or Holling type Ⅱ functional response. It was found that the fear effect has no impact on the dynamic behaviors of model (1.2). However, the dynamic behavior of model (1.2) with Holling type Ⅱ functional responses can be affected by the fear effect. Chen et al. [19] considered the influence of the fear effect and Leslie-Gower function on the dynamic behavior of the predator-prey model, and they demonstrated that there are many types of bifurcation phenomena, including transcritical bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation. Zhang et al. [20] studied a delayed diffusive predator-prey model with spatial memory and a nonlocal fear effect, and they investigated the stability, Hopf bifurcation and Turing-Hopf bifurcation of the system. Many scholars [21,22,23,24,25] have studied the prey-predator model with a fear effect.
In this paper, we incorporate the fear effect into system (1.1) and obtain the following model:
˙x=x(r01+k0y−d0−ax−c0y)−q0Exm1E+m2x,˙y=s0y(1−yhx),if(x,y)≠(0,0),˙y=0,if(x,y)=(0,0), | (1.3) |
where r0 is the birth rate of the prey and d0 is the natural death rate of the prey. In the ecological sense, it is clear that r0>d0. a represents the intra-species competition, c0 is the maximum predation rate, s0 is the intrinsic growth rates of the predators, h is a measure of the quality of the prey as food for the predator, k0 is the fear parameter, q0 is the catchability coefficient, E is the effort applied to harvest the prey species and m1,m2 are suitable constants. For simplicity, letting
¯x=m2m1Ex,¯y=m2m1Ehy,¯t=am21E2m22xt,¯r=r0m2am1E,¯k=k0m1Ehm2,¯d=d0m2am1E,¯c=c0ha,¯q=q0m2am21E,¯s=s0m2am1E, |
for x,y be positive, and dropping the bars, system (1.3) becomes
˙x=x2(r1+ky−d−x−cy)−qx21+x,˙y=sy(x−y), | (1.4) |
where r>d, and r, k, d, c, q and s are positive constants.
The key aim of this study on prey-predator models is to discuss the impacts of prey fear and prey harvesting on system dynamics. The bifurcation phenomenon that distinguishes system (1.4) from system (1.1) deserves further discussion. In addition, by analyzing the observed bifurcation phenomena, we can elucidate the benefits and drawbacks of prey harvesting on both populations.
The structure of the article is as follows. In Section 2, we obtain the boundedness of solutions and analyze the dynamical behaviors of origin. In Section 3, we discuss the existence of boundary equilibria and positive equilibria. In Section 4, we analyze the stability of equilibria. In Section 5, we show that system (1.4) undergoes saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation. In Section 6, we have a summary of the article.
We show that the positive solutions of system (1.4) are ultimately bounded.
Theorem 2.1. All solutions of system (1.4) are bounded for all t≥0.
Proof. Since system (1.3) is equivalent to system (1.4), we now prove that the solution of system (1.3) is bounded. From the first equation of system (1.3), we have
˙x≤x(r0−d0−ax), |
for t≥0, which immediately implies that lim supt→∞x(t)≤r0−d0a≜M. Then from the second equation of system (1.3), it follows that
˙y≤s0y(1−yhM) |
for t≥0, that is, lim supt→∞y(t)≤hM. Hence, x(t) and y(t) are bounded. The proof is completed.
Denote
q0=(d+q)k+c+1. |
Next, we show the dynamic behaviors of the origin of system (1.4).
Lemma 2.1. The types of origin in system (1.4) are as follows:
1) if r−d<q or r−d=q≤q0, the origin of system (1.4) is a non-hyperbolic attractor;
2) if r−d>q or r−d=q>q0, the origin of system (1.4) is a non-hyperbolic repeller.
Proof. The Jacobian matrix of system (1.4) at the origin is degenerate; then, we apply the blow-up method to analyze the type of origin. Notice that when x=0, we have that ˙x=0 and ˙y=−sy2<0, which means that system (1.4) has the invariant line x=0. Using the horizontal blow-up
(x,y)=(u,uv)anddτ=udt, |
system (1.4) can be rewritten as
˙u=u(r1+kuv−d−u−cuv−q1+u),˙v=sv(1−v)−v(r1+kuv−d−u−cuv−q1+u). | (2.1) |
(ⅰ) The equilibria of system (2.1) in u=0 are A(0,0) and B(0,1−r−d−qs) when r−d<q+s. The Jacobian matrix at the equilibria A and B are, respectively
JA=(r−d−q00q+s−(r−d)) |
and
JB=(r−d−q0(s−r+d+q)Rs2r−d−(q+s)), |
where
R=−kr2+(ks−c+(d+q)k)r+(c−q+1)s+c(d+q). |
The eigenvalues of matrix JA are λJA1=r−d−q and λJA2=q+s−(r−d)>0, and the eigenvalues of matrix JB are λJB1=r−d−q and λJB2=r−d−(q+s)<0. If r−d<q, that is, if λJA1<0 and λJB1<0, A is a saddle and B is a stable node (see Figure 1(a)). If q<r−d<q+s, that is, if λJA1>0 and λJB1>0, A is an unstable node and B is a saddle (see Figure 1(c)).
If r−d=q, that is, if λJA1=0 and λJB1=0, both A and B are degenerate equilibria. First, we consider the degenerate equilibrium A. Taking the time variable
dτ=sdt, |
we can obtain the Taylor expansion of system (2.1) as follows:
{˙u=q−1su2−qsu3−dk+kq+csu2v+o(|u,v|2),˙v=v−q−1suv−v2+dk+kq+csuv2+qsu2v+o(|u,v|2). |
By Theorem 7.1 in [26], the degenerate equilibrium A is a saddle-node if q≠1. If q=1, we have that the equilibrium A is a degenerate saddle from −qs<0.
Next, for the degenerate equilibrium B, make the following transformation:
(u,v)=(X,Y+1). |
System (2.1) becomes
{˙X=−[(k−1)q+dk+c+1]X2+[(k2−1)q+dk2]X3+o(|X,Y|2),˙Y=[(k−1)q+dk+c+1]X−sY−[(k2−1)q+dk2]X2+[(2k−1)q+2dk+2c+1]XY−sY2+o(|X,Y|2). | (2.2) |
Using (X,Y)=(sX1,[(k−1)q+dk+c+1]X1+Y1) and dτ=−sdt, system (2.2) becomes
{˙X1=α20X21+α30X31+α21X21Y1+o(|X1,Y1|2),˙Y1=Y1+β20X21+β11X1Y1+β02Y21+o(|X1,Y1|2), | (2.3) |
where
α20=q0−q,α21=(d+q)k+c,α30=(k2−k)q2+(2dk2−k2s+2ck−dk−c+k+s)q+d2k2−dk2s+2cdk+c2+dk+c,β20=(−2k2+3k−1)q2+(−4dk2+k2s−4ck+3dk+3c−3k−s+2)q−2d2k2+dk2s−4cdk−2c2−3dk−3c−1,β11=1−q,β02=1. |
If α20>0 (or α20<0), that is, if q<q0 (or q>q0), B a saddle-node with a stable parabolic sector on the right (or left). If α20=0, which implies that k<1, we get that q=dk+c+11−k. Next, substituting q=dk+c+11−k into the coefficients of the X3 term of system (2.3), we have
α30=s[dk+(1+c)(1+k)]>0. |
Explicitly, from Theorem 7.1 in [26] we know that the degenerate equilibrium B is a stable node if q=q0.
In summary, when r−d=q≤1, A and B are as shown in Figure 1(a). When 1<r−d=q≤q0, A and B are as shown in Figure 1(b). When r−d=q>q0, A and B are as shown in Figure 1(c).
After a blow-down, the origin in system (1.4) is an attractor when r−d<q, r−d=q≤1 (see Figure 2(a)) or 1<r−d=q≤q0 (see Figure 2(b)). The origin is a repeller when r−d=q>q0 or q<r−d<q+s (see Figure 2(c)).
(ⅱ) When r−d=q+s, system (2.1) has only one equilibrium, A(0,0) at u=0, whose Jacobian matrix is
JA=(s000). |
Expanding system (2.1) in a Taylor series and taking a time variable dτ=sdt, we have
{˙u=u+q−1su2+o(|u,v|2),˙v=1−qsuv−v2+o(|u,v|2). |
The coefficient of v2 is −1<0; from Theorem 7.1 in [26], we know that the equilibrium A is a saddle-node (see Figure 1(d)). After a blow-down, we see that the origin is a repeller for system (1.4) (see Figure 2(d)).
(ⅲ) System (2.1) has a unique equilibrium A(0,0) when r−d>q+s. The Jacobian matrix at the equilibrium A is
JA=(r−d−q00d−r+q+s). |
Obviously, the two eigenvalues are r−d−q>0 and d−r+q+s<0. It shows that the equilibrium A is a saddle (see Figure 1(d)). Similarly, the origin is a repeller for system (1.4) (see Figure 2(d)). The proof is completed.
In this section, we will discuss the existence of the boundary equilibria and the positive equilibria of system (1.4).
First, we analyze the existence of boundary equilibria. When y=0, the first equation of system (1.4) can be simplified into
˙x=x2(r−d−x)−qx21+x. |
We have
f(x)=x2−(r−d−1)x−(r−d−q), |
and the discriminant of f(x) is as follows:
Δ1=4(q∗−q), |
where
q∗=(r−d+1)24. |
The two roots of f(x)=0 can be expressed as
x1=r−d−1−√Δ12,x2=r−d−1+√Δ12. |
If r−d>q, f(x)=0 has only one positive root x2. If r−d=q, f(x)=0 has only one positive root q−1 if q>1.
Assume that r−d<q. When r−d≤1, obviously, f(x)=0 has no positive roots. When r−d>1, obviously, q∗>r−d. If q>q∗, then f(x)=0 has no positive roots. If q=q∗, then f(x)=0 has only one positive root r−d−12. If q<q∗, then f(x)=0 has two positive roots x1 and x2.
To summarize, we have the following lemma.
Lemma 3.1. The following claims regarding the existence of the boundary equilibria of system (1.4) are true.
1) If r−d>q, system (1.4) has a unique boundary equilibrium E2(x2,0).
2) If 1<r−d=q, system (1.4) has a unique boundary equilibrium E2(q−1,0).
3) If 1<r−d<q, we have the following three cases:
(a) if 1<r−d<q∗<q, system (1.4) has no boundary equilibrium;
(b) if 1<r−d<q=q∗, system (1.4) has a unique boundary equilibrium ¯E(r−d−12,0);
(c) if 1<r−d<q<q∗, system (1.4) has two boundary equilibria E1(x1,0), E2(x2,0).
Next, we will discuss the positive equilibria E(x,y) of system (1.4). Letting ˙x=˙y=0 in system (1.4), we have
{r1+ky−d−x−cy−q1+x=0,x−y=0. |
We denote
F(x)=k(c+1)x3+((k+1)(c+1)+dk)x2+(q0−(r−d))x−(r−d−q) |
and
F′(x)=3k(c+1)x2+2((k+1)(c+1)+dk)x+q0−(r−d), |
where the discriminant of F′(x) is
Δ2=4((k+1)(c+1)+dk)2−12k(c+1)(q0−(r−d)). |
Define
x∗=−2((k+1)(c+1)+dk)+√Δ26k(1+c),y∗=x∗. |
From F(x)=0, we have
q=−(1+x)[k(c+1)x2+(dk+c+1)x+d−r]kx+1. | (3.1) |
The Jacobian matrix of system (1.4) at E(x,y) is
JE=(x2[q−(x+1)2](1+x)2−x2[rk+c(1+kx)2](1+kx)2sx−sx), |
and
DetJE=[x2[rk+c(1+kx)2](1+kx)2−x2[q−(x+1)2](1+x)2]sx, |
TrJE=x2[q−(x+1)2](1+x)2−sx. |
Substituting (3.1) into DetJE and F′(x), we have
DetJE=(x+s)(1+x)2xF′(x). | (3.2) |
When r−d>q, it is easy to get that the equation F(x)=0 has a unique positive root x4 (see Figure 3(a)).
When q0<r−d=q, Δ2>0. We obtain that F(x)=0 has only one positive root x4 (see Figure 3(b)). When q0≥r−d=q, we find that F(x)=0 has no positive roots.
When r−d<q, obviously, F(x)=0 has no positive roots if q0≥r−d. If q0<r−d, we obtain that F(x)=0 has no positive roots when F(x∗)>0 (see Figure 3(c)). When F(x∗)=0, F(x)=0 has only one positive root x∗ (see Figure 3(d)). When F(x∗)<0, F(x)=0 has two positive roots x3 and x4 (see Figure 3(e)).
To summarize, we have the following lemma.
Lemma 3.2. The following claims regarding the existence of the boundary equilibria of system (1.4) are true.
1) If r−d>q, system (1.4) has a unique positive equilibrium E4(x4,y4).
2) If q0<r−d=q, system (1.4) has a unique positive equilibrium E4(x4,y4).
3) If q0<r−d<q, we obtain the following results:
(a) if F(x∗)>0, system (1.4) has no positive equilibrium;
(b) if F(x∗)=0, system (1.4) has a unique positive equilibrium E∗(x∗,y∗);
(c) if F(x∗)<0, system (1.4) has two positive equilibria E3(x3,y3) and E4(x4,y4).
In this section, we will discuss the stability of the equilibria.
Theorem 4.1. 1) If r−d>q, the unique boundary equilibrium E2 is a saddle.
2) If 1<r−d=q, the unique boundary equilibrium E2 is a saddle.
3) If 1<r−d<q=q∗, the unique boundary equilibrium ¯E is a saddle-node.
4) If 1<r−d<q<q∗, E1 is unstable and E2 is a saddle.
Proof. 1) The Jacobian matrix of system (1.4) at E2 is
JE2=(−√Δ1x221+x2−(rk+c)x220sx2). |
Obviously, the equilibrium E2 is a saddle.
2) The Jacobian matrix of system (1.4) at E2(q−1,0) is
JE2(q−1,0)=(−(q−1)3q−(rk+c)(q−1)20s(q−1)), |
which implies that E2 is a saddle.
3) The Jacobian matrix of system (1.4) at ¯E is
J¯E=(0−(rk+c)(r−d−1)240s(r−d−1)2), |
which means that ¯E is a degenerate equilibrium. First, we transform ¯E to the origin by letting X=x−r−d−12,Y=y. Then, system (1.4) can be rewritten as
{˙X=a01Y+a20X2+a11XY+a02Y2+o(|X,Y|2),˙Y=b01Y+b11XY+b02Y2+o(|X,Y|2), | (4.1) |
where
a01=−(r−d−1)2(kr+c)4,a20=(r−d−1)22(−r+d−1),a11=−(r−d−1)(kr+c),a02=(r−d−1)2rk24,b01=s(r−d−1)2,b11=s,b02=−s. |
Next, applying the following transformation:
X=u+v,Y=−2s(kr+c)(r−d−1)v,dτ=r−d−12dt, |
system (4.1) becomes
{˙u=c20u2+c11uv+c02v2+o(|u,v|2),˙v=v+d11uv+d02v2+o(|u,v|2), |
where
c20=−r−d−1s(r−d+1),c11=2[(r−d−1)2−s(r−d+1)](r−d+1)(r−d−1)s,c02=k2r5+(−3dk2−2k2s+2ck−3k2)r4+h1r3+h2r2+h3r+h4(r−d+1)(kr+c)2(r−d−1)2s,d11=−2r−d−1,d02=−2[−kr2+(dk−c+k)r+cd+c−2s](r−d−1)2(kr+c),h1=3d2k2+4dk2s−2k2s2−6cdk−4cks+6dk2+c2−6ck+3k2,h2=(−d3−2sd2+4s2d−3d2−3d+2s−1)k2+(6cd2+8cds+12cd+4s2+6c)k−3c2d−2c2s−3c2,h3=(3d2+4sd+6d+3)c2+(−2d3k−4sd2k−6d2k−6dk+4sk+4s2−2k)c−2d2k2s2−4dks2+2k2s2+4ks2,h4=(−d3−2sd2−3d2−3d+2s−1)c2+(−4ds2+4s2)c. |
Since c20<0, we get that ¯E is a saddle-node from Theorem 7.1 in [26].
4) The Jacobian matrices of system (1.4) at E1(x1,0) and E2(x2,0), respectively, are
JE1=(x21√Δ11+x21−x21(rk+c)0sx1) |
and
JE2=(−x22√Δ11+x22−x22(rk+c)0sx2). |
This proves that E1 is unstable and E2 is a saddle point. The proof is completed.
If E3 exists, from F′(x3)<0 and (3.2), we obtain that DetJE3<0. That is, E3 is a saddle if it exists. Hence, in the following discussion, we only study the stability of E4.
Define
s∗=x4[q−(1+x4)2](1+x4)2. |
Theorem 4.2. When r−d>q or r−d=q>q0, system (1.4) has a positive equilibrium E4. In addition, the following statements are true.
1) If s∗≤0 or 0<s∗<s, E4 is locally asymptotically stable.
2) If s<s∗, E4 is an unstable node or focus.
3) If s=s∗>0, E4 is a center or weak focus.
Proof. The Jacobian matrix at the equilibrium E4(x4,y4) is
JE4=(x24[q−(1+x4)2](1+x4)2−x24[rk+c(1+kx4)2](1+kx4)2sx4−sx4). |
Notice that F′(x4)>0; so, from (3.2), we know that DetJE4>0.
Obviously,
TrJE4=x4(s∗−s). |
Hence, it is easy to see that E4 is locally asymptotically stable if s∗≤0 or 0<s∗<s, an unstable node or focus if s<s∗ and a center or weak focus if s=s∗>0. The proof is completed.
Define
q1=−18s2+52s+1+18√s(s+8)3. |
Theorem 4.3. If E4 is locally asymptotically stable and 0<q<q1, then E4 is globally asymptotically stable.
Proof. Noting that r−d>q or r−d=q>q0, and according to Lemmas 3.1 and 3.2, system (1.4) has a boundary equilibrium E2 and positive equilibrium E4. By Lemma 2.1 and Theorem 4.1, both the origin and E2 are unstable. We assume that E4 is locally asymptotically stable. Now, we want to show that there is no limit cycle around E4. Hence, taking the Dulac function Φ(x,y)=1x2y2, we have
∂(ΦF)∂x+∂(ΦG)∂y=−x3+(s+2)x2+(2s+1−q)x+sx(1+x)2y2, |
where F=x2(r1+ky−d−x−cy)−qx21+x and G=sy(x−y).
Define
H=x3+(s+2)x2+(2s+1−q)x+s. |
Obviously, H>0 for q≤1+2s. In what follows, we only consider that q>1+2s. By calculation, the discriminant of H is
Δ3=−q[4q2+(s2−20s−8)q−4(s−1)3]108. |
By calculation, we obtain that q1>1+2s and
Δ3|q=1+2s=s(1+2s)(2s2+11s+32)108>0. |
Thus, we have that Δ3>0 for 1+2s<q<q1, which means that H=0 has no positive roots. Then, H>0.
To sum up, H>0 when 0<q<q1. Then, we have
∂(ΦF)∂x+∂(ΦG)∂y<0, |
which implies that E4 is globally asymptotically stable. The proof is completed.
Remark 4.1. If q<min{1,r−d}, from Theorem 4.2, E4 is locally asymptotically stable. Therefore, by Theorem 4.3, E4 is globally asymptotically stable. That is, when the intrinsic growth rate of the prey is high and the catchability coefficient for prey is low, the prey and predator will reach a steady state. Hence, a small catchability coefficient for prey will not lead to the extinction of prey and predator.
Remark 4.2. Assume that r−d=q≤q0. From Lemma 2.1, the origin is a attractor. By Lemma 3.1 and Theorem 4.1, E2 is unstable if it exists. From Lemma 3.2, system (1.4) has no positive equilibrium. Note that system (1.4) has no limit cycle. Therefore, the origin is globally asymptotically stable. That is, when the intrinsic growth rate of the prey and the catchability coefficient for prey are low, the prey and predator will become extinct.
Lemma 4.1 ([27]). The system
{˙x=y+Ax2+Bxy+Cy2+o(|x,y|2),˙y=Dx2+Exy+Fy2+o(|x,y|2), |
is equivalent to the system
{˙x=y,˙y=Dx2+(E+2A)xy+o(|x,y|2) |
in some small neighborhood of (0,0) after changes to the coordinates.
Lemma 4.2 ([27]). The system given by
˙x=y,˙y=x2+a30x3+a40x4+y(a21x2+a31x3)+y2(a12x+a22x2)+o(|x,y|4), |
is equivalent to the system given by
˙x=y,˙y=x2+Gx3y+o(|x,y|4) |
by some nonsingular transformations in the neighborhood of (0,0), where G=a31−a30a21.
By computation, from F(x∗)=F′(x∗)=0, we can express k and r in terms of x∗, c, d, s and q, as follows:
k=q−(1+x∗)2(1+c)(2c+2)x3∗+(4c+d+4)x2∗+(2c+2d+2)x∗+d+q,r=[(c+1)x2∗+(c+d+1)x∗+d+q]2(2c+2)x3∗+(4c+d+4)x2∗+(2c+2d+2)x∗+d+q, | (4.2) |
where q>(1+x∗)2(1+c) because k and r are positive.
Notice that
q−(r−d)=x2∗[(2(x∗+1)(c+1)+d)q−(x∗+1)2(c+1)2](x∗+1)2(2cx∗+d+2x∗)+q,(r−d)−q0=x∗[((3x∗+4)(c+1)+2d)q+(x∗−2)(x∗+1)2(c+1)2](x∗+1)2(2cx∗+d+2x∗)+q; |
clearly, q−(r−d)>0 and (r−d)−q0>0 when q>(1+x∗)2(1+c). Thus, q0<r−d<q.
By computation, we have
(1+x∗)2(1+c)−(x∗+s)(1+x∗)2x∗=(1+x∗)2(cx∗−s)x∗. |
From the above discussions, we can obtain the following theorem.
Theorem 4.4. Assume that (4.2) and q>(1+x∗)2(1+c) hold; system (1.4) has a unique positive equilibrium E∗(x∗,y∗).
1) If s≤cx∗ and q>(1+x∗)2(1+c), or if s>cx∗ and q>(x∗+s)(1+x∗)2x∗, then E∗ is a saddle-node with an unstable parabolic sector.
2) If s>cx∗ and (1+x∗)2(1+c)<q<(x∗+s)(1+x∗)2x∗, then E∗ is a saddle-node with a stable parabolic sector.
Proof. Obviously, we have that DetJE∗=0 by (3.2). Then, the type of E∗ depends on the sign of TrJE∗, as follows:
TrJE∗=x2∗(1+x∗)2(q−(x∗+s)(1+x∗)2x∗). |
First, moving E∗(x∗,y∗) to the origin by the transformation (x,y)=(X+x∗,Y+y∗), it follows that system (1.4) becomes
{˙X=ˆa10X+ˆa01Y+ˆa20X2+ˆa11XY+ˆa02Y2+o(|X,Y|2),˙Y=ˆb10X+ˆb01Y+ˆb20X2+ˆb11XY+ˆb02Y2+o(|X,Y|2), | (4.3) |
where
ˆa10=x2∗[q−(1+x∗)2](1+x∗)2,ˆa01=−x2∗[q−(1+x∗)2](1+x∗)2,ˆa11=−2x∗[q−(1+x∗)2](1+x∗)2,ˆa20=x∗[q(2+x∗)−2(1+x∗)3](1+x∗)3,ˆa02=x2∗[(1+c)(1+x∗)2−q]2[(1+x∗)(cx∗+d+x∗)+q](1+x∗)3,ˆb10=sx∗,ˆb01=−sx∗,ˆb20=0,ˆb11=s,ˆb02=−s. |
The eigenvalues of the Jacobian matrix at point E∗ are λ1=0 and λ2=ˆa10+ˆb01. If q≠(x∗+s)(1+x∗)2x∗, then λ2≠0.
Next, taking the transformation
(XY)=(ˆa01ˆa10−ˆa10ˆb10)(uv), |
and, by introducing the new time variable
dτ=TrJE∗dt, |
system (4.3) is rewritten as
{˙ˆu=ˆc20ˆu2+ˆc11ˆuˆv+ˆc02ˆv2+o(|ˆu,ˆv|2),˙ˆv=ˆv+ˆd20ˆu2+ˆd11ˆuˆv+ˆd02ˆv2+o(|ˆu,ˆv|2), |
where
ˆc20=s(ˆa201ˆa20x∗−ˆa01ˆa10ˆa11x∗+ˆa02ˆa210x∗)(ˆa01sx∗+ˆa210)TrJE∗,ˆc11=s[(−s3x2∗+ˆa10s2x∗)ˆa01+ˆa02ˆa10s2x2∗+ˆa210ˆa11sx∗+ˆa310ˆa20](ˆa01sx∗+ˆa210)TrJE∗,ˆc02=sx∗[(ˆa20−s)ˆa210+(ˆa11sx∗+s2x∗)ˆa10+ˆa02s2x2∗](ˆa01sx∗+ˆa210)TrJE∗,ˆd20=ˆa10(ˆa201ˆa20−ˆa01ˆa10ˆa11+ˆa02ˆa210)(ˆa01sx∗+ˆa210)TrJE∗,ˆd11=ˆa201s2x∗+(ˆa10ˆa11sx∗+2ˆa10s2x∗+2ˆa210ˆa20−ˆa210s)ˆa01−2ˆa02ˆa210sx∗−ˆa310ˆa11(ˆa01sx∗+ˆa210)TrJE∗,ˆd02=−ˆa01s3x2∗+ˆa02ˆa10s2x2∗+ˆa01ˆa10s2x∗+ˆa210ˆa11sx∗+ˆa310ˆa20(ˆa01sx∗+ˆa210)TrJE∗. |
By a simple calculation, we get
ˆc20=sx4∗[q−(1+x∗)2]M(1+x∗)2TrJE∗[(1+c)x2∗+(c+d+1)x∗+d+q], |
where
M=−[(3x∗+2)(1+c)+d]q+(1+c)2(1+x∗)3. |
Note that q>(1+x∗)2(1+c); then, q>(1+x∗)2. By computation, we can obtain
M|q=(1+x∗)2(1+c)=−((2x∗+1)(c+1)+d)(1+x∗)2(1+c), |
which implies that M<0 for q>(1+x∗)2(1+c). Therefore, the sign of ˆc20 is determined by TrJE∗. Considering the time transformation, and by using Theorem 7.1 in [26], if s>cx∗ and (1+x∗)2(1+c)<q<(x∗+s)(1+x∗)2x∗, that is, if TrJE∗<0, then E∗ is a saddle-node with a stable parabolic sector (see Figure 4(a)). If s≤cx∗ and q>(1+x∗)2(1+c), or if s>cx∗ and q>(x∗+s)(1+x∗)2x∗, that is, if TrJE∗>0, then E∗ is a saddle-node with an unstable parabolic sector (see Figure 4(b)). The proof is completed.
From F(x∗)=F′(x∗)=TrJE∗=0, we can express k, r and q in terms of x∗, c, d and s, as follows:
k=s−cx∗(2c+2)x2∗+(d+1)x∗+s,r=[(c+2)x2∗+(d+s+1)x∗+s]2x∗[(2c+2)x2∗+(d+1)x∗+s],q=(x∗+s)(1+x∗)2x∗, | (4.4) |
where s>cx∗.
Theorem 4.5. Assume that (4.4) and s>cx∗ hold.
1) If one of the following conditions holds: (1.1) x∗≥1; (1.2) 0<x∗≤cc+2; (1.3) cc+2<x∗<1,s≠2x2∗1−x∗, E∗ is a cusp of codimension two.
2) If cc+2<x∗<1 and s=2x2∗1−x∗ hold, E∗ is a cusp of codimension three.
Proof. 1) Let X=x−x∗ and Y=y−y∗; then, system (1.4) can be rewritten as follows:
{˙X=sx∗X−sx∗Y+sx∗−x2∗+2sx∗+1X2−2sXY+x∗(cx∗−s)2(c+2)x2∗+(d+s+1)x∗+sY2+o(|X,Y|2),˙Y=sx∗X−sx∗Y+sXY−sY2+o(|X,Y|2). | (4.5) |
Applying the transformation (u,v)=(−1sx∗X,−X+Y), system (4.5) becomes
{˙u=v+e20u2+e11uv+e02v2+o(|u,v|2),˙v=f20u2+f11uv+f02v2+o(|u,v|2), | (4.6) |
where
e20=−sx3∗T1T2(1+x∗),e11=2[c2x3∗−s(3c+2)x2∗−s(d+1)x∗−s2]T2,e02=−(cx∗−s)2T2,f20=−s2x4∗T1T2(1+x∗),f11=sx∗[2c2x3∗−s(5c+2)x2∗+s(s−1−d)x∗−s2]T2,f02=−c2x3∗+s(2−c)x2∗+s(d+2s+1)x∗+s2T2,T1=−(3cx∗+2c+d+3x∗+2)s+c2x2∗+c2x∗−cx2∗−dx∗−2x2∗−x∗,T2=(c+2)x2∗+(d+s+1)x∗+s. |
By Lemma 4.1, system (4.6) is equivalent to the following:
{˙x=y,˙y=D20x2+D11xy+o(|x,y|2), |
where
D20=f20,D11=f11+2e20=sx∗(2x2∗−s(1−x∗))1+x∗. |
Substituting s=cx∗ into T1, we get
T1|s=cx∗=−x∗(c+1)(2cx∗+c+d+2x∗+1)<0. |
So, T1<0 for s>cx∗, that is, D20≠0.
Obviously, if x∗≥1, we have that D11>0, that is, E∗ is a cusp of codimension two by the result in [28] (see Figure 4(c)). When x∗<1, from D11=0, we have that s=2x2∗1−x∗. Noting that s>cx∗, we have
s−cx∗|s=2x2∗1−x∗=x∗((c+2)x∗−c). |
Therefore, if 0<x∗≤cc+2 or cc+2<x∗<1,s≠2x2∗1−x∗ holds and E∗ is a cusp of codimension two.
2) If cc+2<x∗<1 and s=2x2∗1−x∗ hold, that is, D11=0; we will show that E∗ is a cusp of codimension three. When cc+2<x∗<1,s=2x2∗1−x∗ and (4.4) reduces to the following:
k=−cx∗−c+2x∗2(c+1)x2∗+(−2c+d−3)x∗−d−1,r=(cx2∗+(−c+d−3)x∗−d−1)2[2(c+1)x2∗+(−2c+d−3)x∗−d−1](x∗−1),q=(1+x∗)31−x∗. | (4.7) |
Note that cc+2<x∗<1; then, k,r,q in (4.7) are positive.
Then, system (4.5) becomes
{˙x1=g10x1+g01y1+g20x21+g11x1y1+g02y21+g30x31+g21x21y1+g12x1y21+g03y31+g40x41+g22x21y21+g13x1y31+g04y41+o(|x1,y1|4),˙y1=h10x1+h01y1+h11x1y1+h02y21+o(|x1,y1|4), | (4.8) |
where
g10=−2x3∗x∗−1,g01=2x3∗x∗−1,g20=−3x2∗x∗−1,g11=4x2∗x∗−1,g21=2x∗x∗−1,g02=(cx∗−c+2x∗)2x2∗(cx2∗+(−c+d−3)x∗−d−1)(x∗−1),g30=−x2∗(x∗−1)(1+x∗),g12=2(cx∗−c+2x∗)2x∗(cx2∗+(−c+d−3)x∗−d−1)(x∗−1),g40=1(x∗−1)(1+x∗)2,g03=(cx∗−c+2x∗)3x2∗(cx2∗+(−c+d−3)x∗−d−1)2(x∗−1),h10=−2x3∗x∗−1,g22=(cx∗−c+2x∗)2(cx2∗+(−c+d−3)x∗−d−1)(x∗−1),h01=2x3∗x∗−1,g13=2(cx∗−c+2x∗)3x∗(cx2∗+(−c+d−3)x∗−d−1)2(x∗−1),h11=−2x2∗x∗−1,g04=(cx∗−c+2x∗)4x2∗(cx2∗+(−c+d−3)x∗−d−1)3(x∗−1),h02=2x2∗x∗−1. |
Let x2=y1 and y2=˙y1; then, system (4.8) becomes
{˙x2=y2,˙y2=i20x22+i02y22+i30x32+i21x22y2+i12x2y22+i03y32+i40x42+i31x32y2+i22x22y22+i13x2y32+i04y42+o(|x2,y2|4), | (4.9) |
where
i20=−2((c+4)(c+1)x∗−c2+d+1)x5∗(cx2∗+(−c+d−3)x∗−d−1)(x∗−1),i02=52x∗,i30=−2x4∗Q1(cx2∗+(−c+d−3)x∗−d−1)2(1+x∗)(x∗−1),i21=x∗Q2(cx2∗+(−c+d−3)x∗−d−1)(1+x∗)(x∗−1),i12=−7+4x∗2x2∗(1+x∗),i03=−x∗−14x4∗(1+x∗),i40=−2x3∗Q3(cx2∗+(−c+d−3)x∗−d−1)3(x∗−1)(1+x∗)2,i31=2Q4(cx2∗+(−c+d−3)x∗−d−1)2(1+x∗)2,i13=(x2∗+x∗+2)(x∗−1)2x6∗(1+x∗)2,i22=−Q52x3∗(cx2∗+(−c+d−3)x∗−d−1)(1+x∗)2,i04=−(x∗−1)28x9∗(1+x∗)2,Q1=(4c3+20c2+24c+8)x4∗+(−4c3+3c2d−8c2+16cd−24c+12d−20)x3∗+(−4c3−3c2d−24c2+2cd+2d2−70c−50)x2∗+(4c3−3c2d+9c2−18cd+d2−18c−22d−23)x∗+3c2d+3c2−3d2−6d−3,Q2=(2c2+9c+8)x3∗+(−2c2+3c+d+5)x2∗+(−2c2−12c+3d−13)x∗+2c2−4d−4, |
Q3=(7c4+39c3+72c2+56c+16)x7∗+(−7c4+9c3d−26c3+45c2d−63c2+60cd−68c+24d−24)x6∗+(−14c4−9c3d+3c2d2−87c3−15c2d+15cd2−234c2−30cd+12d2−229c−24d−68)x5∗+(14c4−18c3d−3c2d2+51c3−96c2d+3cd2+d3+51c2−198cd+3d2+87c−117d+73)x4∗+(7c4+18c3d−6c2d2+59c3+27c2d−30cd2+d3+189c2−24cd−33d2+342c−21d+205)x3∗+(−7c4+9c3d+6c2d2−27c3+63c2d−9cd2−2d3+9c2+150cd−18d2+159c+114d+130)x2∗+(−9c3d+3c2d2−9c3−18c2d+21cd2−3d3−21c2+42cd+27d2+21c+63d+33)x∗−3c2d2−6c2d+3d3−3c2+9d2+9d+3,Q4=(2c3+10c2+16c+8)x5∗+(c2d+7c2+4cd+20c+4d+12)x4∗−(4c3+12c2−4cd+12c−8d+8)x3∗+(−2c2d−10c2−32c+4d−28)x2∗+(2c3+4c2−8cd+2d2−8c−12d−14)x∗+c2d+c2−2d2−4d−2,Q5=(c2+4)x4∗+(−3c−4d+20)x3∗+(−2c2+6c−7d+41)x2∗+(−3c+10d+14)x∗+c2+d+1. |
Let x3=x2 and y3=(1−i02x2)y2; then, system (4.9) becomes
{˙x3=y3,˙y3=j20x23+j30x33+j21x23y3+j12x3y23+j03y33+j40x43+j31x33y3+j22x23y23+j13x3y33+j04y43+o(|x3,y3|4), | (4.10) |
where
j20=i20,j30=2i02i20+i30,j21=i21,j12=−i202+i12,j03=i03,j40=i202i20−2i02i30i40,j31=i21i02+i31,j22=−i302+i22,j13=i03i02+i13,j04=i04. |
To delete the y33-term, x3y33-term and y43-term in system (4.10), we do the following two transformations:
x3=x4+j032x24y4+j136x34y4+j042x24y24,y3=(1+j03x4y4+j132x24y4+j04x4y24)y4; |
x4=x5,y4=y5+12j03j20x45. |
Hence, system (4.10) becomes
{˙x5=y5,˙y5=m20x25+m30x35+m21x25y5+m12x5y25+m40x45+m31x35y5+m22x25y25+o(|x5,y5|4), | (4.11) |
where
m20=j20,m30=j30,m21=j21,m12=j12,m40=j40,m31=j31−3j20j03,m22=j22. |
Using cc+2<x∗<1, we have that m20=−2((c+4)(c+1)x∗−c2+d+1)x5∗(cx2∗+(−c+d−3)x∗−d−1)(x∗−1)<0. Letting
x6=−x5,y6=y5−√−m20,τ=√−m20t, |
system (4.11) becomes (still denoting τ by t)
{˙x6=y6,˙y6=x26+n30x36+n21x26y5+n12x6y26+n40x46+n31x36y6+n22x26y26+o(|x6,y6|4), | (4.12) |
where
n30=m30m20,n21=−m21√−m20,n12=m12,n40=m40m20,n31=m31√−m20,n22=−m22. |
By Lemma 4.2, system (4.12) is equivalent to the following system:
{˙X=Y,˙Y=X2+GX3Y+o(|X,Y|4), |
where
G=√2((c+4)(c+1)x∗−c2+d+1)x∗(cx2∗+(−c+d−3)x∗−d−1)(x∗−1)δ(x∗)4((c+4)(c+1)x∗−c2+d+1)2(cx2∗+(−c+d−3)x∗−d−1)2(1+x∗)2x3∗ |
and
δ(x∗)=9∑i=0Pixi∗; |
here, the coefficients of Pi,i=0,⋯,9 are given in Appendix A.
Using cc+2<x∗<1, the sign of G is determined by δ(x∗). By computation, we have that δ(cc+2)=(16(29c2+80c+56))(3c2+cd+5c+2d+2)4(c+2)9>0 and δ(1)=384c+128d+384>0. Using Lemma 3.1 in [29], the number of roots for δ(x∗) in cc+2<x∗<1 is equal to that of positive roots for
μ(x∗)=(1+x∗)9δ(cx∗+c+2(c+2)(1+x∗))=16(c+2)99∑i=0Mixi∗ |
in cc+2<x∗<1, and the coefficients of Mi,i=0,⋯,9 are given in Appendix A. Obviously, Mi,i=0,⋯,9 are positive. Hence, μ(x∗)>0 in cc+2<x∗<1, which implies that δ(x∗) has no positive zeros in cc+2<x∗<1. Then, δ(x∗)≠0, that is G≠0 in cc+2<x∗<1, which means that E∗ is a cusp of codimension three (see Figure 4(d)). The proof is completed.
Theorem 4.6. Assume that q0<r−d<q and F(x∗)<0; system (1.4) has two positive equilibria E3(x3,y3) and E4(x4,y4), where E3 is always a saddle point. Moveover,
1) if q<(x4+s)(1+x4)2x4, E4 is a stable node or focus;
2) if q>(x4+s)(1+x4)2x4, E4 is an unstable node or focus;
3) if q=(x4+s)(1+x4)2x4, E4 is a center or weak focus.
Proof. It is clear that F′(x3)<0 and F′(x4)>0 (see Figure 3(e)). Combining (3.2) again, we have that DetJE3<0 and DetJE4>0. Thus, E3 is a saddle point.
In what follows, we consider that
TrJE4=x24[q−(x4+1)2](1+x4)2−sx4. |
When TrJE4<0, E4 is a stable node or focus; when TrJE4>0, E4 is an unstable node or focus; when TrJE4=0, E4 is a center or weak focus. The proof is completed.
We will analyze the bifurcations of system (1.4) in this section, including saddle-node bifurcation, Hopf bifurcation and Bogdanov-Takens bifurcation.
From Lemma 3.1, when 1<r−d<q<q∗, system (1.4) has two boundary equilibria E1(x1,0) and E2(x2,0). However, when q=qSN=q∗, only the boundary equilibrium ¯E exists. Therefore, according to Sotomayor's theorem [28], system (1.4) will produce a saddle-node bifurcation at ¯E.
Theorem 5.1. Assume that 1<r−d<q, with q=qSN being the bifurcation parameter; then, system (1.4) will undergo saddle-node bifurcation at ¯E.
Proof. The eigenvalues of J¯E are λ1=0 and λ2=s(r−d−1)2. Denote the eigenvectors of J¯E and JT¯E as
V=( 1 0) |
and
W=( 1(rk+c)(r−d−1)2), |
respectively.
Denote
F(x,y)=(F1(x,y)F2(x,y))=(x2(r1+ky−d−x−cy)−qx21+xsy(x−y)). |
Then,
Fq(¯E;qSN)=( −(r−d−1)22(r−d+1) 0), |
D2F(¯E;qSN)(V,V)=( −(r−d−1)2r−d+1 0). |
It is easy to know that
WTFq(¯E;qSN)=−(r−d−1)22(r−d+1)≠0,WT[D2F(¯E;qSN)(V,V)]=−(r−d−1)2r−d+1≠0. |
This proves the transversality conditions, which means that system (1.4) will undergo saddle-node bifurcation at ¯E. The proof is completed.
Similarly, it follows from Lemma 3.2 that a saddle-node bifurcation occurs at the positive equilibrium E∗.
Theorem 5.2. Assume that q0<r−d<q and F(x∗)=0; system (1.4) will undergo saddle-node bifurcation at E∗.
It follows from Theorem 4.6 that, if q=(x4+s)(1+x4)2x4, then TrJE4=0. Noting that DetJE4>0, the Jacobian matrix of E4 has a pair of purely imaginary eigenvalues. Thus, system (1.4) may undergo Hopf bifurcation at E4.
For simplicity, similar to the analyses of Dai et al. [30] and Lu et al. [31], we prove the Hopf bifurcation. Letting
˜x=xx4,˜y=yy4,˜t=x24t,˜r=rx4,˜k=kx4,˜d=dx4,˜c=c,˜α=1x4,˜q=qx24,˜s=sx4, |
and by dropping the tilde, system (1.4) becomes
˙x=x2(r1+ky−d−x−cy)−qx2α+x,˙y=sy(x−y), | (5.1) |
where r>d and the other parameters are positive.
Clearly, ˜E4(1,1) is an equilibrium of system (5.1), which implies that
r=[(c+d+1)(α+1)+q](k+1)α+1. |
Define
˜q0=α[(k+1)(c+1)+dk](α+1)(1−αk),˜q1=(α+1)2[(2k+1)(c+1)+dk](1−αk). |
Assume that system (5.1) has another positive equilibrium ˜E3(˜x3,˜y3). By computation, ˜x3 satisfies the following equation:
(x−1)Φ(x)=0, |
where
Φ(x)=k(c+1)(α+1)x2+(α+1)[(c+1)(αk+k+1)+dk]x+(1−αk)(˜q0−q). |
Note that ˜x3<1 is a unique positive root of Φ(x), which implies that αk<1 and ~q0<q. Also, substituting x=1 into Φ(x), we have
Φ(1)=(1−αk)(˜q1−q)>0, |
that is, \alpha k < 1 and q < \tilde{q}_1 .
The Jacobian matrix of system (5.1) at \tilde{E_4} is
J_{\tilde{E}_4} = \left( \begin{array}{ccc} -1+ \frac{q}{(\alpha+1)^{2}}& \frac{\big[(-2 c-d-1) \alpha-2 c-d-q-1\big] k-c (\alpha+1)}{(k+1) (\alpha+1)}\\ s&-s\\ \end{array} \right), |
and
DetJ_{\tilde{E}_4} = \frac{s(1-\alpha k)(\tilde{q}_1- q)}{(\alpha+1)^{2}(k+1)},\,\,\,\,\, TrJ_{\tilde{E}_4} = \frac{q-\tilde{q}}{(\alpha+1)^{2}}, |
where
\tilde{q} = (s+1) (\alpha+1)^{2}. |
We have the following results.
Theorem 5.3. Assuming that \alpha k < 1 and \tilde{q_0} < q < \tilde{q_1} , system (5.1) has the equilibrium \tilde{E}_4(1, 1) . Moveover,
1) \tilde{E_4}(1, 1) is a stable hyperbolic node or focus if q < \tilde{q} ;
2) \tilde{E_4}(1, 1) is an unstable hyperbolic node or focus if q > \tilde{q} ;
3) \tilde{E_4}(1, 1) is a fine focus or center if q = \tilde{q} .
Now, we will study the Hopf bifurcation around \tilde{E_4} in system (5.1). Obviously, the transversality condition
{\frac{\mbox{d}TrJ_{\tilde{E_4}}}{\mbox{d}q}}\Big|_{q = \tilde{q}} = \frac{1}{(\alpha+1)^{2}}\neq 0 |
holds. Then, we can determine the stability of the limit cycle around \tilde{E_4} by calculating the first Lyapunov number. First, using the transformation (\tilde{x}, \tilde{y}) = (x-1, y-1) , the Taylor expansion of system (5.1) at the origin takes the following form:
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{\tilde{x}} = \tilde{a}_{10}\tilde{x}+\tilde{a}_{01}\tilde{y} +\tilde{a}_{20}\tilde{x}^2+\tilde{a}_{11}\tilde{x}\tilde{y}+\tilde{a}_{02}\tilde{y}^2+\tilde{a}_{30}\tilde{x}^3 +\tilde{a}_{21}\tilde{x}^2\tilde{y}+\tilde{a}_{12}\tilde{x}\tilde{y}^2 +\tilde{a}_{03}\tilde{y}^3+o(|\tilde{x},\tilde{y}|^2), \\ \dot{\tilde{y}} = \tilde{b}_{10}\tilde{x}+\tilde{b}_{01}\tilde{y} +\tilde{b}_{20}\tilde{x}^2+\tilde{b}_{11}\tilde{x}\tilde{y}+\tilde{b}_{02}\tilde{y}^2, \end{array} \right. \end{array} \end{equation} | (5.2) |
where
\begin{array}{lcl} \tilde{a}_{10} = s,\,\,\,\, \tilde{a}_{01} = \frac{\big[(-s-1) \alpha-2 c-d-s-2\big] k-c}{k+1},\,\,\,\, \tilde{a}_{20} = \frac{2 s \alpha+s-1}{\alpha+1},\,\,\,\,\\ \tilde{a}_{11} = \frac{\big[(-2 s-2) \alpha-4 c-2 d-2 s-4\big] k-2 c}{k+1},\,\,\,\, \tilde{a}_{02} = \frac{(s \alpha+\alpha+c+d+s+2) k^{2}}{(k+1)^{2}}\\ \tilde{a}_{30} = \frac{\alpha^{2} s-2 \alpha-1}{(\alpha+1)^{2}},\,\,\,\,\, \tilde{a}_{21} = \frac{\big[(-s-1) \alpha-2 c-d-s-2\big] k-c}{k+1},\,\,\,\,\,\\ \tilde{a}_{12} = \frac{2 (s \alpha+\alpha+c+d+s+2) k^{2}}{(k+1)^{2}},\,\,\,\, \tilde{a}_{03} = - \frac{(s \alpha+\alpha+c+d+s+2) k^{3}}{(k+1)^{3}},\,\,\,\,\\ \tilde{b}_{10} = s,\,\,\,\, \tilde{b}_{01} = -s,\,\,\,\, \tilde{b}_{20} = 0,\,\,\,\, \tilde{b}_{11} = s,\,\,\,\, \tilde{b}_{02} = -s.\\ \end{array} |
Next, using the transformation (\tilde{u}, \tilde{v}) = \left(-\tilde{x}, \frac{\tilde{a}_{10}\tilde{x}+\tilde{a}_{01}\tilde{y}}{\sqrt{D}}\right) , where D = \tilde{a}_{10}\tilde{b}_{01}-\tilde{a}_{01}\tilde{b}_{10} = DetJ_{\tilde{E_4}} > 0 , system (5.2) becomes
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{\tilde{u}} = -\sqrt{D}\tilde{v}+\tilde{c}_{20}\tilde{u}^2+\tilde{c}_{11}\tilde{u}\tilde{v}+\tilde{c}_{02}\tilde{v}^2+\tilde{c}_{30}\tilde{u}^3 +\tilde{c}_{21}\tilde{u}^2\tilde{v}+\tilde{c}_{12}\tilde{u}\tilde{v}^2+\tilde{c}_{03}\tilde{v}^3+o(|\tilde{u},\tilde{v}|^3),\\ \dot{\tilde{v}} = \sqrt{D}\tilde{u}+\tilde{d}_{20}\tilde{u}^2+\tilde{d}_{11}\tilde{u}\tilde{v}+\tilde{d}_{02}\tilde{v}^2+\tilde{d}_{30}\tilde{u}^3 +\tilde{d}_{21}\tilde{u}^2\tilde{v}+\tilde{d}_{12}\tilde{u}\tilde{v}^2+\tilde{d}_{03}\tilde{v}^3+o(|\tilde{u},\tilde{v}|^3), \end{array} \right. \end{array} \end{equation} | (5.3) |
where the coefficients are given in Appendix B.
According to the results of [26], the first-order Lyapunov number can be written as
\begin{array}{lcl} l_{1} = \frac{\gamma_1 k^{2}+\gamma_2k+\gamma_3}{8\big[(\alpha s+\alpha+2 c+d+s+2) k+c\big](\alpha+1)^{2} (k+1)D}, \end{array} |
where
\begin{array}{lcl} \gamma_1 = (s^{3}+2 s^{2}+s) \alpha^{4}+(4 c \,s^{2}+2 d \,s^{2}+4 c s+2 d s-2 s^{2}-6 s-4) \alpha^{3}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,+(4 c^{2} s+4 c d s-5 c \,s^{2}+d^{2} s-2 d \,s^{2}-2 s^{3}-16 c s-8 d s-15 s^{2}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,-16 c-8 d-29 s-17) \alpha^{2}+(-8 c^{2} s-8 c d s-6 c \,s^{2}-2 d^{2} s-3 d \,s^{2}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,-s^{3}-16 c^{2}-16 c d-32 c s-4 d^{2}-17 d s-9 s^{2}-36 c-18 d-24 s-18) \alpha\\ \,\,\,\,\,\,\,\,\,\,\,\,\,+c \,s^{2}-4 c^{2}-4 c d+2 c s-d^{2}+s^{2}-4 c-2 d+2 s,\\ \gamma_2 = (2 c \,s^{2}-s^{3}+2 c s-s^{2}) \alpha^{3}+(4 c^{2} s+2 c d s-5 c \,s^{2}-d \,s^{2}-8 c s+2 s^{2}-8 c+3 s) \alpha^{2}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,+(-8 c^{2} s-4 c d s-c \,s^{2}+d \,s^{2}+2 s^{3}-16 c^{2}-8 c d-9 c s+3 d s+10 s^{2}-18 c+12 s\\ \,\,\,\,\,\,\,\,\,\,\,\,\,+2) \alpha+3 c \,s^{2}+d \,s^{2}+s^{3}-4 c^{2}-2 c d+10 c s+4 d s+6 s^{2}+2 c+2 d+10 s+4,\\ \gamma_3 = (c^{2} s-c \,s^{2}) \alpha^{2}+(-2 c^{2} s+c \,s^{2}-4 c^{2}+3 c s) \alpha+c \,s^{2}-c^{2}+4 c s+2 c. \end{array} |
Thus, we can obtain the following theorem about the Hopf bifurcation.
Theorem 5.4. If \alpha k < 1 , \tilde{q}_0 < q < \tilde{q}_1 and q = \tilde q , then the following statements hold.
1) If l_1 > 0 , then system (5.1) undergoes subcritical Hopf bifurcation and an unstable limit cycle comes out around \tilde{E_4} .
2) If l_1 < 0 , then system (5.1) undergoes supercritical Hopf bifurcation and a stable limit cycle appears around \tilde{E_4} .
3) If l_1 = 0 , then system (5.1) undergoes a degenerate Hopf bifurcation and multiple limit cycles may appear around \tilde{E_4} .
By numerical simulation, we show the existence of limit cycles. Letting k = 0.1, \alpha = 1, d = 1, c = 1, s = 1, q = 8 and r = 7.7 , we have that l_{1} = 0.001984126984 . We perturb q to q = 8-0.005 ; then, there exists an unstable limit cycle around \tilde{E_4} (see Figure 5(a), (b)). On the other hand, letting k = 0.1, \alpha = 1, d = 1, c = 1, s = 0.7, q = 6.8 and r = 7.04 , we obtain that l_{1} = - 0.06095323795 . We perturb q to q = 6.8+0.03 ; then, there exists a stable limit cycle around \tilde{E_4} (see Figure 5(c), (d)).
Now, we give an example to illustrate the existence of two limit cycles. The parameters are given as follows:
d = 1, c = 1, s = 1, \alpha = \frac{1}{2}, k = \frac{7}{116}+\frac{3\sqrt{57}}{116}, r = \frac{369}{58}+\frac{9\sqrt{57}}{58}, q = \frac{9}{2}, |
where l_1 = 0 . We perturb k and q to k = \frac{7}{116}+\frac{3\sqrt{57}}{116}+0.03 and q = \frac{9}{2}+0.01 . Hence, system (5.1) undergoes a degenerate Hopf bifurcation and has two limit cycles (the inner one is stable and the outer is unstable) around \tilde{E_4} (Figure 5(e), (f)).
Remark 5.1. In Figure 5(a), (b), the origin is a stable node and the boundary equilibria are unstable. In addition, system (5.1) has two positive equilibria, where \tilde{E}_3 is a saddle point and \tilde{E}_4 is a stable point, and an unstable limit cycle appears around \tilde{E}_4 . The orbits of the phase portraits reveal that the prey and predator tend to coexistent in steady states only when the initial values of system (5.1) lie inside the unstable limit cycle; otherwise, the prey and predator become extinct.
In Figure 5(c), (d), in addition to the origin being stable, the other two boundary equilibria and the two positive equilibria are unstable. System (5.1) has a stable limit cycle that appears around \tilde{E}_4 . When the initial values lie to the right of the two stable invariant manifolds of the saddle, the prey and predator tend to coexist in periodic orbits. In addition, when the initial values lie to the left of the two stable invariant manifolds of the saddle, the prey and predator tend to go extinct.
Figure 5(e), (f) show that system (5.1) undergoes a degenerate Hopf bifurcation and has two limit cycles (the inner one is stable and the outer is unstable) around \tilde{E}_4 . Prey and predator will oscillate and coexist if the initial values lie inside of the unstable limit cycle, while the prey and predator will become extinct if the initial values lie outside of the unstable limit cycle.
From Theorem 4.5(1), the unique positive equilibrium E_{*} of system (1.4) is a cusp of codimension two, which means that a Bogdanov-Takens bifurcation of codimension two may occur. Hence, using q and s as the bifurcation parameters, system (1.4) becomes
\begin{equation} \begin{array}{rcl} \dot{x}& = & x^2\left( \frac{r}{1+ky}-d-x-cy\right)- \frac{(q+\lambda_1)x^2}{1+x}, \\ \dot{y}& = & (s+\lambda_2)y\left(x-y\right), \end{array} \end{equation} | (5.4) |
where \lambda = (\lambda_1, \lambda_2) is a parameter vector in a small neighborhood of the origin.
Theorem 5.5. Assuming that the conditions of Theorem 4.5 (1) hold, system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two around E_{*} .
Proof. First, by initiating the transformation x_1 = x-x_{*} and y_1 = y-y_{*} to move the positive equilibrium E_{*} to the origin, system (5.4) becomes
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_1 = g_{00}+ g_{10}x_1+g_{01}y_1+g_{20}x_1^2+g_{11}x_1y_1+g_{02}y_2^2+o(|x_1,y_1 |^2),\\ \dot{y}_1 = h_{00}+h_{10}x_1+h_{01}y_1+h_{20}x_1^2+h_{11}x_1y_1+h_{02}y_1^2+o(|x_1,y_1 |^2), \end{array} \right. \end{array} \end{equation} | (5.5) |
where
\begin{array}{lcl} g_{00} = - \frac{x_{*}^{2} \lambda_\mathrm{1}}{1+x_{*}},\,\,\,\, g_{10} = \frac{\big[s x_{*}^{2}+(2 s-\lambda_\mathrm{1}) x_{*}+s-2 \lambda_\mathrm{1}\big] x_{*}}{(1+x_{*})^{2}},\,\,\,\, g_{01} = s x_{*},\,\,\,\,\\ g_{20} = \frac{-x_{*}^{4}+(s-2) x_{*}^{3}+(4 s-1) x_{*}^{2}+5 s x_{*}+2 s-\lambda_\mathrm{1}}{(1+x_{*})^{3}}, g_{11} = -2 s,\,\,\,\,\\ g_{02} = \frac{(c x_{*}-s)^{2} x_{*}}{\big[(c+2) x_{*}^{2}+(d+s+1)x_{*}+s\big]},\,\,\,\, h_{00} = 0,\,\,\,\, h_{10} = (s+\lambda_\mathrm{2}) x_{*},\,\,\,\,\\ h_{01} = -(s+\lambda_\mathrm{2}) x_{*},\,\,\,\, h_{20} = 0,\,\,\,\, h_{11} = s+\lambda_\mathrm{2},\,\,\,\, h_{02} = -s-\lambda_\mathrm{2}. \end{array} |
Second, letting
\begin{array}{lcl} x_2 = y_1,\\ y_2 = h_{10}x_1+h_{01}y_1+h_{20}x_1^2+h_{11}x_1y_1+h_{02}y_1^2, \end{array} |
system (5.5) can be written as follows:
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_2 = y_2,\\ \dot{y}_2 = j_{00}+j_{10}x_2+j_{01}y_2+j_{20}x_2^2+j_{11}x_2y_2+j_{02}y_2^2+o(|x_2,y_2 |^2), \end{array} \right. \end{array} \end{equation} | (5.6) |
where
\begin{array}{lcl} j_{00} = g_{00}h_{00},\,\,\,\, j_{10} = g_{00}h_{11}+g_{01}h_{10}-g_{10}h_{01},\,\,\,\, j_{01} = g_{10}h_{01},\,\,\,\,\\ j_{20} = \frac{g_{01}h_{10}h_{11}+g_{02}h_{10}^{2}-g_{10}h_{02} h_{10}-g_{11} h_{01}h_{10}+g_{20}h_{01}^{2}}{h_{10}},\\ j_{11} = \frac{g_{11}h_{10}-2g_{20}h_{01}-h_{11}h_{01}+2h_{10} h_{02}}{h_{10}},\,\,\,\, j_{02} = \frac{g_{20}+h_{11}}{h_{10}}. \end{array} |
Taking a new time variable \tau with \mbox{d}t = (1-j_{02}x_2)\mbox{d}\tau and x_3 = x_2, y_3 = (1-j_{02}x_2)y_2 , system (5.6) becomes
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_3 = y_3,\\ \dot{y}_3 = k_{00}+k_{10}x_3+k_{01}y_3+k_{20}x_3^2+k_{11}x_3y_3+k_{02}y_3^2+o(|x_3,y_3 |^2), \end{array} \right. \end{array} \end{equation} | (5.7) |
where
\begin{array}{lcl} k_{00} = j_{00},\,\,\,\, k_{10} = -2j_{00}j_{02}+j_{10},\,\,\,\, k_{01} = j_{01},\,\,\,\,\\ k_{20} = j_{00}j_{02}^{2}-2j_{02}j_{10}+j_{20},\,\,\,\, k_{11} = -j_{11},\,\,\,\, k_{02} = \frac{g_{20}+h_{11}}{h_{10}}. \end{array} |
From the proof of Theorem 4.4, we have
k_{20}\Big|_{\lambda_1 = \lambda_2 = 0} = \frac{s x_{*}^{3}T_1} {\big[(c+2) x_{*}^{2}+(d+s+1)x_{*}+s\big](1+x_{*})} < 0, |
where T_1 is defined in Theorem 4.4. Letting
x_4 = x_3,\,\,\,\,\,\,\,\,y_4 = \frac{y_3}{\sqrt{-k}_{20}},\,\,\,\,\,\,\,\,\tau = \sqrt{-k}_{20}t, |
system (5.7) becomes
\begin{equation} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_4 = y_4,\\ \dot{y}_4 = m_{00}+m_{10}x_4+m_{01}y_4-x_4^2+m_{11}x_4y_4+o(|x_4,y_4 |^2), \end{array} \right. \end{array} \end{equation} | (5.8) |
where
\begin{array}{lcl} m_{00} = - \frac{k_{00}}{k_{20}},\,\,\,\,\,\,\,\, m_{10} = - \frac{k_{10}}{k_{20}},\,\,\,\,\,\,\,\, m_{01} = \frac{k_{01}}{\sqrt{-k_{20}}},\,\,\,\,\,\,\,\, m_{11} = \frac{k_{11}}{\sqrt{-k_{20}}}. \end{array} |
Next, letting x_5 = x_4- \frac{m_{10}}{2} and y_5 = y_4 , system (5.8) is equivalent to the following system:
\begin{equation*} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_5 = y_5,\\ \dot{y}_5 = n_{00}+m_{01}y_5-x_5^2+n_{11}x_5y_5+o(|x_5,y_5 |^2), \end{array} \right. \end{array} \end{equation*} |
where
\begin{array}{lcl} n_{00} = m_{00}+ \frac{m_{10}^{2}}{4},\,\,\,\,\,\,\,\, n_{01} = m_{01}+ \frac{m_{11}m_{10}}{2},\,\,\,\,\,\,\,\, n_{11} = m_{11}. \end{array} |
From the proof of Theorem 4.4, we obtain
n_{11}\Big|_{\lambda_1 = \lambda_2 = 0} = \sqrt{- \frac{(2 x_{*}^{2}+s x_{*}-s)^2\big[(c+2) x_{*}^{2}+(d+s+1)x_{*}+s\big]} {T_1(1+x_{*})sx_{*}^3}}\neq 0. |
Finally, letting
x_6 = -n_{11}^2x_5,\,\,\,\,\,\,y_6 = -n_{11}^3y_{5},\,\,\,\,\,\,\tau = - \frac{1}{n_{11}}t, |
we obtain the universal unfolding of system (5.4) as follows:
\begin{equation*} \begin{array}{lcl} \left\{ \begin{array}{ll} \dot{x}_6 = y_6,\\ \dot{y}_6 = \mu_{1}+\mu_{2}y_6+x_6^2+x_6y_6+o(|x_6,y_6 |^2), \end{array} \right. \end{array} \end{equation*} |
where
\mu_{1} = -n_{00}n_{11}^4,\,\,\,\,\,\, \mu_2 = -n_{01}n_{11}. |
Using Maple software, we have
\left|\frac{\partial(\mu_1, \mu_2)}{\partial(\lambda_1, \lambda_2)}\right|_{\lambda_1 = \lambda_2 = 0} = - \frac{\big[(c+2) x_{*}^{2}+(d+s+1)x_{*}+s\big]^{4} (2 x_{*}^{2}+s x_{*}-s)^{5}}{s^{3} x_{*}^{8} (1+x_{*})^{2} T_1^{4}}\neq0. |
By the results in [28], system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two. The proof is completed.
The local expression of the bifurcation curves are given in [28] as follows:
(ⅰ) The saddle-node bifurcation curve
SN = \{(\lambda_{1}, \lambda_{2}): \mu_{1}(\lambda_{1},\lambda_{2}) = 0, \mu_{2}(\lambda_{1}, \lambda_{2})\neq0\}; |
(ⅱ) The Hopf bifurcation curve
H = \left\{(\lambda_{1}, \lambda_{2}): \mu_{1}(\lambda_{1},\lambda_{2}) < 0, \mu_{2}(\lambda_{1}, \lambda_{2}) = \sqrt{-\mu_{1}(\lambda_{1}, \lambda_{2})}\right\}; |
(ⅲ) The homoclinic curve
HL = \left\{(\lambda_{1}, \lambda_{2}): \mu_{1}(\lambda_{1}, \lambda_{2}) < 0, \mu_{2}(\lambda_{1}, \lambda_{2}) = \frac{5}{7}\sqrt{-\mu_{1}(\lambda_{1}, \lambda_{2})}\right\}. |
In what follows, we present the phase diagrams of system (5.4), as obtained by some numerical simulations. Choosing c = 2, d = 1, s = 4, q = 20, k = \frac{1}{6} and r = \frac{49}{3} , and from Theorem 4.5(1), E_{*}(1, 1) is a cusp of codimension two. Figure 6 shows that system (1.4) undergoes a Bogdanov-Takens bifurcation of codimension two.
In this paper, we consider a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Fear of predator and nonlinear harvesting are the main factors affecting the dynamic behavior of system (1.4). Via numerical simulations, we show the influences of the fear effect and nonlinear harvesting on the dynamic behavior of system (1.4).
First, let c = 0.1, d = 0.1, s = 0.04, q = 2 and r = 2 . When k = 0.06 , system (1.4) has no positive equilibrium (see Figure 7(a)). In a biological sense, the prey and predator will become extinct when the fear effect is large. When k = 0.053 , system (1.4) has two positive equilibria, where E_3 is a saddle and E_4 is an unstable node (see Figure 7(b)). Hence, in this case, the prey and predator are still extinct. When k = 0.05 , E_4 becomes a stable node, and an unstable limit cycle appears around E_4 (see Figure 7(c)). Then, for system (1.4), a bistable phenomenon occurs, in which the prey and predator tend to steady states (or extinction), depending on the initial values lying inside (or outside) of the unstable limit cycle. When k = 0.03 , E_4 is still a stable node and an unstable limit cycle disappears (see Figure 7(d)). Then, the prey and predator will survive or become extinct depending on the two stable manifolds of the saddle that act as a separatrix curve. When k = 0 , that is, without the fear effect, the dynamic behavior of system (1.4) is similar to that shown in Figure 7(d) (see Figure 7(e)). Figure 7 shows that the prey and predator may survive or become extinct when the fear effect is small. With the increase of the fear effect, the survival area of species decreases, until finally, the prey and predator will become extinct if the fear effect is strong enough. Hence, a strong fear effect is not conducive to the survival of the species.
Second, we consider the impact of nonlinear harvesting on system (1.4). Let c = 0.1, d = 0.1, s = 0.1, k = 0.1 and r = 3 . When q = 3.5 , system (1.4) has no positive equilibrium and the origin is globally asymptotically stable (see Figure 8(a)). When q = 3.3 , system (1.4) has two positive equilibria, where E_3 is a saddle and E_4 is an unstable node (see Figure 8(b)). In this case, the prey and predator are still extinct. When q = 3.287 , system (1.4) has an unstable limit cycle and there is a bistable phenomenon (see Figure 8(c)). That is, an unstable limit cycle acts as a separatrix curve, where the prey and predator will become extinct or survive. When q = 3.283 , there exists an unstable homoclinic loop in system (1.4) (see Figure 8(d)). When q = 3.26 , the unstable limit cycle and homoclinic loop disappear. Hence, the prey and predator will tend to steady states (or extinction) if the initial values lies to the right (or left) of the two stable manifolds of the saddle (see Figure 8(e)). When q = 0 , that is, without nonlinear harvesting, system (1.4) has only one positive equilibrium, which is globally asymptotically stable (see Figure 8(f)). This shows that overfishing can lead to the extinction of the predator and prey, so, maintaining proper harvesting can help the survival of the prey and predator.
By conducting numerical simulations, we were able to clearly observe that, when k\le0.05 and q\le3.281 , the prey and predator tend to coexist around the stable positive equilibrium E_4 . In other words, by effectively controlling the harvesting, we can ensure that the prey's fear of being caught remains within a smaller range, which benefits the survival of both populations. That is, weaker fear effects and less capture are beneficial to the survival of both predator and prey. We have conducted a theoretical analysis of system (1.4) and obtained some conclusions. However, when it comes to solving practical problems, there are many external factors. The actual application of the model may be difficult to achieve in the short term.
When r-d > q , the origin is a repeller, the only boundary equilibrium is a saddle point and the unique positive equilibrium may be stable or unstable. From Remark 4.1, the unique positive equilibrium is globally asymptotically stable if q < 1 . Then, the prey and predator will tend to a positive coexistent steady state if the birth rate of the prey is high and the catchability coefficient is small. When r-d = q\leq q_0 , from Remark 4.2, the origin is globally asymptotically stable, which implies that the prey and predator will become extinct. When r-d = q > q_0 or r-d < q , system (1.4) may have zero, one or two positive equilibria, and these equilibria may be stable or unstable. We show that the unique equilibrium E_{*} is a saddle-node or a cusp of codimension two (or three). Moveover, system (1.4) undergoes saddle-node bifurcation and Bogdanov-Takens bifurcation around E_{*} . Also, system (1.4) undergoes a degenerate Hopf bifurcation and multiple limit cycles may appear around \tilde{E_4} . In Figure 5, we show that system (1.4) has two limit cycles (the inner one is stable and the outer is unstable) around \tilde{E}_4 , which implies the bistable phenomenon. That is a large amount of fear and prey harvesting are detrimental to the survival of the prey and predator. Additionally, the prey and predator will reach a steady state if the intrinsic growth rate of the prey is high and the catchability coefficient for the prey is low. However, the prey and predator will become extinct if the intrinsic growth rate for the prey and the catchability coefficient for the prey are small.
In [9], the authors studied the stability of the equilibria and demonstrated that there exists a limit cycle in system (1.1). Considering the Holling type Ⅱ functional response, [10] showed that a unique equilibrium is a cusp of codimension two and a limit cycle appears. Unlike [9,10], we show that a unique equilibrium is a cusp of codimension three and confirm the occurrence of Bogdanov-Takens bifurcation. We have found that system (1.4) has two limit cycles (the inner one is stable and the outer is unstable), which exhibit the bistable phenomenon. Also, we have proven that the origin and equilibrium are globally asymptotically stable under some conditions. The strong fear effect and nonlinear harvesting are not conducive to the survival of the species. These indicate that the dynamic behavior of system (1.4) is more complex than that of the systems in [9,10].
The authors declare that they have not used artificial intelligence tools in the creation of this article.
This work was supported by the Natural Science Foundation of Fujian Province (2021J01613, 2021J011032) and the Scientific Research Foundation of Minjiang University (MJY22027).
The authors declare that there is no conflict of interest.
\begin{equation*} \begin{split} &P_0 = 2 (d+1)^2 (13 c^2-14 d-14) (c^2-d-1),\\ &P_1 = (d+1) (76 c^5-81 c^4 d+127 c^4-390 c^3 d+147 c^2 d^2-390 c^3-354 c^2 d+328 c d^2-60 d^3-501 c^2\\ &\qquad+656 c d+268 d^2+328 c+716 d+388),\\ &P_2 = 58 c^6-234 c^5 d+38 c^4 d^2+70 c^5-1230 c^4 d+915 c^3 d^2-79 c^2 d^3-852 c^4-1290 c^3 d+2439 c^2 d^2\\ &\qquad -612 c d^3+18 d^4-2205 c^3+2523 c^2 d+2100 c d^2-656 d^3+5 c^2+6036 c d+612 d^2+3324 c\\ &\qquad +3264 d+1978,\\ &P_3 = -177 c^6+100 c^5 d+115 c^4 d^2-1158 c^5+1979 c^4 d-63 c^3 d^2-74 c^2 d^3-2064 c^4+7346 c^3 d\\ &\qquad -2310 c^2 d^2+16 c d^3+26 d^4+1169 c^3+9954 c^2 d-5080 c d^2+36 d^3+8734 c^2+5536 c d-3024 d^2\\ \end{split} \end{equation*} |
\begin{equation*} \begin{split} &\qquad +10632 c+1116 d+4150,\\ &P_4 = 70 c^6+350 c^5 d-130 c^4 d^2+1427 c^5+883 c^4 d-1098 c^3 d^2+36 c^2 d^3+6845 c^4-2440 c^3 d\\ &\qquad -2088 c^2 d^2+344 c d^3-6 d^4+13906 c^3-12212 c^2 d-1692 c d^2+380 d^3+13912 c^2-16064 c d\\ &\qquad -696 d^2+6964 c-6780 d+1470,\\ &P_5 = 275 c^6-380 c^5 d-11 c^4 d^2+1227 c^5-3244 c^4 d+372 c^3 d^2+55 c^2 d^3+395 c^4-9320 c^3 d\\ &\qquad +2241 c^2 d^2+72 c d^3-6 d^4-7492 c^3-12879 c^2 d+3884 c d^2+44 d^3-18425 c^2-9696 c d\\ &\qquad +1992 d^2-17220 c-3276 d-5730,\\ &P_6 = -290 c^6-46 c^5 d+66 c^4 d^2-2724 c^5+672 c^4 d+567 c^3 d^2-31 c^2 d^3-9888 c^4+4870 c^3 d\\ &\qquad +1203 c^2 d^2-148 c d^3-18717 c^3+11647 c^2 d+1224 c d^2-124 d^3-21411 c^2+11340 c d\\ &\qquad +540 d^2-14416 c+3820 d-4268,\\ &P_7 = -43 c^6+204 c^5 d-23 c^4 d^2+312 c^5+1603 c^4 d-303 c^3 d^2+3204 c^4+4450 c^3 d-1116 c^2 d^2\\ &\qquad +9837 c^3+6564 c^2 d-1420 c d^2+13344 c^2+5192 c d-576 d^2+8084 c+1664 d+1728,\\ &P_8 = 162 c^6-70 c^5 d+1355 c^5-761 c^4 d+4549 c^4-2836 c^3 d+8456 c^3-4580 c^2 d+9068 c^2\\ &\qquad -3328 c d+5120 c-896 d+1152,\\ &P_9 = -(c+1) (55 c^5+530 c^4+1812 c^3+2752 c^2+1920 c+512),\\ &M_0 = 8 (c+2)^9 (3 c+d+3),\\ &M_1 = 4 (c+2)^8 (533 c^2+184 c d+15 d^2+1106 c+196 d+573),\\ &M_2 = 2 (c+2)^7 (9109 c^3+4629 c^2 d+751 c d^2+39 d^3+28867 c^2+10354 c d+891 d^2+30371 c+5737 d\\ &\qquad +10613),\\ &M_3 = (c+2)^6 (69015 c^4+45510 c^3 d+10460 c^2 d^2+970 c d^3+29 d^4+289652 c^3+154958 c^2 d+25896 c d^2\\ &\qquad +1334 d^3+454154 c^2+174818 c d+15852 d^2+315428 c+65354 d+81911),\\ &M_4 = 2 (c+2)^5 (73321 c^5+58214 c^4 d+16620 c^3 d^2+2014 c^2 d^3+87 c d^4+378112 c^4+263852 c^3 d\\ &\qquad +62784 c^2 d^2+5788 c d^3+156 d^4+776396 c^3+445684 c^2 d+78358 c d^2+4098 d^3+793634 c^2\\ &\qquad +332644 c d+32310 d^2+403923 c+92598 d+81894),\\ &M_5 = (c+2)^4 (189851 c^6+173592 c^5 d+57866 c^4 d^2+8320 c^3 d^3+435 c^2 d^4+1150784 c^5+976160 c^4 d\\ &\qquad +293568 c^3 d^2+36592 c^2 d^3+1560 c d^4+2891094 c^4+2177892 c^3 d+552334 c^2 d^2+52892 c d^3\\ &\qquad +1392 d^4+3853380 c^3+2410536 c^2 d+457164 c d^2+25160 d^3+2873687 c^2+1324004 c d\\ &\qquad +140568 d^2+1136796 c+288792 d+186328),\\ &M_6 = 2 (c+2)^3 (3 c^2+c d+5 c+2 d+2) (25505 c^5+17632 c^4 d+3965 c^3 d^2+290 c^2 d^3+134140 c^4\\ &\qquad +83604 c^3 d+16324 c^2 d^2+980 c d^3+280723 c^3+147364 c^2 d+22110 c d^2+824 d^3+292258 c^2\\ &\qquad +114464 c d+9852 d^2+151382 c+33072 d+31212)\\ &M_7 = (c+2)^2 (8325 c^4+3848 c^3 d+435 c^2 d^2+36898 c^3+14820 c^2 d+1380 c d^2+61025 c^2+18844 c d\\ &\qquad +1092 d^2+44632 c+7896 d+12180) (3 c^2+c d+5 c+2 d+2)^2,\\ &M_8 = 2 (c+2) (377 c^3+87 c^2 d+1355 c^2+258 c d+1618 c+192 d+640)(3 c^2+c d+5 c+2 d+2)^3,\\ &M_9 = (29 c^2+80 c+56) (3 c^2+c d+5 c+2 d+2)^4.\\ \end{split} \end{equation*} |
\begin{equation*} \begin{split} &\tilde{c}_{20} = \frac{n_1k^2+n_2 k+c^{2} s+c^{2}}{(\alpha+1)(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{c}_{11} = \frac{2\sqrt{D}\big[n_3k^2+(2 \alpha c s+2 \alpha c+4 c^{2}+2 c d+2 c s+4 c) k+c^{2}\big]}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{c}_{02} = - \frac{(\alpha s+\alpha+c+d+s+2) k^{2}D}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{c}_{30} = - \frac{n_4k^3+n_5k^2+n_6k+c^{3} (s+1) (2 \alpha+1)}{(\alpha+1)^{2}(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{c}_{21} = - \frac{\sqrt{D}\big[n_7k^3+n_8 k^2+3 c^{2} (\alpha s+\alpha+2 c+d+s+2) k+c^{3}\big]}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{c}_{12} = \frac{(s \alpha+\alpha+c+d+s+2) k^{2} (2 \alpha k s+2 \alpha k+4 c k+2 d k-k s+2 c+4 k)D}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3} },\\ &\tilde{c}_{03} = - \frac{(s \alpha+\alpha+c+d+s+2) k^{3} D^{\frac{3}{2}}}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{d}_{20} = - \frac{n_9k^2+n_{10}k+\alpha c^{2} s-\alpha c \,s^{2}+2 c^{2} s-c \,s^{2}+c^{2}}{\sqrt{D}\, (\alpha+1) (\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{d}_{11} = - \frac{[n_{11}k^2+n_{12} k+3 c^{2}-2 c s] s}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{d}_{02} = \frac{\sqrt{D}s[(2 \alpha s+2 \alpha+3 c+2 d+2 s+4) k^{2}+(\alpha s+\alpha+3 c+d+s+2) k+c]}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{2}},\\ &\tilde{d}_{30} = \frac{s\big[n_4k^3+n_5k^2+n_6k+c^{3} (s+1) (2 \alpha+1)\big]}{\sqrt{D}(\alpha+1)^{2} (\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{d}_{21} = \frac{s\big[n_{13}k^3+n_{14}k^2+3 c^{2} (\alpha s+\alpha+2 c+d+s+2)k+c^{3}\big]}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{d}_{12} = - \frac{s (s \alpha+\alpha+c+d+s+2) k^{2}(2 \alpha k s+2 \alpha k+4 c k+2 d k-k s+2 c+4 k)\sqrt{D}}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &\tilde{d}_{03} = \frac{s (s \alpha+\alpha+c+d+s+2) k^{3}D}{(\alpha k s+\alpha k+2 c k+d k+k s+c+2 k)^{3}},\\ &n_1 = (2 s^{2}+3 s+1) \alpha^{2}+(3 s^{2} c+d \,s^{2}+8 s c+4 d s+5 s^{2}+4 c+2 d+10 s+4) \alpha+4 c^{2} s+4 c d s+3 s^{2} c\\ &\quad\quad +d^{2} s+d \,s^{2}+4 c^{2}+4 d c+12 s c+d^{2}+6 d s+3 s^{2}+8 c+4 d+8 s+4,\\ &n_2 = (2 s^{2} c+4 s c+2 c) \alpha+4 c^{2} s+2 c d s+2 s^{2} c+4 c^{2}+2 d c+6 s c+4 c,\\ \end{split} \end{equation*} |
\begin{equation*} \begin{split} &n_3 = (s^{2}+2 s+1) \alpha^{2}+(4 c s+2 d s+s^{2}+4 c+2 d+5 s+4) \alpha+4 c^{2}+4 d c+3 c s+d^{2}+d s+8 c+4 d\\ &\quad\quad+2 s+4,\\ &n_4 = (4 s^{3}+10 s^{2}+8 s+2) \alpha^{4}+w_1 \alpha^{3}+w_2 \alpha^{2}+w_3 \alpha+w_4,\\ &n_5 = 2 c (s+1) (2 s^{2}+6 s+3)\alpha^{3}+w_5\alpha^{2}+w_6 \alpha+c \,s^{3}+(10 c^{2}+4 c d+11 c) s^{2}+(12 c^{3}+12 c^{2} d+3 c \,d^{2}\\ &\quad\quad +36 c^{2}+18 c d+24 c) s+12 c^{3}+12 c^{2} d+3 c \,d^{2}+24 c^{2}+12 c d+12 c,\\ &n_6 = 6 c^{2} (s+1)^{2}\alpha^{2}+3 c^{2} (s+1) (3 s+4 c+2 d+5)\alpha+3 c^{2} (s+1) (s+2 c+d+2),\\ &n_7 = (s+1)^{3} \alpha^{3}+(s+1)^{2} (6 c+3 d-s+6) \alpha^{2}+(s+1) (12 c^{2}+12 c d+3 d^{2}-2 d s-2 s^{2}+24 c+12 d\\ &\quad\quad -4 s+12) \alpha+8 c^{3}+(12 d+4 s+24) c^{2}+(6 d^{2}-3 s^{2}+24 d+24) c+(d+2) (d+s+2) (d-2 s+2),\\ &n_8 = 3 c (s+1)^{2}\alpha^{2}+2 c (s+1) (s+6 c+3 d+6)\alpha+12 c^{3}+(12 d+8 s+24) c^{2}+(d+s+2) (3 d-s+6) c,\\ &n_9 = s (s+1)^{2} \alpha^{3}+(s+1) (4 c s+2 d s+2 s^{2}+7 s+1) \alpha^{2}+(s^{3}+(9 c+4 d+12) s^{2}+(4 c^{2}+4 c d+d^{2}\\ &\quad\quad +20 c+10 d+18) s+4 c+2 d+4) \alpha+(5 c+2 d+5) s^{2}+2 (2 c+d+3) (2 c+d+2)s+(2 c+d+2)^{2},\\ &n_{10} = s (s+1) (-s+2 c) \alpha^{2}+(-2 s^{3}+(3 c-d-3) s^{2}+(4 c^{2}+2 c d+10 c) s+2 c) \alpha+(8 s+4) c^{2}+(4 d s\\ &\quad\quad +s^{2}+2 d+10 s+4) c-s^{2} (d+s+2),\\ &n_{11} = 3 (s+1)^{2}\alpha^{2}+2 (s+1) (s+6 c+3 d+6) \alpha+12 c^{2}+(12 d+6 s+24) c+(d+s+2) (3 d-s+6),\\ &n_{12} = 2 (s+1) (-s+3 c)\alpha+12 c^{2}+(6 d+12) c-2 s (d+s+2),\\ &n_{13} = (s+1)^{3} \alpha^{3}+(s+1)^{2} (6 c+3 d-s+6) \alpha^{2}+w_{7}\alpha+w_{8},\\ &n_{14} = 3 c (s+1)^{2} \alpha^{2}+2 c (s+1) (s+6 c+3 d+6) \alpha+12 c^{3}+(12 d+8 s+24) c^{2}+(d+s+2) (3 d-s+6)c,\\ &w_1 = (6 c+2 d+15)s^{3}+(30 c+14 d+48)s^{2}+(36 c+18 d+46) s+12 c+6 d+13,\\ &w_2 = (13 c+4 d+19) s^{3}+(20 c^{2}+18 c d+4 d^{2}+90 c+41 d+79) s^{2}+(48 c^{2}+48 c d+12 d^{2}+138 c\\ &\quad\quad +69 d+93) s+3 (4 c+2 d+5) (2 c+d+2),\\ &w_3 = (8 c+2 d+9) s^{3}+(28 c^{2}+24 c d+5 d^{2}+78 c+34 d+51) s^{2}+(16 c^{3}+24 c^{2} d+12 c \,d^{2}+2 d^{3}\\ &\quad\quad +96 c^{2}+96 c d+24 d^{2}+156 c+78 d+76) s+(4 c+2 d+7) (2 c+d+2)^{2},\\ &w_4 = (c+1) s^{3}+(8 c^{2}+6 c d+d^{2}+18 c+7 d+10) s^{2}+(8 c^{3}+12 c^{2} d+6 c \,d^{2}+d^{3}+36 c^{2}+36 c d\\ &\quad\quad +9 d^{2}+48 c+24 d+20) s+(2 c+d+2)^{3},\\ &w_5 = 9 c \,s^{3}+(22 c^{2}+10 c d+49 c) s^{2}+(48 c^{2}+24 c d+69 c) s+3 c (8 c+4 d+9),\\ &w_6 = 6 c \,s^{3}+(32 c^{2}+14 c d+44 c) s^{2}+(24 c^{3}+24 c^{2} d+6 c \,d^{2}+96 c^{2}+48 c d+78 c) s\\ &\quad\quad +6 c (2 c+d+3) (2 c+d+2),\\ &w_{7} = (s+1) (12 c^{2}+12 c d+3 d^{2}-2 d s-2 s^{2}+24 c+12 d-4 s+12),\\ &w_{8} = (-3 c-2 d-4) s^{2}+(4 c^{2}-d^{2}-4 d-4) s+(2 c+d+2)^{3}.\\ \end{split} \end{equation*} |
[1] |
J. R. Beddington, J. G. Cooke, Harvesting from a prey-predator complex, Ecol. Model., 14 (1982), 155–177. https://doi.org/10.1016/0304-3800(82)90016-3 doi: 10.1016/0304-3800(82)90016-3
![]() |
[2] |
D. M. Xiao, L. S. Jennings, Bifurcations of a ratio-dependent predator-prey system with constant rate harvesting, SIAM J. Appl. Math., 65 (2005), 737–753. https://doi.org/10.1137/s0036139903428719 doi: 10.1137/s0036139903428719
![]() |
[3] |
P. H. Leslie, J. C. Gower, The properties of a stochastic model for the predator-prey type of interaction between two species, Biometrika, 47 (1960), 219–234. https://doi.org/10.1093/biomet/47.3-4.219 doi: 10.1093/biomet/47.3-4.219
![]() |
[4] |
S. Rana, S. Bhattacharya, S. Samanta, Spatiotemporal dynamics of Leslie-Gower predator-prey model with Allee effect on both populations, Math. Comput. Simul., 200 (2022), 32–49. https://doi.org/10.1016/j.matcom.2022.04.011 doi: 10.1016/j.matcom.2022.04.011
![]() |
[5] |
M. X. He, Z. Li, Global dynamics of a Leslie-Gower predator-prey model with square root response function, Appl. Math. Lett., 140 (2023), 108561. https://doi.org/10.1016/j.aml.2022.108561 doi: 10.1016/j.aml.2022.108561
![]() |
[6] |
X. Q. Wang, Y. P. Tan, Y. L. Cai, W. M. Wang, Impact of the fear effect on the stability and bifurcation of a Leslie-Gower predator-prey model, Int. J. Bifurcation Chaos, 30 (2020), 2050210. https://doi.org/10.1142/S0218127420502107 doi: 10.1142/S0218127420502107
![]() |
[7] |
C. Arancibia-Ibarra, J. Flores, Dynamics of a Leslie-Gower predator-prey model with Holling type Ⅱ functional response, Allee effect and a generalist predator, Math. Comput. Simul., 188 (2021), 1–22. https://doi.org/10.1016/j.matcom.2021.03.035 doi: 10.1016/j.matcom.2021.03.035
![]() |
[8] |
J. Huang, Y. Gong, S. Ruan, Bifurcation analysis in a predator-prey model with constant-yield predator harvesting, Discrete Continuous Dyn. Syst. Ser. B, 18 (2013), 2101–2121. https://doi.org/10.3934/dcdsb.2013.18.2101 doi: 10.3934/dcdsb.2013.18.2101
![]() |
[9] |
R. P. Gupta, M. Banerjee, P. Chandra, Bifurcation analysis and control of Leslie-Gower predator-prey model with Michaelis-Menten type prey-harvesting, Differ. Equ. Dyn. Syst., 20 (2012), 339–366. https://doi.org/10.1007/s12591-012-0142-6 doi: 10.1007/s12591-012-0142-6
![]() |
[10] |
R. P. Gupta, P. Chandra, Bifurcation analysis of modified Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting, J. Math. Anal. Appl., 398 (2013), 278–295. https://doi.org/10.1016/j.jmaa.2012.08.057 doi: 10.1016/j.jmaa.2012.08.057
![]() |
[11] |
S. Kumar, H. Kharbanda, Chaotic behavior of predator-prey model with group defense and non-linear harvesting in prey, Chaos, Solitons Fractals, 119 (2019), 19–28. https://doi.org/10.1016/j.chaos.2018.12.011 doi: 10.1016/j.chaos.2018.12.011
![]() |
[12] |
T. Caraballo Garrido, R. Colucci, L. Guerrini, On a predator prey model with nonlinear harvesting and distributed delay, Commun. Pure Appl. Anal., 17 (2018), 2703–2727. https://doi.org/10.3934/cpaa.2018128 doi: 10.3934/cpaa.2018128
![]() |
[13] |
D. Hu, H. Cao, Stability and bifurcation analysis in a predator-prey system with Michaelis-Menten type predator harvesting, Nonlinear Anal. Real World Appl., 33 (2017), 58–82. https://doi.org/10.1016/j.nonrwa.2016.05.010 doi: 10.1016/j.nonrwa.2016.05.010
![]() |
[14] |
C. Zhu, L. Kong, Bifurcations analysis of Leslie-Gower predator-prey models with nonlinear predator-harvesting, Discrete Continuous Dyn. Syst. Ser. S, 10 (2017), 1187–1206. https://doi.org/10.3934/dcdss.2017065 doi: 10.3934/dcdss.2017065
![]() |
[15] |
R. Cristiano, M. M. Henao, D. J. Pagano, Global stability of a Lotka-Volterra piecewise-smooth system with harvesting actions and two predators competing for one prey, J. Math. Anal. Appl., 522 (2023), 126998. https://doi.org/10.1016/j.jmaa.2023.126998 doi: 10.1016/j.jmaa.2023.126998
![]() |
[16] |
R. Sivasamy, K. Sathiyanathan, K. Balachandran, Dynamics of a modified Leslie-Gower model with gestation effect and nonlinear harvesting, J. Appl. Anal. Comput., 9 (2019), 747–764. https://doi.org/10.11948/2156-907x.20180165 doi: 10.11948/2156-907x.20180165
![]() |
[17] |
X. Yan, C. Zhang, Global stability of a delayed diffusive predator-prey model with prey harvesting of Michaelis-Menten type, Appl. Math. Lett., 114 (2021), 106904. https://doi.org/10.1016/j.aml.2020.106904 doi: 10.1016/j.aml.2020.106904
![]() |
[18] |
X. Wang, L. Zanette, X. Zou, Modelling the fear effect in predator-prey interactions, J. Math. Biol., 73 (2016), 1179–1204. https://doi.org/10.1007/s00285-016-0989-1 doi: 10.1007/s00285-016-0989-1
![]() |
[19] |
M. M. Chen, Y. Takeuchi, J. F. Zhang, Dynamic complexity of a modified Leslie-Gower predator-prey system with fear effect, Commun. Nonlinear Sci. Numer. Simul., 119 (2023), 107109. https://doi.org/10.1016/j.cnsns.2023.107109 doi: 10.1016/j.cnsns.2023.107109
![]() |
[20] |
X. B. Zhang, H. L. Hu, Q. An, Dynamics analysis of a diffusive predator-prey model with spatial memory and nonlocal fear effect, J. Math. Anal. Appl., 525 (2023), 127123. https://doi.org/10.1016/j.jmaa.2023.127123 doi: 10.1016/j.jmaa.2023.127123
![]() |
[21] |
C. M. Zhang, S. L. Liu, J. H. Huang, W. M. Wang, Stability and Hopf bifurcation in an eco-epidemiological system with the cost of anti-predator behaviors, Math. Biosci. Eng., 20 (2023), 8146–8161. https://doi.org/10.3934/mbe.2023354 doi: 10.3934/mbe.2023354
![]() |
[22] |
Y. J. Li, M. X. He, Z. Li, Dynamics of a ratio-dependent Leslie-Gower predator-prey model with Allee effect and fear effect, Math. Comput. Simul., 201 (2022), 417–439. https://doi.org/10.1016/j.matcom.2022.05.017 doi: 10.1016/j.matcom.2022.05.017
![]() |
[23] |
J. X. Zhao, Y. F. Shao, Bifurcations of a prey-predator system with fear, refuge and additional food, Math. Biosci. Eng., 20 (2023), 3700–3720. https://doi.org/10.3934/mbe.2023173 doi: 10.3934/mbe.2023173
![]() |
[24] |
M. He, Z. Li, Stability of a fear effect predator-prey model with mutual interference or group defense, J. Biol. Dyn., 16 (2022), 480–498. https://doi.org/10.1080/17513758.2022.2091800 doi: 10.1080/17513758.2022.2091800
![]() |
[25] |
D. Pal, D. Kesh, D. Mukherjee, Qualitative study of cross-diffusion and pattern formation in Leslie-Gower predator-prey model with fear and Allee effects, Chaos, Solitons Fractals, 167 (2023), 113033. https://doi.org/10.1016/j.chaos.2022.113033 doi: 10.1016/j.chaos.2022.113033
![]() |
[26] | Z. Zhang, T. Ding, W. Huang, Z. Dong, Qualitative Theory of Differential Equations, Translations of Mathematical Monographs, American Mathematical Society, 1992. |
[27] |
J. C. Huang, Y. J. Gong, J. Chen, Multiple bifurcations in a predator-prey system of Holling and Leslie type with constant-yield prey harvesting, Int. J. Bifurcation Chaos, 23 (2013), 1350164. https://doi.org/10.1142/s0218127413501642 doi: 10.1142/s0218127413501642
![]() |
[28] | L. Perko, Differential Equations and Dynamical Systems, Springer, New York, 1996. https://doi.org/10.1007/978-1-4684-0392-3 |
[29] |
L. Yang, Recent advances on determining the number of real roots of parametric polynomials, J. Symb. Comput., 28 (1999), 225–242. https://doi.org/10.1006/jsco.1998.0274 doi: 10.1006/jsco.1998.0274
![]() |
[30] |
Y. Dai, Y. Zhao, B. Sang, Four limit cycles in a predator-prey system of Leslie type with generalized Holling type Ⅲ functional response, Nonlinear Anal. Real World Appl., 50 (2019), 218–239. https://doi.org/10.1016/j.nonrwa.2019.04.003 doi: 10.1016/j.nonrwa.2019.04.003
![]() |
[31] |
M. Lu, J. Huang, S. Ruan, P. Yu, Bifurcation analysis of an SIRS epidemic model with a generalized nonmonotone and saturated incidence rate, J. Differ. Equations, 267 (2019), 1859–1898. https://doi.org/10.1016/j.jde.2019.03.005 doi: 10.1016/j.jde.2019.03.005
![]() |
1. | Hongqiuxue Wu, Zhong Li, Mengxin He, Bifurcation Analysis of a Holling–Tanner Model with Generalist Predator and Constant-Yield Harvesting, 2024, 34, 0218-1274, 10.1142/S0218127424500767 | |
2. | Manoj Kumar Singh, Arushi Sharma, Luis M. Sánchez-Ruiz, Dynamical Complexity of Modified Leslie–Gower Predator–Prey Model Incorporating Double Allee Effect and Fear Effect, 2024, 16, 2073-8994, 1552, 10.3390/sym16111552 | |
3. | Xia Qi, Dongpo Hu, Ming Zhao, Huijing Sun, Ming Liu, The impact of Michaelis–Menten predator harvesting and weak Allee effect on prey in a Leslie–Gower predator–prey system: a comprehensive dynamics study, 2025, 0924-090X, 10.1007/s11071-025-11208-x | |
4. | Manisha Yadav, Pradeep Malik, Md. Jasim Uddin, Nehad Ali Shah, Modeling and managing chaos: a fear-affected discretized prey–predator system with Allee Effect, 2025, 140, 2190-5444, 10.1140/epjp/s13360-025-06354-5 |