# The Effects of Harvesting on the Dynamics of a Leslie–Gower Model

Authors: Jingli Xie; Hanyan Liu; Danfeng Luo
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

## Abstract

In this paper, we study a Leslie–Gower predator-prey model with harvesting effects. We carry out local bifurcation analysis and stability analysis. Under certain conditions, the model is shown to undergo a supercritical Hopf bifurcation resulting in a stable limit cycle. Numerical simulations are presented to illustrate our theoretic results.


---

## Body

## 1. Introduction

In this paper, we consider Leslie–Gower predator-prey model with harvesting effect,(1)dx1dτ=rx11−x1k−a1x1x2n+x1−e1x1,dx2dτ=sx21−a2x2n+x1−e2x2,with x10≥0 and x20≥0, where x1τ and x2τ are the prey and predator population densities, respectively, r,s,a1,a2,n,e1,e2>0, and τ is the time.Note thata2x2/n+x1 is Leslie–Gower term in which the carrying capacity of the predator’s environment is a linear function of the prey size x1/a2+n/a2. a1x1x2/n+x1 is the number of prey consumed by the predator in unit time which shows that when the number of the prey x1 is severe scarcity, and the predators can switch over to other populations as food. Constants rands are the intrinsic growth rate of the prey and predator, respectively, and e1ande2 denote the harvesting efforts for the prey and predator, respectively.Since the first prey-predator dynamical models which is the Lotka–Volterra model was built in the 1920s by Mathematician Lotka and Volterra, more and more researchers are interested in such issues, and they start from different angles to think the problem and many important results have been obtained [1–10]. In particular, in 2003, Aziz-Alaoui and Daher Okiye [11] considered the following Leslie–Gower predator-prey model:(2)dx1dτ=rx11−x1k−a1x1x2n1+x1,dx2dτ=sx21−a2x2n2+x1,where x1 is the numbers of prey and x2 is the numbers of predators. Existence and stability of the fixed points were studied by using the Lyapunov function. In 2006, Lin and Ho [12] discussed the local and global stability for system (2) by using Poincaré–Bendixson theorem and Dulac’s criterion.Harvesting is an effective way for humans to control the size of predators and prey so that the population has continued to develop healthily and produced good economic benefits [13–16]. Academically, researchers often only consider the harvesting of prey in order to control the size of the population. In 2010, Zhu and Lan [17] investigated the Leslie–Gower predator-prey systems:(3)dx1dτ=r1x11−x1k−ax1x2−h,dx2dτ=r2x21−x2bx1.In 2013, Gupta and Chandra [18] discussed the following Leslie–Gower predator-prey model with harvesting on the prey and the environment providing the same protection to both the predator and prey:(4)dx1dτ=rx11−x1k−a1x1x2n+x1−qEx1m1E+m2x1,dx2dτ=sx21−a2x2n+x1.For ecological balance and healthy economic development, for fisheries, wildlife resources, etc., we not only need to consider the harvesting of prey, but also the predator. Therefore, in this paper, we study Leslie–Gower predator-prey model (1) with harvesting on the prey and predator.In order to make the system dimensionless, we define new dependent variables byx=x1/k, y=a1x2/k and a new time variable by t=rτ. Moreover, let α=1/r, β=a2/a1, m=n/k, ρ=s/r, γ=e1/r, and δ=e2/ra1. System (1) becomes(5)dxdt=x1−x−αxym+x−γx,dydt=ρy1−βym+x−δy,where α,β,γ,ρ,δ,andm are all positive parameters.Lemma 1 (see19).
Ifa>0,b>0, andx˙≤xb−axα, whereαis a positive constant,t≥0, andx0>0, we have(6)xt≤ba1/α1+bx−α0a−1e−bαt−1/α.Lemma 2 (see [20, 21]).
Consider systemX˙=fX,α and suppose that fX0,α0=0, n×n Jacobian matrix J≡DfX0,α0 has a simple eigenvalue s=0 with eigenvector V, and the transpose of the Jacobian matrix JT has an eigenvector W to the eigenvalue s=0. Then, the system X˙=fX,α experiences a transcritical bifurcation at the equilibrium point X0 as the control parameter α passes through the bifurcation value α=α0 if the following conditions are satisfied:(7)WTDfαX0,α0=0,WTDfαX0,α0V≠0,WTD2fX0,α0V,V=0.The rest of this paper is organized as follows. In Section2, we study boundary of solutions. In Section 3, we discuss existence of equilibria points. In Section 4, we discuss stability of the equilibrium points.


## 2. Boundedness of Solutions

In this section, we prove that every solution of system (5) is positive and uniformly bounded with initial conditions x0=x0>0 and y0=y0>0. Denote ℝ+2=xt,yt∈ℝ2|xt>0,yt>0,t>0.Theorem 1.
Consider system (5). For any given initial conditionsx0,y0∈ℝ+2, the solutionxt,ytof system (5) exists and is unique and positive and ultimate bounded.Proof.
Let functionfxt,yt=xt1−xt−αxtyt/m+xt−γxt,ρyt1−βyt/m+xt−δyt. Obviously, function fxt,yt is continuous differentiable on x,y∈ℝ+2, so for any given the initial conditions x0,y0∈ℝ+2, the solution xt,yt of the system (5) exists and is unique. Furthermore, x−axis and y−axis are the solutions of the system (5); by the uniqueness of the solution, the solutions xt,yt of the system (5) with the initial value x0>0,y0>0 cannot cross with x−axis and y−axis.
Next, we show the solutionsxt,yt of system (5) with the initial value x0>0andy0>0 which is ultimate bounded.
From system (5), we have(8)dxdt≤x1−x.
Combining Lemma1, we have(9)xt≤1+x−10−1e−t−1≤M0,M0=max1,x0.
So, the following inequality is established:(10)dydt≤ρy1−βym+x≤ρy1−βym+1=yρ−ρβym+1.
By Lemma1, we have(11)yt≤M1,t≥0.where M1=maxm+M0/β,y0.


## 3. Existence of Equilibria

In order to find the equilibrium points of system (5), we let fxt,yt=0, i.e.,(12)x1−x−αxym+x−γx=0,ρy1−βym+x−δy=0.It is clear that equation (12) has a trivial solution E0≔0,0. Furthermore, by calculation, we find other solutions of equation (12):(13)1−γ,0,0,ρ−δmρβ,1−γρβ−ρ−δαρβ,ρ−δ1+m−γρβ−ρ−δαργ2.Remark 1.
Whenγ<1,δ<ρ and ρ−δα<1−γρβ, the position of each point in (13) is shown in Figure 1.Figure 1
The position of each point in (13) with ρ=0.1, α=1, δ=0.05, γ=0.05, β=5, and m=0.1.Therefore, we have the following results.wherex∞=1−γρβ−ρ−δα/ρβ,y∞=ρ−δ1+m−γρβ−ρ−δα/ρβ2.Theorem 2.
Consider system (5) admitsx−axial only andy−axial equilibria under following conditions.(i)
Thex−axial equilibrium,(14)E1≔K,0,withK=1−γ,is a boundary equilibrium of system (5) if and only if(15)1−γ>0.(ii)
They−axial equilibrium,(16)E2≔0,M,withM=ρ−δmρβ,is a boundary equilibrium of system (5) if and only if(17)ρ−δ>0.Theorem 3.
Consider system (5) admits a unique positive equilibrium,(18)E3≔x∞,y∞,if only if(19)1−γα>ρ−δρβ>0,Remark 2.
Existence regions for equilibrium points of system (12) is shown in Figure 2. Equilibrium pointsE1andE2exist, but positive equilibrium pointE3does not in region II and equilibrium pointsE1,E2, andE3coexist in region I. Moreover, if1−γ/α=ρ−δ/ρβ>0, then positive equilibriumE3becomes boundary equilibriumE2.Figure 2
Regions for equilibrium points of system (12).


## 4. Stability of of Equilibria

Theorem 4.
If1−γ>0andρ−δ>0, then the equilibriaE0andE1are unstable.Proof.
Firstly, we show the equilibriaE0 is a unstable equilibrium. The Jacobian matrix about E0 is given by(20)JE0=1−γ00ρ−δ,with trace 1−γ+ρ−δ and determinant 1−γρ−δ. From 1−γ>0 and ρ−δ>0, we have 1−γ+ρ−δ>0 and 1−γρ−δ>0. On the contrary, through direct calculation, we obtain 1−γ+ρ−δ2−41−γρ−δ=1−γ−ρ+δ2>0. Therefore, E0 is an unstable node point.
Secondly, we show the equilibriaE1 is a unstable equilibrium. The Jacobian matrix about E1 is given by(21)JE1=−1−γ−α1−γm+1−γ0ρ−δ,with determinant −1−γρ−δ<0. Therefore, E1 is a saddle point which is unstable.Theorem 5.
(1) If1−γ/α>ρ−δ/ρβ>0, then the boundary equilibriumE2is a unstable equilibrium. (2) If0<1−γ/α<ρ−δ/ρβ, then the boundary equilibriumE2is locally asymptotically stable.Proof.
The Jacobian matrix aboutE2 is given by(22)JE2=1−γ−αρ−δρβ0ρ−δ2ρβ−ρ−δ,with trace 1−γ−αρ−δ/ρβ−ρ−δ and determinant −ρ−δ1−γ−αρ−δ/ρβ.(1)
When1−γ/α>ρ−δ/ρβ>0, we have −ρ−δρβ1−γ−αρ−δ/ρβ<0, so E2 is a saddle point which is an unstable.(2)
When0<1−γ/α<ρ−δ/ρβ, we have 1−γ−αρ−δ/ρβ−ρ−δ=ρβ1−γ−αρ−δ/ρβ−ρ−δ<0 and −ρ−δρβ1−γ−αρ−δ/ρβ>0. Furthermore, a direct calculations gives(23)1−γ−αρ−δρβ−ρ−δ2+4ρ−δ1−γ−αρ−δρβ=1−γ−αρ−δρβ+ρ−δ2>0.Therefore,E2 is locally asymptotically stable.Theorem 6.
(1) If0<1−γ/α<ρ−δ/ρβandm<−2x∞2+1−γ−ρ+δx∞/x∞+ρ−δhold, then the positive equilibriumE3is unstable. (2) If0<1−γ/α<ρ−δ/ρβandm>−2x∞2+1−γ−ρ+δx∞/x∞+ρ−δhold, then the positive equilibriumE3is locally asymptotically stable. Furthermore, assume that there isα<ρβ; then, the positive equilibriumE3is globally asymptotically stable.Proof.
The Jacobian matrix aboutE3 is given by(24)JE3=x∞−1+1−γ−x∞m+x∞−αx∞m+x∞ρ−δ2ρβ−ρ−δ,with trace x∞−1+1−γ−x∞/m+x∞−ρ−δ and determinant −x∞ρ−δ−1+1−γ−x∞/m+x∞+αx∞/m+x∞ρ−δ2/ρβ.
A direct calculation gives(25)x∞−1+1−γ−x∞m+x∞−ρ−δ=x∞−m+1−γ−2x∞−mρ−δ−x∞ρ−δm+x∞=−2x∞2−mx∞+ρ−δ+1−γ−ρ+δx∞m+x∞.(1)
From conditions0<1−γ/α<ρ−δ/ρβ and m<−2x∞2+1−γ−ρ+δx∞/x∞+ρ−δ, we obtain(26)x∞−1+1−γ−x∞m+x∞−ρ−δ>0.Therefore, the positive equilibriumE3 is unstable.(2)
From conditions0<1−γ/α<ρ−δ/ρβ and m>−2x∞2+1−γ−ρ+δx∞/x∞+ρ−δ, we have(27)x∞−1+1−γ−x∞m+x∞−ρ−δ<0,−x∞ρ−δ−1+1−γ−x∞m+x∞+αx∞m+x∞ρ−δ2ρβ=−x∞ρ−δ−1+1−γm+x∞−x∞m+x∞−αρ−δρβm+x∞=x∞ρ−δ1−1−γm+x∞+x∞m+x∞+αρ−δρβm+x∞>x∞ρ−δ1+x∞m+x∞>0.So, the positive equilibriumE3 is locally asymptotically stable.Letβx,y=1/xy,Fx,y=dx/dt=x1−x−αxy/m+x−γx,Gx,y=dy/dt=ρy1−βy/m+x−δy. By calculations, we have(28)∂βx,yFx,y∂x+∂βx,yGx,y∂y=−1y+αm+x2−ρβxm+x≤−1y+αxm+x−ρβxm+x<−1y<0.Appling Dulac’s criterion and Theorem1, we know that the positive equilibrium E3 is globally asymptotically stable.


## 5. Local Bifurcation

hold, then system (5) undergoes a Hopf bifurcation with respect to bifurcation parameter m around the equilibrium point E3=x∞,y∞. Furthermore, the direction of the Hopf bifurcation is subcritical and the bifurcation periodic solutions are orbitally asymptotically stable if(41)k0+k1+k2+k3<0.Theorem 7.
If1−γ/α=ρ−δ/ρβ>0, then system (5) undergoes a transcritical bifurcation aroundE2.Proof.
From1−γ/α=ρ−δ/ρβ>0, we have(29)α=1−γρβρ−δ.
Therefore,(30)fαx,y=−xym+x0.
We have(31)fαx2,y2=00.
The eigenvectorsV and W associated to zero eigenvalues of matrixes J and JT, respectively, are(32)V=ρβρ−δ1,(33)W=10.
So,WTfαx2,y2=0.
In addition,Dfαx2,y2V equals(34)Dfαx2,y2V=∂−xy/m+x∂xv1+∂−xy/m+x∂yv20=−10.
Combining (33) and (34), we have(35)WTDfαx2,y2V=10−10=−1≠0.
Then, we computeD2fx2,y2,αV,V as follows:(36)D2fx2,y2,αV,V=∂2f1∂x2v1v1+2∂2f1∂x∂yv1v2+∂2f1∂y2v2v2∂2f2∂x2v1v1+2∂2f2∂x∂yv1v2+∂2f2∂y2v2v2=−2−2ρβm.
Combining (33) and (36), we have(37)WTD2fαx2,y2V,V=10−2−2ρβm=−2≠0.
Since all three conditions of Lemma2 are satisfied, system (5) undergoes a transcritical bifurcation around E2 if 1−γ/α=ρ−δ/ρβ>0.
Let(38)k0=ρ−δ8m+x∞2m−1m+x∞−3αmρβm+x∞2+2ρβα,k1=18mx∞m+x∞+m2ρ−δx∞m+x∞2−αm2ρ−δρβx∞m+x∞3,k2=18ρβρ−δm+x∞2m+x∞α−m+2ρβ1−ρβαx∞+αmρ−δm+x∞2m+x∞α−1ρβ−αρβρ−δm+x∞2m+x∞α−1ρβ2,k3=141+αmρ−δm+x∞2m+x∞α−1ρβ1+αmρ−δm+x∞2m+x∞α−1ρβ−αmρ−δm+x∞2m+x∞α−1ρβ2,and m is mH of (40).Theorem 8.
Assume that0<1−γ/α<ρ−δ/ρβis satisfied. If(39)ρβ1−γ<2α−ρβρ−δ,(40)mH=−2x∞2+1−γ−ρ+δx∞x∞+ρ−δ.The direction of the Hopf bifurcation is supercritical and the bifurcation periodic solutions are unstable if(42)k0+k1+k2+k3>0.Proof.
From (39), we have(43)2x∞<1−γ−ρ+δ.
So,mH>0.
The Jacobian matrix of system (5) evaluated at the point E3 is given by(44)JE3=x∞−1+1−γ−x∞m+x∞−αx∞m+x∞ρ−δ2ρβ−ρ−δ.
The tracez=TrJE3 and the determinant D=DetJE3 of Jacobian matrix JE3 are given by(45)z=−x∞+x∞1−γ−x∞m+x∞−ρ−δ=−2x∞2+1−γ−mx∞m+x∞−ρ−δ=−2x∞2+1−γ−ρ+δx∞−mx∞+ρ−δm+x∞,D=−x∞ρ−δ−1+1−γ−x∞m+x∞+ρ−δ2ρβαx∞m+x∞=x∞ρ−δρβm−ρ−δα+1−γρβρβm+x∞=x∞ρ−δ.
By0<1−γ/α<ρ−δ/ρβ and (40), we have z|mH=0 and D|mH>0.
In addition, we have(46)∂z∂m=−x∞+ρ−δm+x∞−−2x∞2+1−γ−ρ+δx∞−mx∞+ρ−δm+x∞2=x∞2−1−γx∞m+x∞2=−x∞1−γ−x∞m+x∞2=−αx∞ρ−δρβm+x∞2<0.
So,∂z/∂m|mH<0. Therefore, this guarantees the existence of Hopf bifurcation around E3.
We translate the equilibriumE3 to the origin by the translation x˜=x−x∞andy˜=y−y∞. For the sake of convenience, we still denote x˜ and y˜ by x and y, respectively. So, the system (5) becomes(47)x˙=x+x∞1−x−x∞−αx+x∞y+y∞m+x+x∞−γx+x∞,y˙=ρy+y∞1−βy+y∞m+x+x∞−δy+y∞.
Rewrite system (47) to(48)x˙y˙=JE3xy+f1x,y,mf2x,y,m,where(49)f1x,y,m=x+x∞1−x−x∞−αx+x∞y+y∞m+x+x∞−γx+x∞−x∞−1+1−γ−x∞m+x∞x+αx∞m+x∞y,f2x,y,m=ρy+y∞1−βy+y∞m+x+x∞−δy+y∞−ρ−δ2ρβx+ρ−δy.
Denote the eigenvalues ofJE3 by φ+iω with φ=z/2 and ω=4D−z2/2.
Define matrix(50)P=10NG,where G=ωm+x∞/αx∞ and N=m+x∞x∞−1+1−γ−x∞/m+x∞−φ/αx∞.
Obviously,(51)P−1=10−NG1G,and when m=mH, i.e., −2x∞2+1−γ−mHx∞=ρ−δx∞+mH, we have(52)ω0≔ωmH=ρ−δx∞,N0≔N|mH=mH+x∞−1+1−γ−x∞/mH+x∞α=mH+x∞ρ−δα,G0≔G|mH=ω0mH+x∞αx∞,N0G0=ρ−δx∞ω0=ω0,1G0=αx∞ω0mH+x∞.
By the transformation,(53)xy=Puv=uNu+Gv,system (48) becomes(54)u˙v˙=φ−ωωφuv+F1u,v,mF2u,v,m,where(55)F1u,v,m=A20u2+A11uv+A21u2v+A30u3+ou4,u4v,F2u,v,m=B20u2+B11uv+B02v2+B21u2v+ou4,u4v,with(56)A20=−1−αmNx∞+Nm−y∞m+x∞3,A11=−αGmm+x∞2,A21=αGmm+x3,A30=αmNx∞+Nm−y∞m+x∞4,B20=NG1+αmN∞+Nm−y∞m+x∞3−1GρβN2x∞2+2N2x∞m−2Nx∞y∞+Nm−y∞2m+x∞3,B11=αNmm+x∞2−2ρβNx∞+Nm−y∞m+x∞2,B02=−ρβGm+x∞,B21=2ρβNx∞+Nm−y∞−αNmm+x∞3.
In order to determine the stability of the periodic solution, we need to calculate the sign of the coefficientbmH, which is given by(57)bmH≔116Fuuu1+Fuvv1+Fuuv2+Fvvv2+116ωmHFuv1Fuu1+Fvv1−Fuv2Fuu2+Fvv2−Fuu1Fuu2+Fvv1Fvv2.where all partial derivatives are evaluated at the bifurcation point 0,0,mH.
CombiningF1u,v, F2u,v, ω0, G0, and N0, we have(58)Fuuu10,0=6mρ−δm+x∞2−6αmρ−δρβm+x∞3,Fuvv10,0=0,Fuv10,0=−mω0x∞m+x∞,Fuu10,0=−2−2αmρ−δαm+x∞−ρ−δρβm+x∞2,Fvv10,0=0,Fvvv20,0=0,Fuuv20,0=4ρβρ−δαm+x∞−2m+4ρ−δm+x∞2,Fuv20,0=m+2ρ−δm+x∞−2ρβρ−δα,Fuu20,0=2ω01+αmρ−δm+x∞2m+x∞α−1ρβ−αρβρ−δm+x∞2m+x∞α−1ρβ2,Fvv20,0=−2ρβω0αx∞,where m=mH and(59)bmH=ρ−δ8m+x∞2m−1m+x∞−3αmρβm+x∞2+2ρβα+18mx∞m+x∞+m2ρ−δx∞m+x∞2−αm2ρ−δρβx∞m+x∞3+ρβρ−δm+x∞2m+x∞α−m+2ρβ×1−ρβαx∞+αmρ−δm+x∞2m+x∞α−1ρβ−αρβρ−δm+x∞2m+x∞α−1ρβ2+21+αmρ−δm+x∞2m+x∞α−1ρβ1+αmρ−δm+x∞2m+x∞α−1ρβ−αmρ−δm+x∞2m+x∞α−1ρβ2=k0+k1+k2+k3.
Combining (41), we have bmH<0. Therefore, according to Poincare–Andronow’s Hopf bifurcation theory, we have the direction of the Hopf bifurcation is subcritical and the bifurcation periodic solutions are orbitally asymptotically stable.
In addition, combining (42), we have bmH>0. Therefore, according to Poincare–Andronow’s Hopf bifurcation theory, we have the direction of the Hopf bifurcation is supercritical and the bifurcation periodic solutions are unstable.


## 6. Numerical Illustrations

In this section, we perform numerical simulations about system (5). Figure 3 shows that E0 is an unstable node point, E1 is a saddle point, E3 does not exist, and E2 is asymptotically stable and every orbit tends to it. Figure 4 shows that E0 is an unstable node point, E1 is a saddle point, E2 is unstable, E3 is unstable, and there is a limit cycle around E3 to which every orbit tends. Figure 5 shows that E0 is an unstable node point, E1 is unstable and E2 is also unstable, but E3 is asymptotically stable and every orbit approaches this equilibrium.Figure 3
Phase portrait for system (12) with α=β=ρ=m=1, γ=0.5, and δ=0.1.Figure 4
Phase portrait for system (12) with α=ρ=1,β=0.5, m=0.01, γ=0.5, and δ=0.8.Figure 5
Phase portrait for system (12) with α=ρ=1,β=0.5, m=0.3, γ=0.5, and δ=0.8.



---

*Source: 5520758-2021-05-07.xml*