# Michaelis-Menten-Type Prey Harvesting in Discrete Modified Leslie-Gower Predator-Prey Model

Authors: M. Saqib Khan; Mujahid Abbas; Ebenezer Bonyah; Hengxiao Qi
Journal: Journal of Function Spaces
 (2022)

---

## Abstract

A discrete-time Michaelis-Menten-type prey harvesting is discussed in this paper, in the modified Leslie-Gower predator-prey model. Detailed analysis of the topology of nonnegative interior fixed points is given, including their existence and stability dynamics. Also, the conditions for the existence of flip and Neimark-Sacker bifurcations are derived by using the center manifold theorem and bifurcation theory. The numerical simulations are provided, using a computer package, to illustrate the consistency of theoretical results.


---

## Body

## 1. Introduction

The basic interaction between various species on this planet is predation and is widely present in nature. The existence of most species in our ecosystem is based on predator-prey relationships, which makes the rich biodiversity of complex ecosystems possible [1]. These complexities can aptly be explained using mathematical models and the related dynamics. As Clark has observed in [2], the dynamics of the predator-prey models also provide methods to optimally manage renewable resources, apart from coexistence conditions of predators and preys. For this ecosystem, the most basic model was proposed by Lotka [3] and Volterra [4], in early 20th century. This model captured the oscillating behavior in populations of a predator and its prey. The simplicity of this model, however, is unable to address most real-world scenarios. To rectify this, many modifications of Lotka-Volterra type model have been proposed by researchers. In [5], Holling proposed various types of functional responses to better understand the components of predator-prey interactions.An alternative was proposed by Leslie, in [6], in which the number of prey and the carrying capacity of the predator’s environment was proportional. Leslie focused on the fact that the increasing capacity of both predator and prey cannot be unbounded. This fact was not incorporated in the Lotka-Volterra model. This model was shown, in [7], to possess globally asymptotically stable and unique positive equilibrium, for any permissible parameters. Then, May in [8] showed its stability with Holling functional responses. L. Chen and F. Chen proved the global stability of the unique interior equilibrium for a Leslie-Gower predator-prey model with feedback controls [9]. Uniqueness of limit cycles and the Hopf bifurcation for this model was discussed in [10, 11].The primary drawback of this model is that if the prey population is at low densities, it is not possible for the predator to switch to alternative food source [12]. Aziz-Alaoui and Daher [13] rectified this scenario by adding a provisional alternative food source parameter. This model is generally called the modified Leslie-Gower predator-prey model. Caughley [14] modelled the biological control of the prickly-pear cactus by the moth Cactoblastis cactorum, using this system. In [15, 16], the authors used this system to model the predator-prey mite outbreak interactions on fruit trees. Much work has been done on this model. For further details, look at [17–23] and the references therein. In [24], W. C. Allee introduced what was altered termed as the Allee effects, to the Leslie-Gower type models. The stability dynamics and bifurcation analysis of the predator-prey systems subject to Allee effect is discussed in [25–28].Harvesting of different biological species is necessary for a variety of reasons. The ever-increasing use of natural resources, either due to the increase in human population or in the name of economic progress, has caused the ecological system of this planet to trip from the equilibrium it was in, for thousands of years. A common outcome of these practices is the extinction of species from the face of this planet. Harvesting of such species is one way to minimize the damage. Harvesting for species is also practices for purely business reasons. But there is a silver lining in this strategy. The more harvesting is done, the more predator species are likely to become extinct [29]. Thus, the need to reinforce scientific management of harvesting is much needed. Clark [2] discussed the problem of combined or nonselective harvesting of two ecologically independent populations obeying the logistic law of growth. Multispecies harvesting models are studied in detail by Chaudhuri [30, 31], Mesterton-Gibbons [32], Kar and Chaudhuri [33, 34], etc. Nonselective harvesting model of a prey-predator fishery is studied in detail by Chaudhuri and Ray [35], Kar et al. [36], etc.Michaelis-Menten-type predator-prey model, alternatively known as the ratio-dependent predator-prey model, uses the idea that the ratio of prey to predator abundance directly affects the per capita predator growth rate. The effectiveness of this model is backed by a plethora of experimental and observational data, for the predator which must compete for prey. For more details, see ([37–43]). Michaelis-Menten-type predator-prey model has been rigorously studied in ([44–49]) and references therein. This investigation discovered existence of rich dynamics, like stable limit cycle, multiple attractors, and deterministic extinction, to name a few. Existence of hyperbolic, parabolic, and elliptic orbits was revealed near the origin, as well as combinations of such orbits, for various parametric values.The continuous predator-prey models have been successfully used to obtain desired results, and the disadvantages of using continuous systems are also quite apparent. By definition, the continuous systems require the subject species have continuous and overlapping generations. This is not generally true, for example, salmon, which have an annual spawning season and are born at the same time each year. For these types of populations which have nonoverlapping generation characteristics, the discrete time models are more descriptive and suitable than the continuous models [50]. Also, the discrete models should generate richer and reality-based dynamics compared to the continuous time models [51]. In addition, since many continuous models cannot be solved analytically, using difference equations for approximation and finding solution is much practical way to approach the problem. In population biology and complex ecosystems, discrete-time models are used to examine the taxonomic group of organisms and species with the passage of time. These models are best to describe the chaotic behavior of nonlinear dynamics [52, 53].Gupta and Chandra, in [54], proposed and studied modified Leslie-Gower predator-prey model with Michaelis-Menten-type prey harvesting in prey. They proved the permanence and stability and discuss different bifurcations of this model. Their model guaranteed a feasible upper bound of the rate of harvesting for the coexistence of the species. Their model is given in [54] as
(1)dxdt=x1−x−αym+x−hc+x,(2)dydt=ρy1−βym+x,with the initial conditions x0>0, y0>0. Here, xt and yt represent prey and predator densities, respectively, at time t. The details of the parameters are given in their paper. For our purposes, all parameters are positive.Another possible way to understand the complex problem of interaction between prey and predator is by using discrete models. We use the forward Euler method to discretize the above system. The discrete counterpart of (1) can be given as
(3)xn+1=xn+kxn1−xn−αynm+xn−hc+xn,(4)yn+1=yn+kρyn1−βynm+xn.In the remainder of this paper, we present the discrete counterpart of the Michaelis-Menten-type prey harvesting in modified Leslie-Gower predator-prey model, ((4)). Our aim is to illustrate the detailed mathematical analysis of the topology of nonnegative interior fixed points, including their existence in Section 2 and stability dynamics in Section 3. The conditions for the existence of flip and Neimark-Sacker bifurcations will be derived in Section 4, using the center manifold theorem and bifurcation theory. The theory will be verified through numerical examples in Section 5. Finally, we will conclude in Section 6.


## 2. The Fixed Points

In order to find the fixed points, we need to solve the following system.(5)x1−x−αym+x−hc+x=0,(6)y1−βym+x=0.The system has four fixed points on the boundary, namely,(7)E1,E2,E3,4=0,0,0,mβ,1−c±c+12−4h2,0.For positive equilibrium points, we have the following cases.(i)
The system has no positive stationary point, if(8)α−βc+12<4β2h(ii)
Ifα−βc+12=4β2h, then the system has a unique positive stationary point, E∗=x∗,m+x∗/β, where
(9)x∗=121−c−αβ(iii)
Ifα−βc+12>4β2h, the system will have two positive stationary points E5=xm,m+xm/β and E6=xp,m+xp/β, where
(10)xm=−α+β1−c−α−βc+12−4β2h2βxp=−α+β1−c+α−βc+12−4β2h2βAs can be seen in Figure1, the positive interior points may vanish or merge into one another, depending on the values of parameters.Figure 1
The purple graph is forα=0.3,h=0.1,m=0.1,c=0.05,β=1,ρ=2.6, which is showing two positive fixed points. The brown dashed graph is for α=0.42,h=0.1,m=0.1,c=0.05,β=1,ρ=2.6, which is showing unique positive fixed point. The black dotted graph is for α=0.5,h=0.1,m=0.1,c=0.05,β=1,ρ=2.6, which is showing no positive fixed point.


## 3. The Topological Dynamics of the Fixed Points

The respective Jacobian matrices, evaluated at each fixed point, are(11)J1=1+k1−hc001+kρJ2=1+k1−hc−αβ0kρβ1−kρJ3,4=1+kx−1+1−xc+x−αβ−kxαm+x01+kρJ∗,5,6=1+kA1−kA2kA31−kA4,where
(12)A1=xhc+x2+αym+x2−1,A2=αxm+x,A3=ρym+x,A4=ρ.In order to find the topology and properties of the aforementioned fixed points, the following lemma will be used, which is omnipresent in the text books related to discrete dynamical systems.Lemma 1.
Letti=TrJi, deti=Ji, and piz=z2−tiz+deti, i∈1,2,3,4,∗,5,6. Suppose z1 and z2 are the roots of pz. Then,
(i)
z1,2<1 if and only if p−1>0 and det<1(ii)
z1,2>1 if and only if p−1>0 and det>1(iii)
z1<1 and z2>1 (or vice versa) if and only if p−1<0(iv)
z1=−1 and z2≠1 if and only if p−1=0 and t≠−2,0(v)
z1,2 are complex, and z1,2=1 if and only if 4det−t2>0 and det=1Ifz1 and z2 are the eigenvalues of 2×2 Jacobian matrix, then,Lemma 2.
A fixed pointEi, i∈1,2,3,4,∗,5,6 is called
(i)
Sink ifz1<1 and z2<1, so sink is locally asymptotically stable(ii)
Source ifz1>1 and z2>1, so source is locally unstable(iii)
Saddle ifz1<1 and z2>1 (or z1>1 and z2<1)(iv)
Nonhyperbolic if eitherz1=1 or z2=1Using the above lemma, we can easily see thatE1 is saddle if c<h and unstable if c>h. If h=c, then it may undergo transcritical or fold bifurcation, and if h=c1+2/k, it undergoes period-doubling bifurcation. Similarly, the fixed representing the prey free scenario, E2, is saddle if h<c1−α/β and stable if h>c1−α/β. If h=c1−α/β, then it may undergo transcritical or fold bifurcation, and if h=c1−α/β+2/k, it undergoes period-doubling bifurcation. The fixed points, representing the extinction of predator, E3,4, are saddle if h<c+x¯2 and unstable if h>c+x¯2. If h=c+x¯2, then these fixed points may undergo transcritical or fold bifurcation, and if h=1−2/kx¯c+x¯2, the fixed points undergo period-doubling bifurcation. For all these boundary fixed points, Neimark-Sacker bifurcation is not possible.We are much more interested in the behavior of interior fixed points, which are much more richer in comparison and for practical purposes, usable.Theorem 3.
For the fixed pointsE∗, E5, and E6, if β≥A2/A1 then let ρ>0, and if β<A2/A1, then let ρ∉ρ1,ρ2, where
(13)ρ1=−βA1−2A2−2A2A2−A1ββ,ρ2=−βA1−2A2+2A2A2−A1ββ.
If(14)k1=A1−A4−A1+A42−4A2A3A1A4−A2A3,k2=A1−A4+A1+A42−4A2A3A1A4−A2A3.
Then,(1)
The fixed points are locally asymptotically stable ifk∈0,k1(2)
The fixed points are locally unstable ifk∈k2,∞(3)
The fixed points are saddle ifk∈k1,k2(4)
The fixed points are nonhyperbolic with eigenvaluesz1=1 and z2≠1, if and only if h=c+x¯2(5)
The fixed points are nonhyperbolic with eigenvaluesz1=−1 and z2≠1, if k∈k1,k2 and k∉2/A4−A1,4/A4−A1(6)
The fixed points are nonhyperbolic, having complex conjugate eigenvalues, withz1=1=z2 if and only if k=A1−A4/A1A4−A2A3 and k∈0,4/A4−A1Proof.
For any positive interior fixed point,E∗, E5, or E6, the characteristic polynomial of the Jacobian matrix at the given fixed point is
(15)pz=z2−tz+det,where t=2+A1−A4k and det=A2A3−A1A4k2+A1−A4k+1. Then, we can use Lemma 1, since for any k>0,
(16)p1=1−t+det=A2A3−A1A4k2.
Thus,p1=0 if h/c+x2−1=0, or
(17)h=c+x2.
Also,(18)p−1=A2A3−A1A4k2+2A1−A4k+4,which shows that
(19)p−1>0k∈0,k1∪k2,∞,p−1<0k∈k1,k2,p−1=0k∈k1,k2.
Also,(20)det>1k∈A1−A4A1A4−A2A3,∞,det<1k∈0,A1−A4A1A4−A2A3,det=1k=A1−A4A1A4−A2A3.
Then, using Lemmas1 and 2, we have the following results.
(1)
Each fixed point is a sink if and only if(21)k∈0,k1∪k2,∞∩0,A1−A4A1A4−A2A3=0,k1(2)
Each fixed point is a source if and only if(22)k∈0,k1∪k2,∞∩A1−A4A1A4−A2A3,∞=k2,∞(3)
Each fixed point is saddle if and only ifk∈k1,k2(4)
The Jacobian matrix at the given fixed points has eigenvaluesz1=1 and z2≠1 if and only if h=c+x¯2(5)
The Jacobian matrix at the given fixed points has eigenvaluesz1=−1 and z2≠1, if and only if k=k1 or k=k2 and h≠2/A4−A1 and h≠4/A4−A1(6)
The eigenvalues are complex conjugates withz1=1=z2 if and only if k=A1−A4/A1A4−A2A3 and k∈0,4/A4−A1Remark 4.
As used in [55, 56], we know that due to the jury condition, a transcritical or fold bifurcation is obtained when pj1=0, j∈∗,5,6. As shown above, this is only true if h=c+x¯2. For the fixed point E∗, this condition is equivalent to h=1/4α/β−c−12, which we gets automatically, since in this scenario, α/β−c−12−4h=0. Thus, we always have p∗1=0 for the fixed point E∗, i.e., a transcritical or fold bifurcation is always obtained for any parametric values, at the fixed point E∗.


## 4. Bifurcations

For positive fixed points,E5 and E6, define
(23)ΩPD=α,β,h,m,c,ρ,k0∈ℝ+7:k0=k1ork2,k0≠2A4−A1,k0≠4A4−A1,ρ∉ρ1,ρ2,ΩNS=α,β,h,m,c,ρ,k0∈ℝ+7:k0=A1−A4A1A4−A2A3,k∈0,4A4−A1.The procedure is identical for finding the normal form of the bifurcations forE5 and E6, as illustrated below.


### 4.1. Period-Doubling Bifurcation

Letα,β,h,m,c,ρ,k0∈ΩPD and let K denote the perturbation parameter for the mapping (3), i.e., K⋘1. Then, variation of parameters in small neighborhood of ΩPD gives emergence of period doubling bifurcation. Let k0 be the bifurcation parameter. The perturbed mapping is
(24)xy⟶x+k+Kx1−x−αym+x−hc+xy+k+Kρy1−βym+x.DefineX=x−x¯ and Y=y−y¯, to translate the fixed point to 0,0. Here, x¯,y¯ is the positive interior fixed point. Then, the above mapping is converted to
(25)XY⟶X+k+KX+x¯1−X+x¯−αY+y¯m+X+x¯−hc+X+x¯Y+k+KρY+y¯1−βY+y¯m+X+x¯.After Taylor series expansion, the above mapping can be written as(26)XY⟶a11a12a21a22XY+fPDX,Y,KgPDX,Y,K,where a11=1+kA1, a12=−kA2, a21=kA3, and a22=1−kA4, and
(27)fPDX,Y,K=kb1X2+kb2XY+A1KX−A2KY+kb3Y2+kb4X3+kb5X2Y+b1KX2+b2KXY+b3KY2+kb6XY2+kb7Y3+OX+Y+K4,gPDX,Y,K=kb8X2+kb9XY+A3KX−A4KY+kb10Y2+kb11X3+kb12X2Y+b8KX2+b9KXY+b10KY2+kb13XY2+kb14Y3+OX+Y+K4,and
(28)b1=chc+x¯3+αmy¯m+x¯3−1,b2=−αmm+x¯2,b3=0b4=−chc+x¯4−αmy¯m+x¯4,b5=αmm+x¯3,b6=0,b7=0,b8=−ρy¯m+x¯2,b9=2ρm+x¯,b10=−βρm+x¯,b11=ρy¯m+x¯3,b12=−2ρm+x¯2,b13=βρm+x¯2,b14=0.To convert the system (43) in the normal form, let
(29)XY=TPDuv=a12a12z1−a11z2−a22uv=a12a12−1−a11z2−a22uv,where TPD is an invertible matrix. Thus, we can convert the above system, such that we have 0,0 as fixed points, as shown below.
(30)uv⟶−100z2uv+TPD−1fPDu,v,KgPDu,v,K⟶−100z2uv+f~PDu,v,Kg~PDu,v,K,where
(31)f~PDu,v,K=d1u2+d2uv+d3Ku+d4Kv+d5v2+d6u3+d7u2v+d8Ku2+d9Kuv+d10Kv2+d11uv2+d12v3+Ou+v+K4,g~PDu,v,K=d13u2+d14uv+d15Ku+d16Kv+d17v2+d18u3+d19u2v+d20Ku2+d21Kuv+d22Kv2+d23uv2+d24v3+Ou+v+K4,and
(32)d1=c1ka12a11+1b9−b2c2+a122b1c2−b8+a11+12−b10,d2=c1−k−a12a11+a22−z2+1b9−b2c2+2a122b8−b1c2+2a11+1b10a22−z2,d3=c1a12A1c2+A3+a11+1A2c2−A4,d4=c1a22−z2A2c2−A4+a12A1c2+A3,d5=c1ka12a22−z2b9−b2c2+a122b1c2−b8−b10a22−z22,d6=c1ka123b4c2−b11+a11+1a122b12−b5c2−a11+12a12b13,d7=c1ka1222a11+a22−z2+2b12−b5c2−3a123b11−b4c2−a11+1a12b13a11+2a22−2z2+1,d8=c1a12a11+1b9−b2c2+a122b1c2−b8+a11+12−b10,d9=−c1−a12a11+a22−z2+1b9−b2c2+2a122b8−b1c2+2a11+1b10a22−z2,d10=c1a12a22−z2b9−b2c2+a122b1c2−b8−b10a22−z22,d11=c1ka122a11+2a22−2z2+1b12−b5c2−3a123b11−b4c2−a12b13a22−z22a11+a22−z2+2,d12=c1ka122a22−z2b12−b5c2+a123b4c2−b11−a12b13a22−z22,d13=c1ka12a11+1−b2c3−b9+a122b1c3+b8+a11+12b10,d14=c1−ka12a11+a22−z2+1b2c3+b9−2a122b1c3+b8−2a11+1b10a22−z2,d15=c1a12A1c3+A3+a11+1A2c3−A4,d16=c1a22−z2A2c3−A4+a12A1c3+A3,d17=c1k−a12a22−z2b2c3+b9+a122b1c3+b8+b10a22−z22,d18=c1ka123b4c3+b11−a11+1a122b5c3+b12+a11+12a12b13,d19=c1k−a1222a11+a22−z2+2b5c3+b12+3a123b4c3+b11+a11+1a12b13a11+2a22−2z2+1,d20=c1−a12a11+1b2c3+b9+a122b1c3+b8+a11+12b10,d21=−c1a12a11+a22−z2+1b2c3+b9−2a122b1c3+b8−2a11+1b10a22−z2,d22=c1−a12a22−z2b2c3+b9+a122b1c3+b8+b10a22−z22,d23=c1k−a122a11+2a22−2z2+1b5c3+b12+3a123b4c3+b11+a12b13a22−z22a11+a22−z2+2,d24=c1k−a122a22−z2b5c3+b12+a123b4c3+b11+a12b13a22−z22andX=a12u+v,Y=−1+a11u+z2−a22v.In order to implement the center manifold theorem, we assume thatMC be the center manifold of (30), evaluated at 0,0 in a small neighborhood of K=0. It can be approximated as follows.
(33)MC=K,u,v:v=GPDK,u,u<δ1,K<δ2,GPD0,0=0,DGPD0,0=0,where
(34)FPDK,GPDK,u=GPDK,z1u+f~PDK,u,GPDK,u−z2GPDK,u−g~PDK,u,GPDK,u=GPDK,−u+f~PDK,u,GPDK,u−z2GPDK,u−g~PDK,u,GPDK,u.We need to findGPDK,u such that FPDK,GPDK,u=0. Since we need the conditions in MC to be satisfied, we may assume
(35)GPDK,u=m1u2+m2Ku+m3K2+OK+u3.Thus, we can write(36)FPDK,GPDK,u=m11−z2−d13u2−d15+m2z2+1Ku+m31−z2K2+OK+u3.By comparing coefficients in the equationFPDK,GPDK,u=0, we get
(37)m1=d131−z2,m2=−d151+z2,m3=0,and
(38)GPDK,u=d131−z2u2−d151+z2Ku+OK+u3.The dynamics restricted toMC are given locally by the map
(39)FC:u↦z1u+f~PDK,u,GPDK,u=−u+d1u2+d3Ku−d4d151+z2K2u+d8−d4d131−z2−d2d151+z2Ku2+d6+d2d131−z2u3+OK+u4.Furthermore, define(40)Y1=∂2FC∂K∂u+12∂FC∂K∂2FC∂u20,0=d3≠0,Y2=16∂3FC∂u3+12∂2FC∂u220,0=d12+d6+d2d131−z2≠0.Theorem 5.
IfY1≠0 and Y2≠0, then system undergoes period-doubling bifurcation at the unique positive equilibrium point when parameter k varies in small neighborhood of k1 or k2. Moreover, the period-two orbits that bifurcate from positive equilibrium are stable if Y2>0 and unstable if Y2<0.Remark 6.
Since we have already found the center manifold for our system in (39), we can easily show the existence of transcritical or fold bifurcation, when one of the eigenvalue equals 1 and the other not on the unit circle.


### 4.2. Neimark-Sacker Bifurcation

Letα,β,h,m,c,ρ,k0∈ΩNS, i.e., k0=A1−A4/A1A4−A2A3. Then, variation of α,β,h,m,c,ρ,k in small neighborhood of ΩNS gives emergence to Neimark-Sacker bifurcation. Let K denote the perturbation parameter for the mapping (3), i.e., K1. The perturbed mapping is
(41)xy⟶x+k+Kx1−x−αym+x−hc+xy+k+Kρy1−βym+x.DefineX=x−x¯ and Y=y−y¯, to translate the fixed point to 0,0. Here, x¯,y¯ is the positive interior fixed point. Then, the above mapping is converted to
(42)XY⟶X+k+KX+x¯1−X+x¯−αY+y¯m+X+x¯−hc+X+x¯Y+k+KρY+y¯1−βY+y¯m+X+x¯.After Taylor series expansion, the above mapping can be written as(43)XY⟶c11c12c21c22XY+fNSX,YgNSX,Y,where c11=1+k+KA1, c12=−k+KA2, c21=k+KA3, and c22=1−k+KA4, and
(44)fNSX,Y=k+Kb1X2+k+Kb2XY+k+Kb3Y2+kb4X3+kb5X2Y+kb6XY2+kb7Y3+OX+Y4,gNSX,Y=k+Kb8X2+k+Kb9XY+k+Kb10Y2+kb11X3+kb12X2Y+kb13XY2+kb14Y3+OX+Y4,and bi,i∈1,2,⋯,14 are defined as above. Let the characteristic polynomial of matrix in (43) be
(45)pNSz=z2−PKz+QK,where
(46)PH=c11+c22,QH=c11c22−c21c12.Sinceα,β,h,m,c,ρ,k∈ΩNS, therefore the roots of (45) are complex conjugate z1, z2, and z1=z2=1, and it follows that z1,2=PH/2±1/24QH−PH2 and 1=z1=z2=QH. Note that PK2−4QK<0, which implies that PK2<4QK or PK∈−2,2. Also,
(47)dz1,2dKK=0=12ρ−x¯hc+x¯2+αy¯m+x¯2−1.Thus, ifρ≠x¯h/x¯+c2+αy¯/x¯+m2−1, then
(48)dz1,2dKK=0≠0.In order to ensure that the roots of the characteristic polynomial do not lie in the intersection of unit circle of coordinate axis whenK=0, we need to check that z1,2m≠1, for m=1,2,3,4 at K=0. This is equivalent to checking P0≠−2,0,1,2. Since α,β,h,m,c,ρ,k∈ΩNS, we already know that P0≠−2,0,2. Finally,
(49)P0≠1⇔k≠1A4−A1.Thus, the roots of (45) do not lie in the intersection of unit circle of coordinate axis when K=0. Now, to convert the system (43) to the normal form, when K=0, let R=P0/2 and S=1/24Q0−P02 and let
(50)TNS=01SR.Here,TNS is an invertible matrix. Consider the transformation
(51)XY=TNSuv.The system (43) can be written as
(52)uv⟶R−SSRuv+f~NSu,vg~NSu,v,where
(53)f~NSu,v=e1u2+e2uv+e3v2+e4u3+e5u2v+e6uv2+e7v3+Ou+v4,g~NSu,v=e8u2+e8uv+e10v2+e11u3+e12u2v+e13uv2+e14v3+Ou+v4,and
(54)e1=k+K−Rb2S+b1+Sb10S+b9+b8S,e2=b2−R+2b10S+b9Rk+KS,e4=k−Rb5S+b4+Sb13S+b12+b11S,e3=b10R2k+KS,e5=kR2b13S+b12−b5RS,e6=b13kR2S,e7=0,e8=b2S+b1k+K,e9=b2Rk+K,e10=0,e11=kb5S+b4,e12=b5kR,e13=0,e14=0,X=v,Y=Su+Rv.For the map (52) to undergo Neimark-Sacker bifurcation, the following quantity must not zero [22]:
(55)MNS=−Reϖ20ϖ111−2z1z221−z1−12ϖ112−ϖ022+Reϖ21z2K=0,where
(56)ϖ20=18∂2f¯NS∂u2−∂2f¯NS∂v2+2∂2g¯NS∂u∂v+ι∂2g¯NS∂u2−∂2g¯NS∂v2−2∂2f¯NS∂u∂v,ϖ11=14∂2f¯NS∂u2+∂2f¯NS∂v2+ι∂2g¯NS∂u2+∂2g¯NS∂v2,ϖ02=18∂2f¯NS∂u2−∂2f¯NS∂v2−2∂2g¯NS∂u∂v+ι∂2g¯NS∂u2−∂2g¯NS∂v2+2∂2f¯NS∂u∂v,ϖ21=116∂3f¯NS∂u3+∂3f¯NS∂u∂v2+∂3g¯NS∂u2∂v+∂3g¯NS∂v3+ι∂3g¯NS∂u3+∂3g¯NS∂u∂v2−∂3f¯NS∂u2∂v−∂3f¯NS∂v3.Using the result in [57] and the above analysis, we have the following result.Theorem 7.
Ifk≠1/A4−A1 and MNS defined in (55) is nonzero, then the model (5) undergoes Neimark-Sacker bifurcation at the equilibrium point x¯,y¯, provided that parameter k changes in the small neighborhood of k0=A1−A4/A1A4−A2A3 and α,β,h,m,c,ρ,k0∈ΩNS. Moreover, if MNS<0 (resp., MNS>0), then an attracting (resp., repelling) invariant closed curve bifurcate from the fixed point x¯,y¯ for k>k0 (resp., k<k0).


## 4.1. Period-Doubling Bifurcation

Letα,β,h,m,c,ρ,k0∈ΩPD and let K denote the perturbation parameter for the mapping (3), i.e., K⋘1. Then, variation of parameters in small neighborhood of ΩPD gives emergence of period doubling bifurcation. Let k0 be the bifurcation parameter. The perturbed mapping is
(24)xy⟶x+k+Kx1−x−αym+x−hc+xy+k+Kρy1−βym+x.DefineX=x−x¯ and Y=y−y¯, to translate the fixed point to 0,0. Here, x¯,y¯ is the positive interior fixed point. Then, the above mapping is converted to
(25)XY⟶X+k+KX+x¯1−X+x¯−αY+y¯m+X+x¯−hc+X+x¯Y+k+KρY+y¯1−βY+y¯m+X+x¯.After Taylor series expansion, the above mapping can be written as(26)XY⟶a11a12a21a22XY+fPDX,Y,KgPDX,Y,K,where a11=1+kA1, a12=−kA2, a21=kA3, and a22=1−kA4, and
(27)fPDX,Y,K=kb1X2+kb2XY+A1KX−A2KY+kb3Y2+kb4X3+kb5X2Y+b1KX2+b2KXY+b3KY2+kb6XY2+kb7Y3+OX+Y+K4,gPDX,Y,K=kb8X2+kb9XY+A3KX−A4KY+kb10Y2+kb11X3+kb12X2Y+b8KX2+b9KXY+b10KY2+kb13XY2+kb14Y3+OX+Y+K4,and
(28)b1=chc+x¯3+αmy¯m+x¯3−1,b2=−αmm+x¯2,b3=0b4=−chc+x¯4−αmy¯m+x¯4,b5=αmm+x¯3,b6=0,b7=0,b8=−ρy¯m+x¯2,b9=2ρm+x¯,b10=−βρm+x¯,b11=ρy¯m+x¯3,b12=−2ρm+x¯2,b13=βρm+x¯2,b14=0.To convert the system (43) in the normal form, let
(29)XY=TPDuv=a12a12z1−a11z2−a22uv=a12a12−1−a11z2−a22uv,where TPD is an invertible matrix. Thus, we can convert the above system, such that we have 0,0 as fixed points, as shown below.
(30)uv⟶−100z2uv+TPD−1fPDu,v,KgPDu,v,K⟶−100z2uv+f~PDu,v,Kg~PDu,v,K,where
(31)f~PDu,v,K=d1u2+d2uv+d3Ku+d4Kv+d5v2+d6u3+d7u2v+d8Ku2+d9Kuv+d10Kv2+d11uv2+d12v3+Ou+v+K4,g~PDu,v,K=d13u2+d14uv+d15Ku+d16Kv+d17v2+d18u3+d19u2v+d20Ku2+d21Kuv+d22Kv2+d23uv2+d24v3+Ou+v+K4,and
(32)d1=c1ka12a11+1b9−b2c2+a122b1c2−b8+a11+12−b10,d2=c1−k−a12a11+a22−z2+1b9−b2c2+2a122b8−b1c2+2a11+1b10a22−z2,d3=c1a12A1c2+A3+a11+1A2c2−A4,d4=c1a22−z2A2c2−A4+a12A1c2+A3,d5=c1ka12a22−z2b9−b2c2+a122b1c2−b8−b10a22−z22,d6=c1ka123b4c2−b11+a11+1a122b12−b5c2−a11+12a12b13,d7=c1ka1222a11+a22−z2+2b12−b5c2−3a123b11−b4c2−a11+1a12b13a11+2a22−2z2+1,d8=c1a12a11+1b9−b2c2+a122b1c2−b8+a11+12−b10,d9=−c1−a12a11+a22−z2+1b9−b2c2+2a122b8−b1c2+2a11+1b10a22−z2,d10=c1a12a22−z2b9−b2c2+a122b1c2−b8−b10a22−z22,d11=c1ka122a11+2a22−2z2+1b12−b5c2−3a123b11−b4c2−a12b13a22−z22a11+a22−z2+2,d12=c1ka122a22−z2b12−b5c2+a123b4c2−b11−a12b13a22−z22,d13=c1ka12a11+1−b2c3−b9+a122b1c3+b8+a11+12b10,d14=c1−ka12a11+a22−z2+1b2c3+b9−2a122b1c3+b8−2a11+1b10a22−z2,d15=c1a12A1c3+A3+a11+1A2c3−A4,d16=c1a22−z2A2c3−A4+a12A1c3+A3,d17=c1k−a12a22−z2b2c3+b9+a122b1c3+b8+b10a22−z22,d18=c1ka123b4c3+b11−a11+1a122b5c3+b12+a11+12a12b13,d19=c1k−a1222a11+a22−z2+2b5c3+b12+3a123b4c3+b11+a11+1a12b13a11+2a22−2z2+1,d20=c1−a12a11+1b2c3+b9+a122b1c3+b8+a11+12b10,d21=−c1a12a11+a22−z2+1b2c3+b9−2a122b1c3+b8−2a11+1b10a22−z2,d22=c1−a12a22−z2b2c3+b9+a122b1c3+b8+b10a22−z22,d23=c1k−a122a11+2a22−2z2+1b5c3+b12+3a123b4c3+b11+a12b13a22−z22a11+a22−z2+2,d24=c1k−a122a22−z2b5c3+b12+a123b4c3+b11+a12b13a22−z22andX=a12u+v,Y=−1+a11u+z2−a22v.In order to implement the center manifold theorem, we assume thatMC be the center manifold of (30), evaluated at 0,0 in a small neighborhood of K=0. It can be approximated as follows.
(33)MC=K,u,v:v=GPDK,u,u<δ1,K<δ2,GPD0,0=0,DGPD0,0=0,where
(34)FPDK,GPDK,u=GPDK,z1u+f~PDK,u,GPDK,u−z2GPDK,u−g~PDK,u,GPDK,u=GPDK,−u+f~PDK,u,GPDK,u−z2GPDK,u−g~PDK,u,GPDK,u.We need to findGPDK,u such that FPDK,GPDK,u=0. Since we need the conditions in MC to be satisfied, we may assume
(35)GPDK,u=m1u2+m2Ku+m3K2+OK+u3.Thus, we can write(36)FPDK,GPDK,u=m11−z2−d13u2−d15+m2z2+1Ku+m31−z2K2+OK+u3.By comparing coefficients in the equationFPDK,GPDK,u=0, we get
(37)m1=d131−z2,m2=−d151+z2,m3=0,and
(38)GPDK,u=d131−z2u2−d151+z2Ku+OK+u3.The dynamics restricted toMC are given locally by the map
(39)FC:u↦z1u+f~PDK,u,GPDK,u=−u+d1u2+d3Ku−d4d151+z2K2u+d8−d4d131−z2−d2d151+z2Ku2+d6+d2d131−z2u3+OK+u4.Furthermore, define(40)Y1=∂2FC∂K∂u+12∂FC∂K∂2FC∂u20,0=d3≠0,Y2=16∂3FC∂u3+12∂2FC∂u220,0=d12+d6+d2d131−z2≠0.Theorem 5.
IfY1≠0 and Y2≠0, then system undergoes period-doubling bifurcation at the unique positive equilibrium point when parameter k varies in small neighborhood of k1 or k2. Moreover, the period-two orbits that bifurcate from positive equilibrium are stable if Y2>0 and unstable if Y2<0.Remark 6.
Since we have already found the center manifold for our system in (39), we can easily show the existence of transcritical or fold bifurcation, when one of the eigenvalue equals 1 and the other not on the unit circle.


## 4.2. Neimark-Sacker Bifurcation

Letα,β,h,m,c,ρ,k0∈ΩNS, i.e., k0=A1−A4/A1A4−A2A3. Then, variation of α,β,h,m,c,ρ,k in small neighborhood of ΩNS gives emergence to Neimark-Sacker bifurcation. Let K denote the perturbation parameter for the mapping (3), i.e., K1. The perturbed mapping is
(41)xy⟶x+k+Kx1−x−αym+x−hc+xy+k+Kρy1−βym+x.DefineX=x−x¯ and Y=y−y¯, to translate the fixed point to 0,0. Here, x¯,y¯ is the positive interior fixed point. Then, the above mapping is converted to
(42)XY⟶X+k+KX+x¯1−X+x¯−αY+y¯m+X+x¯−hc+X+x¯Y+k+KρY+y¯1−βY+y¯m+X+x¯.After Taylor series expansion, the above mapping can be written as(43)XY⟶c11c12c21c22XY+fNSX,YgNSX,Y,where c11=1+k+KA1, c12=−k+KA2, c21=k+KA3, and c22=1−k+KA4, and
(44)fNSX,Y=k+Kb1X2+k+Kb2XY+k+Kb3Y2+kb4X3+kb5X2Y+kb6XY2+kb7Y3+OX+Y4,gNSX,Y=k+Kb8X2+k+Kb9XY+k+Kb10Y2+kb11X3+kb12X2Y+kb13XY2+kb14Y3+OX+Y4,and bi,i∈1,2,⋯,14 are defined as above. Let the characteristic polynomial of matrix in (43) be
(45)pNSz=z2−PKz+QK,where
(46)PH=c11+c22,QH=c11c22−c21c12.Sinceα,β,h,m,c,ρ,k∈ΩNS, therefore the roots of (45) are complex conjugate z1, z2, and z1=z2=1, and it follows that z1,2=PH/2±1/24QH−PH2 and 1=z1=z2=QH. Note that PK2−4QK<0, which implies that PK2<4QK or PK∈−2,2. Also,
(47)dz1,2dKK=0=12ρ−x¯hc+x¯2+αy¯m+x¯2−1.Thus, ifρ≠x¯h/x¯+c2+αy¯/x¯+m2−1, then
(48)dz1,2dKK=0≠0.In order to ensure that the roots of the characteristic polynomial do not lie in the intersection of unit circle of coordinate axis whenK=0, we need to check that z1,2m≠1, for m=1,2,3,4 at K=0. This is equivalent to checking P0≠−2,0,1,2. Since α,β,h,m,c,ρ,k∈ΩNS, we already know that P0≠−2,0,2. Finally,
(49)P0≠1⇔k≠1A4−A1.Thus, the roots of (45) do not lie in the intersection of unit circle of coordinate axis when K=0. Now, to convert the system (43) to the normal form, when K=0, let R=P0/2 and S=1/24Q0−P02 and let
(50)TNS=01SR.Here,TNS is an invertible matrix. Consider the transformation
(51)XY=TNSuv.The system (43) can be written as
(52)uv⟶R−SSRuv+f~NSu,vg~NSu,v,where
(53)f~NSu,v=e1u2+e2uv+e3v2+e4u3+e5u2v+e6uv2+e7v3+Ou+v4,g~NSu,v=e8u2+e8uv+e10v2+e11u3+e12u2v+e13uv2+e14v3+Ou+v4,and
(54)e1=k+K−Rb2S+b1+Sb10S+b9+b8S,e2=b2−R+2b10S+b9Rk+KS,e4=k−Rb5S+b4+Sb13S+b12+b11S,e3=b10R2k+KS,e5=kR2b13S+b12−b5RS,e6=b13kR2S,e7=0,e8=b2S+b1k+K,e9=b2Rk+K,e10=0,e11=kb5S+b4,e12=b5kR,e13=0,e14=0,X=v,Y=Su+Rv.For the map (52) to undergo Neimark-Sacker bifurcation, the following quantity must not zero [22]:
(55)MNS=−Reϖ20ϖ111−2z1z221−z1−12ϖ112−ϖ022+Reϖ21z2K=0,where
(56)ϖ20=18∂2f¯NS∂u2−∂2f¯NS∂v2+2∂2g¯NS∂u∂v+ι∂2g¯NS∂u2−∂2g¯NS∂v2−2∂2f¯NS∂u∂v,ϖ11=14∂2f¯NS∂u2+∂2f¯NS∂v2+ι∂2g¯NS∂u2+∂2g¯NS∂v2,ϖ02=18∂2f¯NS∂u2−∂2f¯NS∂v2−2∂2g¯NS∂u∂v+ι∂2g¯NS∂u2−∂2g¯NS∂v2+2∂2f¯NS∂u∂v,ϖ21=116∂3f¯NS∂u3+∂3f¯NS∂u∂v2+∂3g¯NS∂u2∂v+∂3g¯NS∂v3+ι∂3g¯NS∂u3+∂3g¯NS∂u∂v2−∂3f¯NS∂u2∂v−∂3f¯NS∂v3.Using the result in [57] and the above analysis, we have the following result.Theorem 7.
Ifk≠1/A4−A1 and MNS defined in (55) is nonzero, then the model (5) undergoes Neimark-Sacker bifurcation at the equilibrium point x¯,y¯, provided that parameter k changes in the small neighborhood of k0=A1−A4/A1A4−A2A3 and α,β,h,m,c,ρ,k0∈ΩNS. Moreover, if MNS<0 (resp., MNS>0), then an attracting (resp., repelling) invariant closed curve bifurcate from the fixed point x¯,y¯ for k>k0 (resp., k<k0).


## 5. Numerical Simulations

In this section, we will present some examples which will show the presence of period-doubling and Neimark-Sacker bifurcations for the system (3), using specific numerical values of its parameters α,β,h,m,c,ρ and taking the step-size, k, as the bifurcation parameter. The illustration will be done using bifurcation diagrams, and it will be ratified by showing that Theorems 5 and 7 are satisfied.Example 1.
Selectα,β,h,m,c,ρ=0.4,1,0.1,0.1,0.05,2.6, k∈0.8,1.14, and the initial value x0,y0=0.8,0.4. The Jacobian matrix is
(57)1+0.179861k−0.311111k2.6k1−2.6k,and z1=1−2.26979k and z2=1−0.150344k are the corresponding eigenvalues. Both eigenvalues are less then 1, for any k>0. Also, z1=−1 and z2∈−1,1, if k=0.881137 and z1<−1 and z2=−1, if k=13.3028. Thus, with these parametric conditions, the fixed point is never unstable or saddle. The fixed point is stable for any k≠0.881137, or k=13.3028; otherwise, the fixed point undergoes period-doubling bifurcation, as k varies in the small neighborhood of k0, where let k0=0.881137. From the phase portrait in Figure 2(b), it is apparent that the fixed point 0.35,0.45 of map (3) is stable for k<0.881137 and the period-doubling bifurcation occurs around 0.881137. Figure 2(c) shows that another stable period-doubling bifurcation at around k=1.06. This can also be verified using the bifurcation diagrams in Figure 3.
After calculations, we get,Y1=0.330205≠0 and Y2=336.675≠0. These values further complement the dynamics of our map, observed above and proves the correctness of Theorem 5. Moreover, Y2>0, which shows that the period-two orbits that bifurcate from positive equilibrium are stable.Example 2.
Selectα,β,h,m,c,ρ=0.01,0.6,0.2,0.4,0.1,1.9, k∈1,1.42, and the initial value x0,y0=0.3,0.6. The Jacobian matrix is
(58)1−0.528241k−0.00651351k3.16667k1−1.9k,and z1=1−1.8848k and z2=1−0.543446k are the corresponding eigenvalues. Both eigenvalues are less then 1, for any k>0. Also, z1=−1 and z2∈−1,1, if k=1.06112 and z1<−1 and z2=−1, if k=3.68022. Thus, with these parametric conditions, the fixed point is never unstable or saddle. The fixed point is stable for any k≠1.06112, or k=3.68022; otherwise, the fixed point undergoes period-doubling bifurcation, as k varies in the small neighborhood of k0, where let k0=1.06112. From the phase portrait in Figure 4(b), it is apparent that the fixed point 0.747285,1.91214 of map (3) is stable for k<1.06112 and the period-doubling bifurcation occurs around 1.06112. Figure 4(c) shows that another stable period-doubling bifurcation at around k=1.31. This can also be verified using the bifurcation diagrams in Figure 5.Figure 2
Phase portraits for the system, with initial conditionsx0,y0=0.8,0.4, at parametric values α,β,h,m,c,ρ=0.4,1,0.1,0.1,0.05,2.6, around the bifurcation parameter k0=0.881137.
(a)
k=0.88(b)
k=0.89(c)
k=1.06(d)
k=1.13Figure 3
Period-Doubling bifurcation diagrams for positive interior points of the model (3), in the interval k∈0.8,1.14, with initial conditions x0,y0=0.8,0.4 and the parametric values α,β,h,m,c,ρ=0.4,1,0.1,0.1,0.05,2.6.
(a)(b)Figure 4
Phase portraits for the system, with initial conditionsx0,y0=0.3,0.6, at parametric values α,β,h,m,c,ρ=0.01,0.6,0.2,0.4,0.1,1.9, around the bifurcation parameter k0=1.06112.
(a)
k=1.06(b)
k=1.062(c)
k=1.31(d)
k=1.38Figure 5
Period-doubling bifurcation diagrams for positive interior points of the model (3), in the interval k∈1,1.42, with initial conditions x0,y0=0.3,0.6 and the parametric values α,β,h,m,c,ρ=0.01,0.6,0.2,0.4,0.1,1.9.
(a)(b)After calculations, we get,Y1=−0.23876≠0 and Y2=8.71736≠0. These values further complement the dynamics of our map, observed above and proved the correctness of Theorem 5. Moreover, Y2>0, which shows that the period-two orbits that bifurcate from positive equilibrium are stable.


### 5.1. Neimark-Sacker Bifurcation

Example 3.
Let us assume the following parameter values,α,β,h,m,c,ρ=0.05,2.3,0.001,0.05,0.5,1.275, k∈1.75,1.85, and the initial value x0,y0=0.955,0.3. The Jacobian matrix is
(59)1−0.956455k−0.0475671k0.554348k1−1.275k,and z1=−0.998393−0.056668ι and z2=−0.998393+0.056668ι are the corresponding eigenvalues, for k=k0=1.79111. Both eigenvalues are on the unit circle. Also, k0≠1/A4−A1=0.448138, k0∈0,4/A4−A1=0,1.79255, dz1/dKK=0=dz2/dKK=0=1.11573≠0 and MNS=−1297652.848483635≠0. Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the fixed point x¯,y¯=0.977584,0.446776. Moreover, since MNS<0, an attracting invariant closed curve bifurcates from the fixed point 0.977584,0.446776 for k>1.79111.From the phase portrait in Figure6(a), the attractor becomes a stable spiral point and as the parameter k increases, and the stable spiral enlarges in size. As k continues to increase, the spiral becomes larger before changing in to an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation. This spiral and edges are plotted in Figure 6(b) at k=1.795. At k=1.81, the attracting invariant closed curves become larger with some rough edges as shown in Figure 6(c). It is clear from the close values of the parameter k how quasi-periodic motions are produced around fixed point due to this type of bifurcation. Increasing k furthermore gives rise to continued invariant closed curves around the fixed point, but their edges begin to vanish as plotted in Figure 6(d).Example 11.
Let us assume that fork∈3.71,3.85 and the initial value x0,y0=0.45,1.1, we have the following parameter values, α,β,h,m,c,ρ=0.5,1.1,0.02,0.9,0.05,0.4. The Jacobian matrix is
(60)1−0.31283k−0.180788k0.363636k1−0.4k,and z1=−0.331058−0.943611ι and z2=−0.331058+0.943611ι are the corresponding eigenvalues, for k=k0=3.734575030140114. Both eigenvalues are on the unit circle. Also, k0≠1/A4−A1=1.40286, k0∈0,4/A4−A1=0,5.61144, dz1/dKK=0=dz2/dKK=0=0.356415≠0, and MNS=−6.56955≠0. Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point x¯,y¯=0.50972255689403,1.2815659608127545. Moreover, since MNS<0, an attracting invariant closed curve bifurcates from the fixed point for k>3.734575030140114.Figure 6
Phase portraits for the parametric valuesα,β,h,m,c,ρ,k0=0.05,2.3,0.001,0.05,0.5,1.275,1.79111 at x¯,y¯=0.977584,0.446776.
(a)
A stable spiral point atk=1.78(b)
An invariant closed curve atk=1.795(c)
Larger attracting invariant closed curves with some rough edges atk=1.81(d)
Larger attracting invariant closed curves with edges starting to vanish atk=1.84From the phase portrait in Figure7(a), we can see that the attractor becomes a stable spiral point and as the parameter k increases, the stable spiral enlarges in size (Figure 7(b)). As k continues to increase, the spiral becomes larger before changing into an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation. This spiral and edges are plotted in Figure 7(c) at k=3.74. At k=3.79, we see continued invariant closed curves around the fixed point, but their edges begin to vanish as plotted in Figure 7(d).Figure 7
Phase portraits for model (3) for the parametric values at x¯,y¯=0.50972255689403,1.2815659608127545.
(a)
A small stable spiral atk=3.66(b)
A large stable spiral atk=3.71(c)
An invariant closed curve with many edges atk=3.74(d)
Larger attracting invariant closed curves with edges starting to vanish atk=3.79


## 5.1. Neimark-Sacker Bifurcation

Example 3.
Let us assume the following parameter values,α,β,h,m,c,ρ=0.05,2.3,0.001,0.05,0.5,1.275, k∈1.75,1.85, and the initial value x0,y0=0.955,0.3. The Jacobian matrix is
(59)1−0.956455k−0.0475671k0.554348k1−1.275k,and z1=−0.998393−0.056668ι and z2=−0.998393+0.056668ι are the corresponding eigenvalues, for k=k0=1.79111. Both eigenvalues are on the unit circle. Also, k0≠1/A4−A1=0.448138, k0∈0,4/A4−A1=0,1.79255, dz1/dKK=0=dz2/dKK=0=1.11573≠0 and MNS=−1297652.848483635≠0. Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the fixed point x¯,y¯=0.977584,0.446776. Moreover, since MNS<0, an attracting invariant closed curve bifurcates from the fixed point 0.977584,0.446776 for k>1.79111.From the phase portrait in Figure6(a), the attractor becomes a stable spiral point and as the parameter k increases, and the stable spiral enlarges in size. As k continues to increase, the spiral becomes larger before changing in to an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation. This spiral and edges are plotted in Figure 6(b) at k=1.795. At k=1.81, the attracting invariant closed curves become larger with some rough edges as shown in Figure 6(c). It is clear from the close values of the parameter k how quasi-periodic motions are produced around fixed point due to this type of bifurcation. Increasing k furthermore gives rise to continued invariant closed curves around the fixed point, but their edges begin to vanish as plotted in Figure 6(d).Example 11.
Let us assume that fork∈3.71,3.85 and the initial value x0,y0=0.45,1.1, we have the following parameter values, α,β,h,m,c,ρ=0.5,1.1,0.02,0.9,0.05,0.4. The Jacobian matrix is
(60)1−0.31283k−0.180788k0.363636k1−0.4k,and z1=−0.331058−0.943611ι and z2=−0.331058+0.943611ι are the corresponding eigenvalues, for k=k0=3.734575030140114. Both eigenvalues are on the unit circle. Also, k0≠1/A4−A1=1.40286, k0∈0,4/A4−A1=0,5.61144, dz1/dKK=0=dz2/dKK=0=0.356415≠0, and MNS=−6.56955≠0. Thus, the fixed point exhibits Neimark-Sacker bifurcation, at the equilibrium point x¯,y¯=0.50972255689403,1.2815659608127545. Moreover, since MNS<0, an attracting invariant closed curve bifurcates from the fixed point for k>3.734575030140114.Figure 6
Phase portraits for the parametric valuesα,β,h,m,c,ρ,k0=0.05,2.3,0.001,0.05,0.5,1.275,1.79111 at x¯,y¯=0.977584,0.446776.
(a)
A stable spiral point atk=1.78(b)
An invariant closed curve atk=1.795(c)
Larger attracting invariant closed curves with some rough edges atk=1.81(d)
Larger attracting invariant closed curves with edges starting to vanish atk=1.84From the phase portrait in Figure7(a), we can see that the attractor becomes a stable spiral point and as the parameter k increases, the stable spiral enlarges in size (Figure 7(b)). As k continues to increase, the spiral becomes larger before changing into an attracting invariant closed curve with rough edges due to Neimark-Sacker bifurcation. This spiral and edges are plotted in Figure 7(c) at k=3.74. At k=3.79, we see continued invariant closed curves around the fixed point, but their edges begin to vanish as plotted in Figure 7(d).Figure 7
Phase portraits for model (3) for the parametric values at x¯,y¯=0.50972255689403,1.2815659608127545.
(a)
A small stable spiral atk=3.66(b)
A large stable spiral atk=3.71(c)
An invariant closed curve with many edges atk=3.74(d)
Larger attracting invariant closed curves with edges starting to vanish atk=3.79


## 6. Conclusion

In this paper, we studied the existence and stability of the positive interior fixed points of system (5) and flip bifurcation and Neimark-Sacker bifurcation, under certain conditions, by using central manifold theorem and bifurcation theory. Our main results are given in Theorem 3 and Theorems 5 and 7 and numerical simulations in Section 5. In the details of the result, when the integral step size k is chosen as a bifurcation parameter, the discrete-time modified Leslie-Gower system with Michaelis-Menten-type prey harvesting displays much richer dynamical behavior. It can be seen from Figures 3–9 that the system displays period-1, 2, 4, and 8 orbits. Moreover, Figures 3 and 5 show that the system exhibits period-doubling bifurcation, and the period-two orbits that bifurcate from positive equilibrium are stable, which shows the correctness of Theorem 5. Similarly, Figures 8 and 9 show that the system exhibits Neimark-Sacker bifurcation, with an attracting invariant closed curve bifurcates from the fixed point, proving the correctness of Theorem 7.Figure 8
Bifurcation diagrams for initial conditionsx0,y0=0.955,0.3 and k∈1.75,1.85, with parameters α,β,h,m,c,ρ=0.05,2.3,0.001,0.05,0.5,1.275, at x¯,y¯=0.977584,0.446776.
(a)(b)Figure 9
Bifurcation diagrams for initial conditionsx0,y0=0.45,1.1 and k∈3.71,3.85, at x¯,y¯=0.50972255689403,1.2815659608127545, with parameters α,β,h,m,c,ρ=0.5,1.1,0.02,0.9,0.05,0.4.
(a)(b)



---

*Source: 9575638-2022-03-12.xml*