Processing math: 76%
Research article

Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey

  • Received: 15 February 2023 Revised: 04 April 2023 Accepted: 10 April 2023 Published: 23 April 2023
  • We incorporate the strong Allee effect and fear effect in prey into a Leslie-Gower model. The origin is an attractor, which implies that the ecological system collapses at low densities. Qualitative analysis reveals that both effects are crucial in determining the dynamical behaviors of the model. There can be different types of bifurcations such as saddle-node bifurcation, non-degenerate Hopf bifurcation with a simple limit cycle, degenerate Hopf bifurcation with multiple limit cycles, Bogdanov-Takens bifurcation, and homoclinic bifurcation.

    Citation: Zhenliang Zhu, Yuming Chen, Zhong Li, Fengde Chen. Dynamic behaviors of a Leslie-Gower model with strong Allee effect and fear effect in prey[J]. Mathematical Biosciences and Engineering, 2023, 20(6): 10977-10999. doi: 10.3934/mbe.2023486

    Related Papers:

    [1] 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
    [2] 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
    [3] Hongqiuxue Wu, Zhong Li, Mengxin He . Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Mathematical Biosciences and Engineering, 2023, 20(10): 18592-18629. doi: 10.3934/mbe.2023825
    [4] 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
    [5] 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
    [6] Kawkab Al Amri, Qamar J. A Khan, David Greenhalgh . Combined impact of fear and Allee effect in predator-prey interaction models on their growth. Mathematical Biosciences and Engineering, 2024, 21(10): 7211-7252. doi: 10.3934/mbe.2024319
    [7] 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
    [8] 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
    [9] Zhenzhen Shi, Huidong Cheng, Yu Liu, Yanhui Wang . Optimization of an integrated feedback control for a pest management predator-prey model. Mathematical Biosciences and Engineering, 2019, 16(6): 7963-7981. doi: 10.3934/mbe.2019401
    [10] Kunlun Huang, Xintian Jia, Cuiping Li . Analysis of modified Holling-Tanner model with strong Allee effect. Mathematical Biosciences and Engineering, 2023, 20(8): 15524-15543. doi: 10.3934/mbe.2023693
  • We incorporate the strong Allee effect and fear effect in prey into a Leslie-Gower model. The origin is an attractor, which implies that the ecological system collapses at low densities. Qualitative analysis reveals that both effects are crucial in determining the dynamical behaviors of the model. There can be different types of bifurcations such as saddle-node bifurcation, non-degenerate Hopf bifurcation with a simple limit cycle, degenerate Hopf bifurcation with multiple limit cycles, Bogdanov-Takens bifurcation, and homoclinic bifurcation.



    Predation is a ubiquitous interaction among species and thus it has been studied extensively. Most predator-prey models are extensions and modifications of the classical Lotka-Volterra model by incorporating many other factors. The obtained theoretical results can guide us to make appropriate policies on ecological protection and ecological sustainable development [1]. One important modification is the Leslie-Gower model described by

    dxdt=rx(1xK)P(x)y,dydt=ys(1ynx), (1.1)

    where x(t) and y(t) stand for the densities of the prey and predator at time t, respectively. Here both species have the logistic growths with respective intrinsic growth rates r and s and with the carrying capacities K for the prey and nx for the predator (which depends on the density of the prey). n is a measure of the quality of the prey and nx is called the Leslie-Gower term [2]. P(x) is the functional response of predators to prey. Model (1.1) has been extensively studied and promoted (see, for example, [3,4]). Recently, the influence of fear effect and Allee effect on the dynamical behaviors of predator-prey models has attracted many researchers' attention.

    Fear effect of predator is a form of indirect predation impacting the prey population [5,6,7]. To some extent, the physiological states of prey are disturbed by the fear effect. As a result, in the long run, the fear effect may lead to a loss in prey population. Meanwhile, the prey suffered from scare usually forage less and then their birth rate decreases which may be not good for the survival of prey's population. There is an example that some birds make anti-predator defenses against the sound of the predator and once perceive danger, they flee from their nests. As far as prey is concerned, frequent exposure to fear and anti-predation behavior will affect their birth rate. Therefore, it is necessary to take into account of fear factor in the predator-prey interaction [8,9,10]. In [11], Wang et al. first modeled the fear effect by introducing a factor f(k,y), where k represents the intensity of anti-predator behaviors of the prey caused by fear. See the paper for more detail on the properties of f(k,y). A recent experiment designed by Elliott et al. [12] demonstrated that Allee effect [13,14,15] can be induced by the fear effect. Moreover, the strong Allee effects [16,17] should be paid more attention, because too large fear intensity will drive the species to extinction.

    Allee effect mainly signifies that individual fitness is directly proportionate to population density. Due to the variety of species, the causes of Allee effect are also different. Such as sperm limitation, reproduction facilitation, cooperative breeding, anti-predator behavior and predator satiation. During the last decades, lots have been done for predator-prey models with Allee effects. For example, Merdan [18] investigated the stabilityof a Lotka-Volterra type predator-prey system involving Allee effect; Liu and Dai [19] considered an impulsive predator-prey model with double Allee effects; Wu et al. [20] tried to understand the Allee effect in a commensal symbiosis mode; studied Guan and Chen [21] studied an amensalism model with Allee effect on the second species; and Naik et al. [22] explored the effect of strong Allee effect in a discrete predator-prey model.

    On the one hand, based on the experimental observations of Elliott et al. [12], Sasmal [15] proposed and investigated the following predator-prey model with prey subject to Allee effect and fear effect,

    dxdt=rx(1xK)(xθ)11+fyaxy,dydt=αaxymy,

    where θ and f are Allee constant and fear factor, respectively. They showed that the cost of fear has no influence on the stability of equilibria. On the other hand, González-Olivares et al. [23] studied the following Leslie-Gower model with Allee effect in prey,

    dxdt=rx(1xK)(xm)qxy,dydt=sy(1ynx),

    where m represents the Allee effect threshold. Without Allee effect, there is a unique equilibrium, which is globally asymptotically stable. However, with the presence of strong Allee effect, the number of equilibria varies with respect to m. It is shown that the system exhibits local bifurcations including Hopf bifurcation and Bogdanov-Takens bifurcation.

    Inspired by these two works [15,23], it is natural and interesting to investigate the combined influence of both Allee effect and fear effect on dynamics of Leslie-Gower predator-prey models. In other words, whether the qualitative structure and bifurcation phenomena will become more complex? The thoretical results will help humans manage ecosystems. Precisely, the model to be studied in this paper is

    dxdt=rx(1xK)(xm)11+fyqxy,dydt=ys(1ynx), (1.2)

    where r, K, f, q, n, s, and m are all positive constants and 0<m<K. For simplicity of analysis, we introduce new variables

    ˉx=xK,ˉy=1nKy,ˉt=rKt,

    and denote

    ˉm=mK,ˉq=nqr,ˉs=srk,ˉf=fnK.

    Then system (1.2) becomes

    dxdt=x(1x)(xm)11+fyqxy,dydt=ys(1yx), (1.3)

    (note that we have dropped the bars here), where m, f, q, and s are positive constants with 0<m<1.

    Though not defined at x=0, the dynamical behavior of (1.3) in the absence of prey is of interest. For this end, we take the time scaling dt=(1+fy)xdτ to get a new system (still label τ as t),

    dxdt=x2(1x)(xm)qx2y(1+fy), (1.4a)
    dydt=ys(xy)(1+fy). (1.4b)

    This new system with the time scale reaches stable equilibria much faster than the original one. As a result, in the following we only focus on (1.4). Due to the biological background, the initial conditions (x(0),y(0))R2+. Obviously, for such an initial condition, (1.4) has a unique global solution which stays in R2+.

    In the following, we first present the results on boundedness of solutions and the attraction of the origin. Followed is the existence and local stability of equilibria. The results indicate possible occurrence of bifurcations. We thus explore the detail on saddle-node bifurcation of codimension 1, Hopf bifurcation of codimension 1, and Bogdanov-Takens bifurcation of codimension 2. A brief conclusion ends the paper.

    We first show the ultimate boundedness of solutions of (1.4).

    Proposition 1. Let (x(t),y(t)) be a solution of (1.4). Then

    lim suptx(t)1andlim supty(t)1.

    Proof. We show lim suptx(t)1 by distinguishing two cases. Firstly, we assume x(t)1 for t0. Then it follows from (1.4a) that x is a decreasing function. Denote limtx(t) by l. It is easy to see that l=1. Next, we assume that x(t0)<1 for some t00. We claim that, for tt0, x(t)1. Otherwise, there exists a t>t0 such that x(t)>1. Let x(ˆt)=max{x(t)|t0tt}. As dx(t)dt<0, we have ˆt(t0,t) and hence dx(ˆt)dt=0. This is impossible as x(ˆt)>1, which gives dx(ˆt)dt<0. The claim is proved. To sum up, lim suptx(t)1 is proved. Now, for ε0>0, there exists ˘t0 such that x(t)1+ε0 for t˘t. It follows that, for t˘t,

    dy(t)dtys(1+ε0y)(1+fy).

    Then arguing similarly as for lim suptx(t)1, we can get lim supty(t)1+ε0. As ε0 is arbitrary, we immediately have lim supty(t)1.

    We start with finding the equilibria to (1.4), which are determined by solutions to

    {x2(1x)(xm)qx2y(1+fy)=0,ys(xy)(1+fy)=0.

    Clearly, (1.4) always admits two positive boundary equilibria P1(m,0) and P2(1,0). Moreover, for a positive equilibrium (x,x) (or coexistence equilibrium), x satisfies

    (1+fq)x2+(mq+1)xm=0. (3.1)

    If mq+10 then (3.1) has no positive roots and hence (1.4) has no positive equilibria. Now assume that mq+1>0, which implies that q<m+1<2. Denote the discriminant of (3.1) by Δ and it can be expressed as a quadratic function of m,

    Δ(m)=m22(2fq+q+1)m+(q1)2.

    Note that Δ(1)=q(q44f)<0 and Δ(q1)=2(1q)(2+fq). It follows that if q1 then Δ(m)<0, which again implies that (1.4) has no positive equilibrium. Finally, if 0<q<1, then Δ(m)>0 only when 0<m<m1 while Δ(m)=0 only when m=m1, where

    m1=2fq+q+12f2q2+fq2+fq+q.

    In the former, (3.1) has two positive roots, x3=mq+1Δ(m)2(1+fq) and x4=mq+1+Δ(m)2(1+fq), which are in (0,1) while in the latter, ˉx=1fq+qfq+1(0,1) is the unique positive root. The above discussion is summarized as follows.

    Theorem 2. (ⅰ) There are always three boundary equilibria (0,0), P1(m,0), and P2(1,0) for system (1.4).

    (ⅱ) Besides the three boundary equilibria, for (1.4) to have positive equilibria it is necessary that 0<q<1. In this case,

    (a) if m1<m<1, there is no positive equilibrium;

    (b) if m=m1, there is a unique positive equilibrium ˉP(ˉx,ˉx), where ˉx=1fq+qfq+1;

    (c) if 0<m<m1, there are distinct positive equilibria P3(x3,x3) and P4(x4,x4), where x3=mq+1Δ(m)2(1+fq) and x4=mq+1+Δ(m)2(1+fq).

    Remarks 3. Considering f=0, we have q+12q:=m0. Obviously, m1<m0, which means that influenced by fear effect, Allee threshold value decreases. In other words, prey are more likely to die out at low density.

    Next, we consider the attractivity of the origin.

    Theorem 4. The origin is a non-hyperbolic attractor for (1.4).

    Proof. Note that the technique of linearization is inapplicable as J(0,0) is the zero matrix for (1.4). The approach here is the blow-up method. With the horizontal blow-up

    (x,y)=(x,ux)anddτ=xdt

    along the invariant line x=0, we rewrite system (1.4) as

    dx(τ)dτ=x[m(m+1)x+qux+x2+fqu2x2],du(τ)dτ=u[m+ssu(m+1)x+(fs+q)ux+x2fsu2x+fqu2x2].

    On the nonnegative u-axis, the above system admits two equilibria C1(0,0) and C2(0,m+ss). The corresponding Jacobian matrices are

    JC1=(m00m+s)andJC2=(m0(m+1)(m+s)sms),

    respectively. Obviously, C1 is a saddle while C2 is an attractor. After blow-down, the origin is an attractor.

    Theorem 5. P1 is a hyperbolic repeller node while P2 is a saddle.

    Proof. At P1 and P2, the Jacobian matrices are respectively

    JP1=(m(1m)00s)andJP2=(m100s).

    respectively. Clearly, since m<1, JP1 has two positive eigenvalues while JP2 has one positive and one negative eigenvalues. Then the results follow immediately.

    The result below indicates that the stability of the positive equilibrium ˉP depends on s.

    Theorem 6. Suppose that 0<q<1 and m=m1, which guarantee the existence of ˉP. Then we have the following statements.

    (i) ˉP is a saddle-node with an attracting parabolic sector when s>s, where s=qˉx(1+2fˉx)fˉx+1;

    (ii) ˉP is a saddle-node with a repelling parabolic sector when s<s;

    (iii) ˉP is a degenerate equilibrium when s=s. Furthermore, if either (0<ˉx<23 and ff1) or 23ˉx<1, then ˉP is a cusp of codimension two; if 0<ˉx<23 and f=f1, then ˉP is a cusp of codimension at least three, where

    f1=2ˉx25ˉx+1+ˉx2ˉx24ˉx+12ˉx2(2ˉx).

    Fig. 1 shows the phase portraits.

    Figure 1.  Phase portraits near ˉP.

    Proof. Note that the Jacobian matrix of (1.4) at ˉP is

    JˉP=(ˉx2(2ˉxm11)qˉx2(1+2fˉx)ˉxs(1+fˉx)ˉxs(1+fˉx)).

    We easily see that

    det(ˉP)=sˉx31+fˉx1+fq(ˉx1+fq+qfq+1)=0

    and hence JˉP has zero as an eigenvalue.

    First of all, we use the transformation

    X=xˉx,Y=yˉy

    to translate ˉP into (0,0) and system (1.4) into

    dXdt=a10X+a01Y+a20X2+a11XY+a02Y2+a30X3+a21X2Y+a12XY2+a03Y3+O(|X,Y|4),dYdt=b10X+b01Y+b20X2+b11XY+b02Y2+b30X3+b21X2Y+b12XY2+b03Y3+O(|X,Y|4), (3.2)

    where

    a10=qˉx2(1+2fˉx),a01=qˉx2(1+2fˉx),a20=5ˉx2+(2m1+2)ˉx,a11=2qˉx(1+2fˉx),a02=qfˉx2,a30=m1+14ˉx,a21=q(1+2fˉx),a12=2fqˉx,a03=0,b10=s(1+fˉx),b01=s(1+fˉx),b20=0,b11=s(1+2fˉx),b02=s(1+2fˉx),b30=b21=0,b12=fs,b03=fs.

    Now, assume that ss. Then the linear transformation

    X=u+v,Y=u+b01a01v

    is nonsingular. This, combined with the time scale dτ=(a10+b01)dt, transforms system (3.2) into (we still denote τ as t)

    dudt=a20u2+a11uv+a02v2+O(|u,v|3),dvdt=v+b20u2+b11uv+b02v2+O(|u,v|3), (3.3)

    where

    a20=sˉx(1+fˉx)m1(a10+b01)2,a11={sˉx(4ˉx4f2+2(f2sfm1f+4)fˉx3+2(2fm1+3fs2m1fq2)fˉx2+(6fm1+5fsq)ˉx+2m1+s)}(1+2fˉx)(a10+b01)2,a02=s(1+fˉx)20f2qˉx5q(1+2fˉx)2(a10+b01)2+s(1+fˉx)(16f3s10f2m110f2+18f)qˉx4q(1+2fˉx)2(a10+b01)2+s(1+fˉx)(3f3s2+4f2m1q+28f2qs9fm1q9fq+4q)ˉx3q(1+2fˉx)2(a10+b01)2+s(1+fˉx)(6f2s2+4fm1q+16fqs2m1q2q)ˉx2q(1+2fˉx)2(a10+b01)2+s(1+fˉx)(4fs2+m1q+3qs)ˉxq(1+2fˉx)2(a10+b01)2+s(1+fˉx)s2q(1+2fˉx)2(a10+b01)2,b20=s(1+fˉx)ˉxm1(a10+b01)2,b11=12fqˉx5(a10+b01)2+(10f2qs6fm1q+8fq26fq+12q)ˉx4(a10+b01)2+(f2s2+2fm1q+6fqsm1qq)ˉx3(a10+b01)2+(3fs2+2m1q+3qs)ˉx2(a10+b01)2+s2ˉx(a10+b01)2,b02=20f2q2ˉx6q(1+fˉx)(a10+b01)2+(8f3s10f2m110f2+18f)q2ˉx5q(1+fˉx)(a10+b01)2+(5f3qs2+4f2m1q2+16f2q2s9fm1q29fq2+4q2)ˉx4q(1+fˉx)(a10+b01)2+(2f3s3+10f2qs2+4fm1q2+10fq2s2m1q22q2)ˉx3q(1+fˉx)(a10+b01)2+(5f2s3+6fqs2+m1q2+2q2s)ˉx2q(1+fˉx)(a10+b01)2+(qs24fs3)ˉxq(1+fˉx)(a10+b01)2s3q(1+fˉx)(a10+b01)2.

    Since a20>0, with the center manifold method [24], the equation near the center manifold can be approximated by

    dudt=s(1+fˉx)ˉxm1(a10+b01)2u2+O(|u|3). (3.4)

    Making use of [25,Theorem 7.1], the degenerate equilibrium ˉP is a saddle-node where the parabolic sector is located at the right half plane. From the time scaling dτ=(a10+b01)dt, we see that if s>s, the parabolic sector is attracting while if s<s the parabolic sector is repelling.

    Next, we assume s=s. It follows from

    ˉx=m1q+12(1+fq)andm1=2fq+q+12f2q2+fq2+fq+q

    that

    m1=ˉx2(f+1)1ˉx2f+2ˉxfandq=q(ˉx1)21ˉx2f+2ˉxf.

    Obviously, m1(0,1), q>0, and s>0 still hold under the assumptions that f>0 and ˉx(0,1). Then (X,Y)=(u,a01a10u+1a10v) transforms (3.2) into

    dudt=v+c20u2+c11uv+c02v2+O(|u,v|3),dvdt=d20u2+d11uv+d02v2+O(|u,v|3), (3.5)

    where

    c20=ˉx2(f+1)ˉx2f2ˉxf1,c11=2(3ˉxf+1)ˉx(2ˉxf+1),c02=f(ˉx2f2ˉxf1)(2ˉxf+1)2ˉx2(ˉx1)2,d20=ˉx4(f+1)(ˉx1)2(2ˉxf+1)(ˉx2f2ˉxf1)2,d11=ˉx(ˉx1)2(2ˉx2f2+4ˉxf+1)(ˉxf+1)(ˉx2f2ˉxf1)2,d02=3ˉx2f2+3ˉxf+1ˉx(2ˉxf+1)(ˉxf+1).

    Thus the C change of coordinates in a small neighborhood of (0,0),

    z1=uc11+d122u2c02uv,z2=v+c20u2d02uv,

    produces the normal form of system (3.5),

    dz1dt=z2+O(|z1,z2|3),dz2dt=f20z21+f11z1z2+O(|z1,z2|3), (3.6)

    where

    f20=ˉx4(f+1)(1+ˉx)2(2ˉxf+1)(ˉx2f2ˉxf1)2,f11=ˉx((2ˉx44ˉx3)f2+(4ˉx310ˉx2+2ˉx)f+ˉx24ˉx+1)(ˉxf+1)(1ˉx2f+2ˉxf).

    Obviously, f20<0 and the sign of f11 relies on the sign of

    φ(f)(2ˉx44ˉx3)f2+(4ˉx310ˉx2+2ˉx)f+ˉx24ˉx+1.

    Note that 2ˉx44ˉx3<0 since 0<ˉx<1. Denote

    Δ1(4ˉx310ˉx2+2ˉx)24(2ˉx44ˉx3)(ˉx24ˉx+1)=4(2ˉx24ˉx+1)ˉx2(ˉx1)2.

    Then

    Δ1{<0 if 222<ˉx<1>0 if 0<ˉx<222=0 if ˉx=222

    Thus φ(f)<0 if 222<ˉx<1. If ˉx=222, then φ(f)=2232(f+1)2<0 since f>0. Now consider the case where 0<ˉx<222. Notice that

    4ˉx310ˉx2+2ˉx{>0 if 0<ˉx<5174,<0 if 5174<ˉx<222,=0 if ˉx=5174,

    and

    ˉx24ˉx+1{>0if0<¯x<23,=0if¯x=23,<0if23<ˉx<1.

    It follows that φ(f)<0 if 23<ˉx<222 and φ(f)=f[(90523)f+(38223)]<0 when ˉx=23 since f>0. If 0<ˉx<23, then

    φ(f){>0if0<f<f1,=0iff=f1,<0iff>f1,

    where

    f1=2ˉx25ˉx+1+(1ˉx)2ˉx24ˉx+12ˉx2(2ˉx).

    Then the results follow and this completes the proof.

    At last, we consider the positive equilibria P3 and P4.

    Theorem 7. Suppose 0<q<1 and 0<m<m1. Then P3 is always a saddle point while P4 is a sink if s>x4(2x4+m+1)fx4+1, is a source if s<x4(2x4+m+1)fx4+1, and a center or fine focus if s=x4(2x4+m+1)fx4+1.

    Proof. After a simple calculation, we can get

    det(JP3)=sx33Δ(m)1+fx31+fq<0anddet(JP4)=sx34Δ(m)1+fx41+fq>0,

    where JP3 and JP4 are the Jacobian matrices of (1.4) at P3 and P4, respectively. It follows immediately that P3 is always a saddle. For the stability of P4, we need the trace of JP4, which is

    Tr(JP4)=x4[2x24+(fsm1)x4+s]{<0ifs>x4(2x4+m+1)fx4+1,>0ifs<x4(2x4+m+1)fx4+1,0ifs=x4(2x4+m+1)fx4+1.

    Then the results follow easily.

    The goal of this section is to analyse the saddle-node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation that may occur in system (1.4).

    From Theorems 6 and 7, it is not difficult to obtain the saddle-node surface,

    SN={(m,q,f,s):m=m1,ss,f>0,s>0,q>0,0<m<1}.

    The phase of the saddle-node bifurcation corresponding to the ecology system has obvious critical Allee value m1. When m>m1, the prey will face the risk of extinction; when m<m1, the dynamic behavior of system (1.4) becomes complex since the saddle and anti-saddle come out in the phase, which means that the density of prey must be large enough for survival.

    From Theorem 7, the local stability of P4 will change as the parameter s varies. Besides, there exists an sH=x4(2x4+m+1)fx4+1 which satisfies Tr(JP4)|sH=0 and then P4 becomes a non-hyperbolic equilibrium which loses its stability. Meanwhile, it is easy to verify the transversality condition,

    ddsTr(JP4)|sH=x4(fx4+1)<0.

    Therefore, there exists a Hopf bifurcation in a small neighbourhood of P4. This is summarized as follows.

    Theorem 8. Let s be the Hopf bifurcation parameter. When s crosses the threshold value sH=x4(2x4+m+1)fx4+1, system (1.4) undergoes a non-degenerate Hopf bifurcation around P4.

    There is at least one limit cycle which appears around P4 due to the occurrence of Hopf bifurcation. The stability of the limit cycle is very significant for ecology systems and we can use the sign of the Lyapunov number to determine it. In order to simplify the computation, we translate P4(x4,x4) to (1,1) with the following state variable scaling and time rescaling:

    ˉx=xx4,ˉy=yx4,τ=x24t,ˉK=1x4,ˉm=mx4,ˉq=q,ˉf=fx4,ˉs=sx4. (4.1)

    After dropping the bars and relabelling τ as t, we transform system (1.4) into

    dxdt=x2(1xK)(xm)qx2y(1+fy),dydt=ys(xy)(1+fy). (4.2)

    Because system (4.2) has an equilibrium ¯P4(1,1) (i.e., P4(x4,x4) of (1.4)), we have q=(K1)(1m)K(1+f)>0. Then we can get another positive equilibrium ¯P3(ˉx3,ˉx3), where ˉx3 and ˉx4 are positive distinct roots of the equation

    (1m)(K1)+(1+f)K(1+f)x2+(m+K)(1+f)(1m)(K1)K(1+f)xm=0.

    By Vieta theorem [26], ˉx3ˉx4=Km(1+f)Kfm+Kf+fm+1<1. From the scalings (4.1), we see that the parameters in system (4.2) satisfy

    K>1,m<1,f>0,s>0,Km(1+f)Kfm+Kf+fm+1<1. (4.3)

    With q being replaced by (K1)(1m)K(1+f), system (4.2) becomes

    dxdt=x2(1xK)(xm)(K1)(1m)K(1+f)x2y(1+fy),dydt=ys(xy)(1+fy), (4.4)

    where the parameters K, m, f, and s satisfy (4.3). Therefore, we study the stability of the weak focus ¯P4(1,1) for system (4.4) to acquire the stability of the weak focus P4(x4,y4) in system (1.4).

    The Jacobian matrix of system (4.4) at ¯E4(1,1) is

    J¯E4=(2+m+KK(K1)(1+m)(1+2f)K(1+f)s(1+f)s(1+f)).

    It follows that

    det(J¯P4)=(2KfmKf+Kmfm1)sK

    and

    Tr(J¯P4)=Kfs+KsKm+2K.

    With the help of (4.3), we see that det(J¯P4)>0. Thus the necessary conditions on parameters for Hopf bifurcation occurring around ¯E4 are

    m+K2K(1+f)>0,K>1,m<1,f>0,s>0,Km(1+f)Kfm+Kf+fm+1<1. (4.5)

    Note that at s=shfm+K2K(1+f), J¯P4 has a pair of conjugate pure imaginary eigenvalues. Moreover,

    ddsTr(J¯P4)|shf=(f+1)<0,

    that is, the transversality condition holds.

    Now translating ¯P4(1,1) to the origin (0,0) by the change (ˆx,ˆy)=(x1,y1) when s=shf, we rewrite (4.4) as

    {dˆxdt=^a10ˆx+^a01ˆy+^a20ˆx2+^a11ˆxˆy+^a02ˆy2+^a30ˆx3+^a21ˆx2ˆy+^a12ˆxˆy2+^a03ˆy3+O(|ˆx,ˆy|4),dˆydt=^b10ˆx+^b01ˆy+^b20ˆx2+^b11ˆxˆy+^b02ˆu2+^b30ˆx3+^b21ˆx2ˆy+^b12ˆxˆy2+^b03ˆy3+O(|ˆx,ˆy|4), (4.6)

    where

    ^a10=2+m+KK,^a01=(K1)(1+m)(1+2f)K(1+f),^a20=2K+2m5K,^a11=2(K1)(1+m)(1+2f)K(1+f),^a02=2(K1)(1+m)(1+2f)K(1+f),^a30=4+K+mK,^a21=(K1)(1+m)(1+2f)K(1+f),^a12=2(K1)(1+m)fK(1+f),^a03=0,^b10=2+m+KK,^b01=2+m+KK,^b20=0,^b11=(2+m+K)(1+2f)K(1+f),^b02=(2+m+K)(1+2f)K(1+f),^b30=0,^b21=0,^b12=(2+m+K)fK(1+f),^b03=(2+m+K)fK(1+f).

    We now make another transformation ˆu=ˆx, ˆv=1D(^a10ˆx+^a01ˆy), and dτ=Ddt to change system (4.6) into

    {dˆudt=ˆv+f(ˆu,ˆv),dˆvdt=ˆu+g(ˆu,ˆv),

    where

    D=(2mK)(2KfmKf+Kmfm1)K2(1+f),f(ˆu,ˆv)=^c20ˆu2+^c11ˆuˆv+^c02ˆv2+^c30ˆu3+^c21ˆu2ˆv+^c12ˆuˆv2+^c03ˆv3+O(|ˆu,ˆv|4),g(ˆu,ˆv)=^d20ˆu2+^d11ˆuˆv+^d02ˆv2+^d30ˆu3+^d21ˆu2ˆv+^d12ˆuˆv2+^d03ˆv3+O(|ˆu,ˆv|4).

    Here we have omitted the expressions of ^cij and ^dij.

    Applying the formal series method in [25] and calculating the first Lyapunov number with the help of MAPLE, we get

    l1=18p1f2+p2f+p3Dp4,

    where

    p1=(2K1)2m3+(4K332K2+33K8)m2(4K1)(K28K+6)m+K38K2+6K,p2=(4K23K)m3+(4K330K2+24K1)m2+(3K3+24K212K6)mK26K+6,p3=(K2K)m3+K(K1)(K7)m2(K1)(K26K2)m2K+2,p4=(1+2f)(1+m)(K1)K((2KmKm)f+Km1).

    Under condition (4.5), we have p4>0. Thus the sign of l1 depends on that of

    l11=p1f2+p2f+p3.

    Theorem 9. If l11>0 then system (4.4) exhibits a subcritical Hopf bifurcation around ¯P4 at s=shf (see Fig. 2(a)) while if l11<0 then it exhibits a supercritical Hopf bifurcation around ¯P4 at s=shf (see Fig. 2(b)).

    Figure 2.  Hopf bifurcation around ¯E4 of (4.4).

    When l11=0, we should consider the second Lyapunov number l2 to determine the stability of the order-two weak focus. The expression of l2 is too complicated, so we just provide a numerical example with (K,f,m,s)=(2,3628121430,25,428657711). For such a set of parameter values, we have l1=0 and l2=0.4025314467>0, which means that system (4.4) undergoes a multiple Hopf bifurcation of codimension two and there are two limit cycles (the inner one is stable and the outer one is unstable) bifurcated from ¯P4(1,1) when the bifurcation parameters (f,s)=(3628121430,428657711) is disturbed (see Fig. 2(c)).

    From Theorem 6, system (1.4) admits that ˉP is a cusp of codimension 2 when m=m1, s=s, q=q, and φ(f)0. In the following, we show that Bogdanov-Takens bifurcation around ˉP can occur.

    Theorem 10. Under the condition that ˉP is a cusp of codimension 2, we choose the parameters m and s as bifurcation parameters. As the parameters (m,s) vary around (m1,s), system (1.4) experiences Bogdanov-Takens bifurcation in a small neighborhood of\ ˉP. Furthermore,

    (i) Given φ(f)>0, system (1.4) exhibits a supercritical Bogdanov -Takens bifurcation, which accompanies with the appearance of a stable limit cycle originated from the supercritical Hopf bifurcation followed by a stable homoclinic loop;

    (ii) Given φ(f)<0, system (1.4) exhibits a subcritical Bogdanov -Takens bifurcation, which accompanies with the appearance of an unstable limit cycle originated from the subcritical Hopf bifurcation followed by an unstable homoclinic loop.

    Proof. We give the unfolding system of system (1.4),

    dxdt=x2(1x)(xm1λ1)qx2y(1+fy),dydt=y(s+λ2)(xy)(1+fy), (4.7)

    where λ1 and λ2 are disturbed in a small neighbourhood of the origin. For convenience, we should find the versal unfolding of system (4.7), which is helpful to study the change of the topological structure. We make some transformations as follows to achieve it.

    Firstly, use the transformation u1=xˉx and u2=yˉx to translate the positive equilibrium ˉP to (0,0) and transform system (4.7) into

    du1dt=˜a00+˜a10u1+˜a01u2+˜a20u21+˜a11u1z2+˜a02u22+O(|u1,u2|3),du2dt=˜b10u1+˜b01u2+˜b20u21+˜b11u1u2+˜b02u22+O(|u1,u2|3), (4.8)

    where

    ˜a00=ˉx3λ1ˉx2λ1,˜a10=ˉx(2ˉx4f+(3fλ14f+1)ˉx3+(8fλ1+2f2)ˉx2+(4fλ1+3λ1+1)ˉx2λ1)2ˉxf+1ˉx2f,˜a01=(ˉx22ˉx+1)ˉx2(2ˉxf+1)ˉx2f2ˉxf1,˜a20=5ˉx4f+(3fλ110f+2)ˉx3+(7fλ1+4f5)ˉx2+(2fλ1+3λ1+2)ˉxλ12ˉxf+1ˉx2f,˜a11=2ˉx(ˉx22ˉx+1)(2ˉxf+1)ˉx2f2ˉxf1,˜a02=(ˉx22ˉx+1)ˉx2fˉx2f2ˉxf1,˜b10=(2ˉx4f+(f2λ24f+1)ˉx3+(2f2λ2fλ2+2f2)ˉx2+(3fλ2+1)ˉx+λ2)ˉx2ˉxf+1ˉx2f,˜b01=˜b10,˜b20=0,˜b11=(2ˉx4f+(f2λ24f+1)ˉx3+(2f2λ2fλ2+2f2)ˉx2(3fλ2+1)ˉx+λ2)(2ˉxf+1)(ˉxf+1)(ˉx2f2ˉxf1),˜b02=˜b11.

    Next, letting u3=u2, u4=du2dt, we transform system (4.8) into

    du3dt=u4,du4dt=˜c00+˜c10u3+˜c01u4+˜c20u23+˜c11u3u4+˜c02u24+O(|u3,u4|3), (4.9)

    where

    ˜c00=˜a00˜b10,˜c10=˜b11˜a00˜a10˜b01+˜a01˜b10,˜c01=˜a10+˜b01,˜c20=˜a10˜b10˜b02˜a01˜b10˜b11˜a20˜b201+˜a11˜b01˜b10˜a02˜b210˜b10,˜c11=2˜a20~b01˜a11˜b102˜b02˜b10+˜b11˜b01˜b10,˜c02=˜a20+˜b11˜b10.

    In order to remove the ˜c02u24 term in system (4.9), we make a new time variable τ by dt=(1˜c02u3)dτ and let u5=u3, u6=u4(1˜c02u3). We still denote τ as t to get from (4.9) that

    \begin{equation} \begin{array}{rcl} \frac{du_5}{dt}& = & u_6, \\ \frac{du_6}{dt}& = & \tilde{d}_{00}+\tilde{d}_{10}u_5+\tilde{d}_{01}u_6 +\tilde{d}_{20}u_5^2 +\tilde{d}_{11}u_5u_6+O(|u_5,u_6|^3), \end{array} \end{equation} (4.10)

    where

    \begin{eqnarray*} &&\tilde{d}_{00} = \tilde{c}_{00}, \qquad \tilde{d}_{10} = \tilde{c}_{10}-2\tilde{c}_{02}\tilde{c}_{00}, \quad \tilde{d}_{01} = \tilde{c}_{01}, \\ && \tilde{d}_{20} = \tilde{c}_{20}-2\tilde{c}_{10}\tilde{c}_{02} +\tilde{c}_{00}\tilde{c}_{02}^2,\quad \tilde{d}_{11} = \tilde{c}_{11}-\tilde{c}_{01}\tilde{c}_{02}. \end{eqnarray*}

    Notice that \tilde{d}_{20}\to -\frac{(\bar{x}-1)^2(2\bar{x}f+1)\bar{x}^4(f+1)} {(\bar{x}^2f-2\bar{x}f-1)^2} < 0 when (\lambda_1, \lambda_2)\to (0, 0) . Then we employ the transformation

    u_7 = u_5,\qquad u_8 = \frac{u_6}{\sqrt{-\tilde{d}_{20}}},\qquad \tau = \sqrt{-\tilde{d}_{20}}t

    to change \tilde{d}_{20} into -1 . Meanwhile with this transform and relabelling \tau as t , system (4.10) becomes

    \begin{equation} \begin{array}{rcl} \frac{du_7}{dt} & = & u_8, \\ \frac{du_8}{dt}& = & \tilde{e}_{00}+\tilde{e}_{10}u_7+\tilde{e}_{01}u_8-u_7^2+\tilde{e}_{11}u_7u_8+O(|u_7,u_8|^3), \end{array} \end{equation} (4.11)

    where

    \tilde{e}_{00} = -\frac{\tilde{d}_{00}}{\tilde{d}_{20}},\quad \tilde{e}_{10} = -\frac{\tilde{d}_{10}}{\tilde{d}_{20}}, \quad \tilde{e}_{01} = \frac{\tilde{d}_{01}}{\sqrt{-\tilde{d}_{20}}}, \quad \tilde{e}_{11} = \frac{\tilde{d}_{11}}{\sqrt{-\tilde{d}_{20}}}.

    Now eliminate u_7 by the transformation u_9 = u_7-\frac{\tilde{e}_{10}}{2} , u_{10} = u_8 to change system (4.11) into

    \begin{equation} \begin{array}{rcl} \frac{du_9}{dt}& = & u_{10}, \\ \frac{du_{10}}{dt}& = & \tilde{f}_{00}+\tilde{f}_{01}u_{10}-u_9^2 +\tilde{f}_{11}u_9u_{10}+O(|u_9,u_{10}|^3), \end{array} \end{equation} (4.12)

    where

    \tilde{f}_{00} = \tilde{e}_{00}+\frac{1}{4}\tilde{e}_{10}^2, \qquad \tilde{f}_{01} = \tilde{e}_{01}+\frac{1}{2}\tilde{e}_{10}\tilde{e}_{11}, \qquad \tilde{f}_{11} = \tilde{e}_{11}.

    Lastly, to obtain the versal unfolding of (4.7), we need to make \tilde{f_{11}} be 1. Since \tilde{f}_{11}\to \frac{\bar{x}\varphi(f)}{(\bar{x}f+1)(1-\bar{x}^2f+2\bar{x}f)} \neq 0 as (\lambda_1, \lambda_2) \to (0, 0) , we can perform the transformation

    \begin{equation} x = -\tilde{f_{11}}^2u_9, \qquad y = \tilde{f}_{11}^3u_{10}, \qquad \tau = -\frac{t}{\tilde{f}_{11}}. \end{equation} (4.13)

    Relabelling \tau as t , we can rewrite (4.12) as

    \begin{equation} \begin{array}{rcl} \frac{dx}{dt}& = &y, \\ \frac{dy}{dt}& = &\tilde{g}_{00}+\tilde{g}_{01}y+x^2+xy+O(|x,y|^3), \end{array} \end{equation} (4.14)

    where

    \tilde{g}_{00} = -\tilde{f_{00}}\tilde{f_{11}}^4, \qquad \tilde{g}_{01} = -\tilde{f_{01}}\tilde{f_{11}}.

    With the assistance of MAPLE,

    \begin{equation} \frac{\partial(\tilde{g}_{00},\tilde{g}_{01})} {\partial(\lambda_1,\lambda_2)}\Big\vert_{\lambda_1 = \lambda_2 = 0} = \frac{\varphi^5(f)(\bar{x}^2f-2\bar{x}f-1)^2} {(f+1)^4(\bar{x}f+1)^4(\bar{x}-1)^5(2\bar{x}f+1)^3\bar{x}^6} \neq 0, \end{equation} (4.15)

    which implies that \tilde{g}_{00} and \tilde{g}_{10} are independent parameters. Therefore, system (4.14) is the versal unfolding of system (4.7) such that both have the same topological structure. As \tilde{g}_{00} and \tilde{g}_{10} vary, the topological structure of system (4.14) will change and system (1.4) will undergo Bogdanov-Takens bifurcation. Noticing that there exists time scale t = -\tilde{f_{11}}\tau in the transformation (4.13), we have the following conclusions:

    (i) if \varphi (f) > 0 , then -\tilde{f}_{11} < 0 . Thus system (1.4) undergoes a supercritical Bogdanov-Takens bifurcation of codimension 2, which consists of saddle-node bifurcation of codimension 1, supercritical Hopf bifurcation of codimension 1, and homoclinic bifurcation of codimension 1.

    (ii) if \varphi (f) < 0 , then -\tilde{f}_{11} > 0 . Thus system (1.4) undergoes a subcritical Bogdanov-Takens bifurcation of codimension 2, which consists of saddle-node bifurcation of codimension 1, subcritical Hopf bifurcation of codimension 1, and homoclinic bifurcation of codimension 1.

    Before drawing the phase portraits, we need to acquire the expression of the bifurcation curves within Bogdanov-Takens bifurcation in the (\lambda_1, \lambda_2) -plane.

    By the conclusion in Bogdanov [27] and Takens [28] (see also Chow et al. [29] and Perko [24]), we can obtain the following local representations of the bifurcation curves.

    (i) The saddle-node bifurcation curve SN = \{(\tilde{g}_{00}, \tilde{g}_{01}):\tilde{g}_{00} = 0, \tilde{g}_{01}\neq 0 \} ;

    (ii) the Hopf bifurcation curve H = \{(\tilde{g}_{00}, \tilde{g}_{01}):\tilde{g}_{01} = \sqrt{-\tilde{g}_{00}}, \tilde{g}_{00} < 0 \} ;

    (iii) the homoclinic bifurcation curve HL = \{(\tilde{g}_{00}, \tilde{g}_{01}):\tilde{g}_{01} = \frac{5}{7}\sqrt{-\tilde{g}_{00}}, \tilde{g}_{00} < 0 \} .

    In the case of repelling Bogdanov-Takens bifurcation ( \varphi (f) > 0 ), we will use \lambda_1 and \lambda_2 to describe the bifurcation curves SN , H , and HL . By (4.15) and the Implicit Function Theorem, we get the expressions of \lambda_1 and \lambda_2 in terms of \tilde{g_{00}} and \tilde{g_{10}} ,

    \begin{equation} \lambda_1 = a_1\tilde{g}_{00}+a_2\tilde{g}_{10}+O(\mid\tilde{g}_{00}, \tilde{g}_{10}\mid), \quad \lambda_2 = b_1\tilde{g}_{00}+b_2\tilde{g}_{10}+O(\mid\tilde{g}_{00}, \tilde{g}_{10}\mid), \end{equation} (4.16)

    where

    \begin{eqnarray*} a_1& = &\frac{\begin{array}{c}\bar{x}^6(\bar{x}^3f^2+\bar{x}^3f+\bar{x}^2f +\bar{x}^2)(\bar{x}f+1)^3(f+1)^2(2\bar{x}^3f -4\bar{x}^2f \\ +\bar{x}^2+2\bar{x}f-2\bar{x}+1)^2 \end{array}} { \begin{array}{c}(\bar{x}^2f-2\bar{x}f-1)(\bar{x}-1)(2\bar{x}^5f^2 -4\bar{x}^4f^2+4\bar{x}^4f \\ -10\bar{x}^3f +\bar{x}^3+2\bar{x}^2f-4\bar{x}^2+\bar{x})^4 \end{array}}, \\ a_2 & = &0, \\ b_1& = & \frac{\bar{b_1}\bar{x}^3(f+1)^2 (f\bar{x}+1)(\bar{x}-1)^2(2f\bar{x}+1)} {\varphi^4(f)(1-f\bar{x}^2+2f\bar{x})}, \\ b_2& = &\frac{\bar{x}^2(\bar{x}-1)^2(f+1)(2\bar{x}f+1)} {(1-\bar{x}^2f+2\bar{x}f)\varphi(f)} > 0, \\ \bar{b_1}& = & 16\bar{x}^8f^4-(68f^4-50f^3)\bar{x}^7 +(76f^4-246f^3+48f^2)\bar{x}^6 \\ && -(12f^4-330f^3+284f^2-17f)\bar{x}^5 \\ &&-(8f^4+114f^3-436f^2+131f-2)\bar{x}^4 \\ &&-(4f^3+188f^2-231f+21) \bar{x}^3+(12f^2-113f+43)\bar{x}^2 \\ &&+(12f-23)\bar{x}+3. \end{eqnarray*}

    For the saddle-node bifurcation curve SN , let \Lambda_1 \triangleq \tilde{g_{00}}(\lambda_1, \lambda_2) = 0 . Due to

    \frac{\partial\Lambda_1}{\partial\lambda_1}\Bigg\vert_{\lambda = 0} = \frac{\varphi^4(f)(\bar{x}^2f-2\bar{x}f-1)} {\bar{x}^4(\bar{x}-1)^3(2\bar{x}f+1)^2(f+1)^3(\bar{x}f+1)^4}\neq0,

    turning to the Implicit Function Theorem, there exists a unique function \lambda_1(\lambda_2) , which satisfies \lambda_1(0) = 0 and \Lambda_1(\lambda_1(\lambda_2), \lambda_2) = 0 . It follows that

    \lambda_1(\lambda_2) = 0.

    Furthermore, on the curve \Lambda_1 = 0 , from (4.16), we know \lambda_2 = b_2\tilde{g_{10}}+O(\mid\tilde{g_{10}}\mid^2) , which implies that \lambda_2 > 0 if \tilde{g_{10}} > 0 and \lambda_2 < 0 if \tilde{g_{10}} < 0 . Therefore, we acquire

    SN^+ = \{(\lambda_1,\lambda_2)\vert\lambda_1 = 0,\lambda_2 < 0\} \qquad {\rm{and}} \qquad SN^- = \{(\lambda_1,\lambda_2)\vert \lambda_1 = 0,\lambda_2 > 0\}.

    For the Hopf bifurcation curve H , let \Lambda_2 \triangleq \tilde{g_{00}}(\lambda_1, \lambda_2) +\tilde{g_{10}}^2(\lambda_1, \lambda_2) = 0 . Due to

    \frac{\partial\Lambda_2}{\partial\lambda_1} \Big\vert_{\lambda = 0} = \frac{\varphi^4(f) (\bar{x}^2f-2\bar{x}f-1)}{\bar{x}^4(\bar{x}-1)^3 (2\bar{x}f+1)^2(f+1)^3(\bar{x}f+1)^4}\neq0,

    there exists a unique function \lambda_1(\lambda_2) , which satisfies \lambda_1(0) = 0 and \Lambda_2(\lambda_1(\lambda_2), \lambda_2) = 0 . Then

    \lambda_1(\lambda_2) = \frac{(f\bar{x}^2-2f\bar{x}-1) (f+1)(f\bar{x}+1)^4}{\varphi^2(f)(1-\bar{x})} \lambda_2^2+O(\mid\lambda_2\mid^2).

    Furthermore, on the curve \Lambda_2 = 0 , from (4.16), we know \lambda_2 = b_2\tilde{g_{10}}+O(\mid\tilde{g_{10}}\mid^2) , which implies that \lambda_2 < 0 if \tilde{g_{10}} < 0 . Therefore, we get

    H = \left\{\left(\lambda_1, \lambda_2\right) \mid \lambda_1 = \frac{\left(f \bar{x}^2-2 f \bar{x}-1\right)(f+1)(f \bar{x}+1)^4}{\varphi^2(f)(1-\bar{x})} \lambda_2^2+O\left(\left|\lambda_2\right|^2\right)\right\}.

    For the Homoclinic bifurcation curve HL , let \Lambda_3 \triangleq \frac{25}{49}\tilde{g_{00}}(\lambda_1, \lambda_2) +\tilde{g_{10}}^2(\lambda_1, \lambda_2) = 0 . Due to

    \frac{\partial\Lambda_3}{\partial\lambda_1} \Bigg\vert_{\lambda = 0} = \frac{25}{49}\frac{\varphi^4(f) (\bar{x}^2f-2\bar{x}f-1)}{\bar{x}^4(\bar{x}-1)^3 (2\bar{x}f+1)^2(f+1)^3(\bar{x}f+1)^4}\neq0,

    there exists a unique function \lambda_1(\lambda_2) , which satisfies \lambda_1(0) = 0 and \Lambda_3(\lambda_1(\lambda_2), \lambda_2) = 0 . Then

    \lambda_1(\lambda_2) = \frac{49}{25} \frac{(f\bar{x}^2-2f\bar{x}-1)(f+1)(f\bar{x}+1)^4} {\varphi^2(f)(1-\bar{x})}\lambda_2^2+O(\mid\lambda_2\mid^2).

    Moreover, on the curve \Lambda_3 = 0 , from (4.16), we know \lambda_2 = b_2\tilde{g_{10}}+O(\mid\tilde{g_{10}}\mid^2) , which implies that \lambda_2 < 0 if \tilde{g_{10}} < 0 . Therefore,

    HL = \left\{(\lambda_1,\lambda_2)\Bigg\vert\lambda_1 = \frac{49}{25}\frac{(f\bar{x}^2-2f\bar{x}-1)(f+1) (f\bar{x}+1)^4}{\varphi^2(f)(1-\bar{x})}\lambda_2^2 +O(\mid\lambda_2\mid^2)\right\}.

    Analogously, the Bogdanov-Takens bifurcation curve in terms of \lambda_1 and \lambda_2 can be described.

    To end this section, we present two sets of numerical simulations to verify the theoretical results with the help of MATLAB.

    First, let \bar{x} = 0.5 and f = 3 . Then m_1 = \frac{4}{13} , q_* = \frac{1}{13} , s_* = \frac{4}{65} , and \varphi(f) = -\frac{57}{8} < 0 . Therefore, system (1.4) undergoes a repelling Bogdanov-Takens bifurcation. With the disturbance of (\lambda_1, \lambda_2) in a small neighborhood of (0, 0) , the repelling Bogdanov-Takens bifurcation diagram and its local phase portraits are shown in Fig 3.

    Figure 3.  Repelling Bogdanov-Taken bifurcation of (4.7).

    Next, let \bar{x} = 0.1 and f = 8 . Then m_1 = \frac{1}{28} , q_* = \frac{9}{28} , s_* = \frac{13}{280} , and \varphi(f) = 1.1988 > 0 . Therefore, system (1.4) undergoes an attracting Bogdanov-Takens bifurcation. With the disturbance of (\lambda_1, \lambda_2) in a small neighborhood of (0, 0) , the attracting Bogdanov-Takens bifurcation diagram and its local phase portraits are illustrated by Fig. 4.

    Figure 4.  Attracting Bogdanov-Takens bifurcation of (4.7).

    In the present study, the dynamical behavior of a Leslie-Gower model with linear (Holling Type I) functional response and with both Allee effect and fear effect in prey is investigated. Although the model is not defined when the density of the prey is zero, the origin is an attractor with the help of the blow-up method and the technique of time scaling. Correspondingly, prey will die out at low densities. Moreover, solutions are bounded. Based on the work of Korobeinikov [30] that system (1.4) admits a unique globally asymptotically stable positive equilibrium when m = 0 and f = 0 , qualitative analysis reveals that Allee effect results in complex dynamics. The number and stability of equilibria in system (1.4) depend on the Allee constant. The model has two boundary equilibria, one an unstable node and the other a saddle which means that the prey population will not fully reach the environmental capacity. There is a cusp coexistence equilibrium of codimension 2 or 3 if Theorem 6 holds. Furthermore, with the help of the formal series method in [25], we calculated the order of the weak focus to determine the stability of limit cycles bifurcated from Hopf bifurcation, which means that there is the possibility of periodic distributions of prey and predator populations. Choosing the Allee constant and intrinsic growth rate of the predator as bifurcation parameters and by means of a series of near-identity transformations, it was demonstrated that system (1.4) undergoes Bogdanov-Takens bifurcation of codimension 2 through calculating the versal unfolding.

    Compared to the work of González-Olivares et al. [23], we focus on the cost of the anti-predation response in the prey reproduction. Our investigation indicates that the cost of fear still does not chang the stability of equilibria. This is similar to the result that the stability of positive equilibrium will not be influenced and the system still undergoes analogous bifurcation.

    Different from the work of Sasmal [15], we considered a Leslie-Gower model where the carrying capacity of the predator relies on the abundance of prey. The model has abundant bifurcation phenomena. Through rigorous mathematical analysis, as the Allee constant varies, the model experiences saddle-node bifurcation and there are at most two equilibria in the phase diagram. Meanwhile, the model can undergo more complex bifurcation, Bogdanov-Takens bifurcation, than the traditional predator-prey model. Further, the model undergoes subcritical or supercritical Hopf bifurcation, which means that there is a stable or unstable limit cycle around a positive equilibrium. The model can even exhibit multiple Hopf bifurcations by a set of numerical simulations when the first Lyapunov number is zero.

    This work was supported partially by the Natural Science Foundation of Fujian Province (Nos. 2020J01499, 2021J01613) and the National Natural Science Foundation of China (No. 12071418).

    The authors declare there is no conflict of interest.



    [1] M. Yavuz, N. Sene, Stability analysis and numerical computation of the fractional predator–prey model with the harvesting rate, Fractal Fract., 4 (2020), 35. https://doi.org/10.3390/fractalfract4030035 doi: 10.3390/fractalfract4030035
    [2] 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
    [3] L. Chen, F. Chen, Global stability of a Leslie-Gower predator-prey model with feedback controls, Appl. Math. Lett., 22 (2009), 1330–1334. https://doi.org/10.1016/j.aml.2009.03.005 doi: 10.1016/j.aml.2009.03.005
    [4] S. Pal, N. Pal, S. Samanta, J. Chattopadhyay, Fear effect in prey and hunting cooperation among predators in a Leslie-Gower model, Math. Biosci. Eng., 16 (2019), 5146–5179. https://doi.org/10.3934/mbe.2019258 doi: 10.3934/mbe.2019258
    [5] J. Chen, X. He, F. Chen, The influence of fear effect to a discrete-time predator-prey system with predator has other food resource, Mathematics, 9 (2021), 1–20. https://doi.org/10.3390/math9080865 doi: 10.3390/math9080865
    [6] P. Cong, M. Fan, X. Zou, Dynamics of a three-species food chain model with fear effect, Commun. Nonlinear Sci. Numer. Simul., 99 (2021), 105809. https://doi.org/10.1016/j.cnsns.2021.105809 doi: 10.1016/j.cnsns.2021.105809
    [7] Y. Shi, J. Wu, Q. Cao, Analysis on a diffusive multiple Allee effects predator-prey model induced by fear factors, Nonlinear Anal. Real World Appl., 59 (2021), 103249. https://doi.org/10.1016/j.nonrwa.2020.103249 doi: 10.1016/j.nonrwa.2020.103249
    [8] Y. Huang, Z. Zhu, Z. Li, Modeling the Allee effect and fear effect in predator-prey system incorporating a prey refuge, Adv. Difference Equ., 2020 (2020), 321. https://doi.org/10.1186/s13662-020-02727-5 doi: 10.1186/s13662-020-02727-5
    [9] X. Wang, X. Tan, Y. Cai, W. Wang, Impact of the fear effect on the stability and bifurcation of a Leslie-Gower predator-prey model, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 30 (2020), 2050210. https://doi.org/10.1142/S0218127420502107 doi: 10.1142/S0218127420502107
    [10] Z. Zhu, R. Wu, L. Lai, X. Yu, The influence of fear effect to the Lotka-Volterra predator-prey system with predator has other food resource, Adv. Difference Equ. 2020 (2020), 237. https://doi.org/10.1186/s13662-020-02612-1 doi: 10.1186/s13662-020-02612-1
    [11] X. Wang, L. Zanette, X. Zou, Modelling the fear effect in predator-prey interactions, J. Math. Biol., 73 (2016), 1179–1204. 10.1007/s00285-016-0989-1 doi: 10.1007/s00285-016-0989-1
    [12] K. H. Elliott, G. S. Betini, D. R. Norris, Fear creates an Allee effect: experimental evidence from seasonal populations, Proc. Biol. Sci., 284 (2017), 20170878. https://doi.org/10.1098/rspb.2017.0878 doi: 10.1098/rspb.2017.0878
    [13] B. Dennis, Allee effects: population growth, critical density, and the chance of extinction, Natur. Resource Model., 3 (1989), 481–538. https://doi.org/10.1111/j.1939-7445.1989.tb00119.x doi: 10.1111/j.1939-7445.1989.tb00119.x
    [14] F. Courchamp, T. Clutton-Brock, B. Grenfell, Inverse density dependence and the Allee effect, Trends Ecol. Evol., 14 (1999), 405–410. https://doi.org/10.1016/S0169-5347(99)01683-3 doi: 10.1016/S0169-5347(99)01683-3
    [15] S. K. Sasmal, Population dynamics with multiple Allee effects induced by fear factors-a mathematical study on prey-predator interactions, Appl. Math. Model., 64 (2018), 1–14. https://doi.org/10.1016/j.apm.2018.07.021 doi: 10.1016/j.apm.2018.07.021
    [16] M. Verma, A. K. Misra, Modeling the effect of prey refuge on a ratio-dependent predator-prey system with the Allee effect, Bull. Math. Biol., 80 (2018), 626–656. 10.1007/s11538-018-0394-6 doi: 10.1007/s11538-018-0394-6
    [17] Z. Zhu, M. He, Z. Li, F. Chen, Stability and bifurcation in a logistic model with Allee effect and feedback control, Internat. J. Bifur. Chaos Appl. Sci. Engrg., 30 (2020), 2050231. https://doi.org/10.1142/S0218127420502314 doi: 10.1142/S0218127420502314
    [18] H. Merdan, Stability analysis of a Lotka-Volterra type predator-prey system involving Allee effect, ANZIAM J., 52 (2010), 139–145. https://doi.org/10.1017/S1446181111000630 doi: 10.1017/S1446181111000630
    [19] X. Liu, B. Dai, Dynamics of a predator-prey model with double Allee effects and impulse, Nonlinear Dynam., 88 (2017), 685–701. https://doi.org/10.1007/s11071-016-3270-7 doi: 10.1007/s11071-016-3270-7
    [20] R. Wu, L. Li, Q. Lin, A Holling type commensal symbiosis model involving Allee effect, Commun. Math. Biol. Neurosci., 2018 (2018), 6. https://doi.org/10.28919/cmbn/3679 doi: 10.28919/cmbn/3679
    [21] X. Guan, F. Chen, Dynamical analysis of a two species amensalism model with Beddington-DeAngelis functional response and Allee effect on the second species, Nonlinear Anal. Real World Appl., 48 (2019), 71–93. https://doi.org/10.1016/j.nonrwa.2019.01.002 doi: 10.1016/j.nonrwa.2019.01.002
    [22] P. A. Naik, Z. Eskandari, M. Yavuz, J. Zu, Complex dynamics of a discrete-time Bazykin–Berezovskaya prey-predator model with a strong Allee effect, J. Comput. Appl. Mathe., 413 (2022), 114401. https://doi.org/10.1016/j.cam.2022.114401 doi: 10.1016/j.cam.2022.114401
    [23] E. González-Olivares, J. Mena-Lorca, A. Rojas-Palma, J. Flores, Dynamical complexities in the Leslie-Gower predator-prey model as consequences of the Allee effect on prey, Appl. Math. Model., 35 (2011), 366–381. https://doi.org/10.1016/j.apm.2010.07.001 doi: 10.1016/j.apm.2010.07.001
    [24] L. Perko, Differential Equations and Dynamical Systems, 3rd ed., Springer-Verlag, New York, 2001.
    [25] Z. F. Zhang, T. R. Ding, W. Z. Huang, Z. X. Dong, Qualitative Theory of Differential Equations (in Chinese), Science Press, Beijing, 1992. English edition, American Mathematical Society, Providence, RI (1992).
    [26] W. Wang, C. Wang, S. Guo, On the walk matrix of the Dynkin graph Dn, Linear Algebra Appl., 653 (2022), 193–206. https://doi.org/10.1016/j.laa.2022.08.015 doi: 10.1016/j.laa.2022.08.015
    [27] R. I. Bogdanov, The versal deformations of a singular point on the plane in the case of zero eigenvalues, Trudy Sem. Petrovsk., 2 (1976), 37–65.
    [28] F. Takens, Forced oscillations and bifurcations. In: Applications of Global Analysis I, pp. 1–59. Comm. Math. Inst. Rijksuniv. Utrecht, Math. Inst. Rijksuniv. Utrecht, 1974.
    [29] S.-N. Chow, C. Z. Li, D. Wang, Normal Forms and Bifurcation of Planar Vector Fields, Cambridge University Press, Cambridge, 1994.
    [30] A. Korobeinikov, A Lyapunov function for Leslie-Gower predator-prey models, Appl. Math. Lett., 14 (2001), 697–699. https://doi.org/10.1016/S0893-9659(01)80029-X doi: 10.1016/S0893-9659(01)80029-X
  • This article has been cited by:

    1. Muhammad Aqib Abbasi, Periodic behavior and dynamical analysis of a prey–predator model incorporating the Allee effect and fear effect, 2024, 139, 2190-5444, 10.1140/epjp/s13360-024-04909-6
    2. Xinhao Huang, Lijuan Chen, Yue Xia, Fengde Chen, Dynamical Analysis of a Predator–Prey Model with Additive Allee Effect and Migration, 2023, 33, 0218-1274, 10.1142/S0218127423501791
    3. Yue Xia, Lijuan Chen, Vaibhava Srivastava, Rana D. Parshad, Stability and bifurcation analysis of a two-patch model with the Allee effect and dispersal, 2023, 20, 1551-0018, 19781, 10.3934/mbe.2023876
    4. A. M. Madian, M. F. Elettreby, M. M. A. El-sheikh, A. A. El-Gaber, Dynamical and Phase Analyses on Chaos Control of Discrete Predator–Prey Systems of Mixed Holling Types with Density-Dependent Birth Rates, 2024, 34, 0218-1274, 10.1142/S0218127424501712
    5. Mohamed Ch-Chaoui, Karima Mokni, A multi-parameter bifurcation analysis of a prey–predator model incorporating the prey Allee effect and predator-induced fear, 2025, 0924-090X, 10.1007/s11071-025-11063-w
    6. A.M. Madian, Rong Zheng, M.M.A. El-sheikh, M.F. Elettreby, A.A. El-Gaber, The influence of fear effect and phase of chaos control of discrete predator–prey system with mixed Holling types, 2025, 127, 11100168, 296, 10.1016/j.aej.2025.04.048
  • Reader Comments
  • © 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1938) PDF downloads(150) Cited by(6)

Figures and Tables

Figures(4)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog