# A Novel Discrete-Time Leslie–Gower Model with the Impact of Allee Effect in Predator Population

Authors: S. Vinoth; R. Sivasamy; K. Sathiyanathan; B. Unyong; R. Vadivel; Nallappan Gunasekaran
Journal: Complexity
 (2022)

---

## Abstract

The discrete-time system has more complex and chaotic dynamical behaviors as compared to the continuous-time system. This paper extends a discrete Leslie–Gower predator-prey system with the Allee effect in the predator’s population, whose dynamics are analyzed and explored. We have determined the equilibrium points and studied their local stability properties. We find that the system undergoes flip bifurcation and Neimark–Sacker bifurcation around the interior equilibrium point by choosing the Allee parameter as a bifurcation parameter. We discuss the stability and direction of both bifurcations with the help of the normal form theory and center manifold theorem. The flip bifurcation and Neimark–Sacker bifurcation are the most common routes to the chaotic orbit in the discrete system. Moreover, we utilize state feedback, pole placement, and hybrid control methods to control the chaos in the system. The work is complete with the numerical simulations to confirm the analytical findings.


---

## Body

## 1. Introduction

The study of interactions between prey and predator has particular interest for many mathematicians and ecologists. Many researchers investigated the dynamic behavior of the predator-prey system and contributed to the development of continuous-time systems. Another possible way to understand the complex problem of two interacting species is by discrete systems [1]. Moreover, it has been observed that the discrete-time systems are more appropriate than the continuous-time systems for the populations with nonoverlapping generations [2]. The discrete-time models are more appropriate and provide efficient results than continuous models for small-size populations [3]. For example, the discrete model is better for studying insect populations since there is only one generation per year and an annual plant population because there are no overlapping generations annually. Jing and Yang [4], Liu and Xiao [5], and Elabbasy et al. [6] showed that the discrete-time system showed richer and more complex dynamics compared to the continuous-time systems. Many researchers formulated and studied the discrete-time predator-prey system by implementing the forward Euler scheme [4, 5, 7], nonstandard finite difference scheme [8, 9], and piecewise constant arguments [10–12]. For instance, the authors in [13] obtained a discrete-time predator-prey system with a crowding effect and predator partially dependent on prey by applying the forward Euler’s scheme from the continuous-time system. They also showed the existence of a cascade of period-doubling bifurcation and Hopf bifurcation in the considered system. On the other hand, Abbasi et al. [12] used the method of piecewise constant arguments and obtained the discrete system. Also, the authors investigated stability, bifurcations, and chaos control analysis. In [14], the authors investigated the discrete-time predator-prey system with hunting cooperation through numerical simulation and observed that the system undergoes both types of bifurcations. In [15], the authors discussed the chaotic dynamics of a discrete prey-predator system with Holling-II type functional response. The other fruitful results on discrete-time models can be found in [16, 17].The Allee effect is a biological phenomenon named after Allee [18], who was the first to write thoroughly about it. It denotes the presence of a positive relationship between population size and per capita growth rate. This effect usually saturates or vanishes as the population expands. There has been a great deal of literature on deriving continuous-time prey-predator models with the Allee effect [19–22]. The impact of the Allee effect in the discrete-time prey-predator models has also been studied in the literature [23–26]. For example, Celik and Duman [23] studied the impact of the Allee effect in the discrete-time predator-prey model with linear interaction. They showed that the Allee effect on prey population changes the unstable equilibrium into a stable state. In [25], the authors studied the discrete-time predator-prey model with strong and weak Allee effects. Also, they showed that the chaotic orbits appear for larger step size values through period-doubled orbits and invariant circle orbits. AlSharawi et al. [26] considered the Allee effect in the discrete-time prey-predator model with a nonmonotonic functional response. They performed the stability analysis and provided conditions for the occurrence of flip and Neimark–Sacker bifurcations by taking the Allee parameter as a bifurcation parameter. Also, they showed that the system has two types of bistability behaviors.Several research works have been done for bifurcation and chaos control in nonlinear systems. It refers to the role of constructing a controller to alter the chaotic and bifurcating properties of a given nonlinear system to achieve some desired dynamical behaviors [27–30]. One can shift the chaotic attractor to any one of a large number of possible attracting periodic motions, as in [31]. Recently, the studies on bifurcation and chaos control in the discrete-time predator-prey models attained much more interest among researchers. For example, Din [32] showed that the model undergoes flip bifurcation (FB) and Neimark–Sacker bifurcation (NSB) for larger values of the growth rate of the prey population. Furthermore, the author implemented three different types of control strategies to control the chaos. Moreover, for some interesting results related to bifurcation and chaos control in the predator-prey models, we refer the readers to [33–35].Motivated by the above-mentioned works, we consider the modified Leslie–Gower prey-predator system introduced by Alaoui and Okiye in [36] by assuming that the growth rate of the predator is affected by the presence of the Allee effect. Also, the prey and predator populations have nonoverlapping generations. The discrete-time model is obtained by the method of piecewise constant arguments for the differential equation [10–12] from its corresponding continuous-time system. We attempt to study the discrete-time system and observe some rich dynamics that the continuous-time system does not have. Feng and Kang [20] investigated the [36] continuous-time model with an Allee effect in both prey and predator. To the best of our knowledge, there has been less research done in discrete predator-prey systems with the Allee effect in the predator population. The main purpose of this paper is to show the rich dynamics of the discrete system in terms of bifurcations (FB and NSB) and chaos by taking the Allee parameter as a bifurcation parameter. Flip and Neimark–Sacker bifurcations are the common routes to chaos in discrete systems. Additionally, if the discrete system is chaotic under certain parametric conditions, we can use various control methods to stabilize the chaotic orbits near the unstable equilibrium point. It is worth mentioning here that the study on stability, bifurcation, and chaos control analysis for the prey-predator model with Allee effect in predator population is different from the whole of the existing works.This paper is organized as follows: in Section2, we describe the formation of a discrete-time system for populations with nonoverlapping generations from the continuous counterpart. In Section 3, we study the existence and local stability of the equilibrium points for the discrete-time system. In Section 4, we illustrate the existence of flip and Neimark–Sacker bifurcation by taking the Allee parameter as a bifurcation parameter. Section 5 is related to the implementation of state feedback, pole-placement, and hybrid control methods to delay the chaos in the system. Lastly, to ensure our analytical results, the various numerical simulations are performed in Section 6.


## 2. Model Formulation

Firstly, we describe the continuous-time predator-prey system of modified Leslie–Gower with Holling type II functional response. It was introduced by Alaoui and Okiye in [36], which is of the following form:(1)dxtdt=xtr1−axt−gytxt+b,dytdt=ytr2−hytxt+c.x0≥0 and y0≥0, where xt and yt represent the population sizes of prey and predator at time t. r1, a, g, b, r2, h, and c are all positive constants. r1 and r2 are the growth rates of x and y. a measures the strength of competition among the individuals of species x. g is the maximum value a per capita reduction rate of x can attain. b (respectively, c) measures the extent to which the environment provides protection to prey x (respectively, to predator y), and h has a similar meaning to g. The assumption that the amount to which the environment provides protection to both the predator and the prey is the same (i.e., b=c) has been considered in [37, 38]. Singh et al. [39] extended system (1) with the death rate of predators, obtained a discrete system by Euler’s method, and also considered step size as a bifurcation parameter. This fact violates the numerical method of discretization. To overcome this, Din in [32] applied the method of piecewise constant arguments to obtain the discrete system. Recently, the discrete version of system (1) was derived using Euler’s method with step size 1 by the authors in [40] and showed bifurcation behavior with other model parameters.Now, we consider that system (1) is subject to the Allee effect in predator population, and it is given by(2)dxtdt=xtr1−axt−gxtytxt+b,dytdt=ytr2ytyt+m−hytxt+c,where yt/m+yt is the Allee effect term and m>0 represents the severity of Allee effect in the predator population. The authors in [41] studied system (2) with ratio-dependent functional response and the fear in prey population.Next, based on the appropriate modifications of overlapping generations, one can get the difference equation for modeling population with nonoverlapping generations. We aim to study the populations that have nonoverlapping generations for system (2). Moreover, it is necessary to obtain the discrete-time system from its continuous counterpart. In this way, the method of piecewise constant argument has been useful. The corresponding discrete-time system (2) is obtained by the method of piecewise constant arguments for the differential equations [10–12, 14], assuming that the populations have no overlap between successive generations and the population growth occurs in discrete steps t∈n,n+1, n=0,1,2,⋯. Let us consider that the variables and constants in (2) change in the regular time intervals and obtain the following modified system:(3)1xtdxtdt=r1−axt−gytxt+b,1ytdytdt=r2ytyt+m−hytxt+c,t≠0,1,2,…,where t is the integer part of t, t∈0,+∞. On any interval of the form n,n+1, n=0,1,2,…, we can integrate (3) and obtain the following:(4)lnxtxn=r1−axn−gynxn+bt−n,lnytyn=r2ytyt+m−hynxn+ct−n,n=0,1,2,….By simplification and lettingt⟶n+1, the corresponding discrete-time system for (2) is obtained by the method of piecewise constant arguments for differential equations can be written as follows:(5)xn+1=xnexpr1−axn−gynxn+b,yn+1=ynexpr2ynyn+m−hynxn+c,where xn+1 and yn+1 denote the populations in generation n+1, and they are related to the sizes xn and yn of the populations in the preceding generation n. Note that in the absence of a predator, system (5) reduces to the one-dimensional system, similar to the Ricker model [16]. Thus, this current article aims to analyze the local stability, bifurcation, and chaos control analyses for the discrete system (5) that models the interaction between populations that have nonoverlapping generations.


## 3. The Existence and Local Stability of the Equilibia

### 3.1. The Equilibria

To find the equilibrium points of system (5), we use direct substitution method to solve the following equations:(6)r1−ax−gyx+b=0,r2y+m−hx+c=0.Form the above equations, we have the following points of equilibria:The originE0=0,0.The predator-free equilibriumE1=x¯,0=r1/a,0.The prey-free equilibriumE2=0,y¯=0,r2c/h−m exists if r2c/h>m.The interior equilibrium point isE∗x∗,y∗, where y∗=r2x∗/h+r2c−mh/h and x∗ is the positive root of the following equation:(7)ahx∗2+abh+gr2−r1hx∗+gr2c−mh−r1hb=0,Ifr2c−mh/h<r1b/g holds, then (6) has at least one positive real root x∗.To guarantee the existence of interior equilibrium pointE∗x∗,y∗, we need the following lemma.Lemma 1.
Let us assume thatr2c−mh/h<r1b/g and r2x∗+c>mh always hold. Then, E∗x∗,y∗ is the unique positive interior equilibrium point of system (5).


### 3.2. Local Stability Analysis

To analyze the local stability properties of the equilibria, we need the Jacobian matrix at an arbitrary equilibriumEx,y, which is given as follows:(8)J=1−ax+gxyx+b2A1−gxx+bA1hy2c+x2A21−hyc+x+r2myy+m2A2,where A1=expr1−ax−gy/x+b and A2=expr2y/y+m−hy/x+c.Then, we have the subsequent lemmas for the local stability analysis of the equilibria.Lemma 2.
For system (5), we have,(1)
E0, E1 are nonhyperbolic points.(2)
(a)
E20,y¯ is a sink if expr1−gy¯/b<1 and 1−hy¯/c+r2my¯/y¯+m2A3<1.(b)
E20,y¯ is a source if expr1−gy¯/b>1 and 1−hy¯/c+r2my¯/y¯+m2A3<1.(c)
E20,y¯ is a saddle if expr1−gy¯/b<1, 1−hy¯/c+r2my¯/y¯+m2A3<1, expr1−gy¯/b<1, and 1−hy¯/c+r2my¯/y¯+m2A3>1.(d)
E20,y¯ is a nonhyperbolic if expr1−gy¯/b=11−hy¯/c+r2my¯/y¯+m2A3=1.Proof.
The Jacobian matrix atE00,0 is as follows:(9)JE0=expr1001,with the eigenvalues λ1,2=expr1,1.
The Jacobian matrix atE1x¯,0 is as follows:(10)JE1=1−r1−gr1ab+r101,with the eigenvalues λ1,2=exp1−r1,1.
The Jacobian matrix atE20,y¯ is as follows:(11)JE2=expr1−gy¯b0hy¯2c2A31−hy¯c+r2my¯y¯+m2A3,with A3=expr2y¯/y¯+m−hy¯/c, and the eigenvalues are λ1=expr1−gy¯/b and λ2=1−hy¯/c+r2my¯/y¯+m2A3.Lemma 3.
The interior equilibrium point(a)
E∗x∗,y∗ is a sink if B2<1 and B1<B2+1.(b)
E∗x∗,y∗ is a source if B2>1 and B1<B2+1 or B1>B2+1.(c)
E∗x∗,y∗ is a saddle if 0<B1+B2+1<2B1.(d)
E∗x∗,y∗ is nonhyperbolic if B1=B2+1, or B2=1, and B1≤2.Proof.
The Jacobian matrix atE∗x∗,y∗ is as follows:(12)JE∗=1−ax∗+gx∗y∗x∗+b2−gx∗x∗+bhy∗2x∗+c21−r2y∗2y∗+m2.
Then, the characteristic polynomial ofJE∗ is as follows:(13)Fλ≔λ2−B1λ+B2=0,where(14)B1=2−ax∗+gx∗y∗x∗+b2−r2y∗2y∗+m2,B2=1−ax∗+gx∗y∗x∗+b21−r2y∗2y∗+m2+ghx∗y∗2x∗+bx∗+c2.


### 3.3. Bistability

Since model (5) is with the Allee effect in predator population, it is possible for the occurrence of bistability criteria, which is for two positive equilibrium points to exist and for both to be stable. From Lemma 2,E0 and E1 are nonhyperbolic points. The only possibility is between E2 and E∗. Also, from Lemma 1,E2 exists if r2c/h>m, and E∗ exists if r2c−mh/h<r1b/g and r2x∗+c>mh holds. If r2c−mh/h<r1b/g, then E2 is unstable form Lemma 2 (a), since one of the eigenvalues of JE2 is expr1−y¯/b>1. Then, we can say that if E∗ exists, then E2 is unstable. Hence, there is less possibility for the occurrence of bistability for model (5).


## 3.1. The Equilibria

To find the equilibrium points of system (5), we use direct substitution method to solve the following equations:(6)r1−ax−gyx+b=0,r2y+m−hx+c=0.Form the above equations, we have the following points of equilibria:The originE0=0,0.The predator-free equilibriumE1=x¯,0=r1/a,0.The prey-free equilibriumE2=0,y¯=0,r2c/h−m exists if r2c/h>m.The interior equilibrium point isE∗x∗,y∗, where y∗=r2x∗/h+r2c−mh/h and x∗ is the positive root of the following equation:(7)ahx∗2+abh+gr2−r1hx∗+gr2c−mh−r1hb=0,Ifr2c−mh/h<r1b/g holds, then (6) has at least one positive real root x∗.To guarantee the existence of interior equilibrium pointE∗x∗,y∗, we need the following lemma.Lemma 1.
Let us assume thatr2c−mh/h<r1b/g and r2x∗+c>mh always hold. Then, E∗x∗,y∗ is the unique positive interior equilibrium point of system (5).


## 3.2. Local Stability Analysis

To analyze the local stability properties of the equilibria, we need the Jacobian matrix at an arbitrary equilibriumEx,y, which is given as follows:(8)J=1−ax+gxyx+b2A1−gxx+bA1hy2c+x2A21−hyc+x+r2myy+m2A2,where A1=expr1−ax−gy/x+b and A2=expr2y/y+m−hy/x+c.Then, we have the subsequent lemmas for the local stability analysis of the equilibria.Lemma 2.
For system (5), we have,(1)
E0, E1 are nonhyperbolic points.(2)
(a)
E20,y¯ is a sink if expr1−gy¯/b<1 and 1−hy¯/c+r2my¯/y¯+m2A3<1.(b)
E20,y¯ is a source if expr1−gy¯/b>1 and 1−hy¯/c+r2my¯/y¯+m2A3<1.(c)
E20,y¯ is a saddle if expr1−gy¯/b<1, 1−hy¯/c+r2my¯/y¯+m2A3<1, expr1−gy¯/b<1, and 1−hy¯/c+r2my¯/y¯+m2A3>1.(d)
E20,y¯ is a nonhyperbolic if expr1−gy¯/b=11−hy¯/c+r2my¯/y¯+m2A3=1.Proof.
The Jacobian matrix atE00,0 is as follows:(9)JE0=expr1001,with the eigenvalues λ1,2=expr1,1.
The Jacobian matrix atE1x¯,0 is as follows:(10)JE1=1−r1−gr1ab+r101,with the eigenvalues λ1,2=exp1−r1,1.
The Jacobian matrix atE20,y¯ is as follows:(11)JE2=expr1−gy¯b0hy¯2c2A31−hy¯c+r2my¯y¯+m2A3,with A3=expr2y¯/y¯+m−hy¯/c, and the eigenvalues are λ1=expr1−gy¯/b and λ2=1−hy¯/c+r2my¯/y¯+m2A3.Lemma 3.
The interior equilibrium point(a)
E∗x∗,y∗ is a sink if B2<1 and B1<B2+1.(b)
E∗x∗,y∗ is a source if B2>1 and B1<B2+1 or B1>B2+1.(c)
E∗x∗,y∗ is a saddle if 0<B1+B2+1<2B1.(d)
E∗x∗,y∗ is nonhyperbolic if B1=B2+1, or B2=1, and B1≤2.Proof.
The Jacobian matrix atE∗x∗,y∗ is as follows:(12)JE∗=1−ax∗+gx∗y∗x∗+b2−gx∗x∗+bhy∗2x∗+c21−r2y∗2y∗+m2.
Then, the characteristic polynomial ofJE∗ is as follows:(13)Fλ≔λ2−B1λ+B2=0,where(14)B1=2−ax∗+gx∗y∗x∗+b2−r2y∗2y∗+m2,B2=1−ax∗+gx∗y∗x∗+b21−r2y∗2y∗+m2+ghx∗y∗2x∗+bx∗+c2.


## 3.3. Bistability

Since model (5) is with the Allee effect in predator population, it is possible for the occurrence of bistability criteria, which is for two positive equilibrium points to exist and for both to be stable. From Lemma 2,E0 and E1 are nonhyperbolic points. The only possibility is between E2 and E∗. Also, from Lemma 1,E2 exists if r2c/h>m, and E∗ exists if r2c−mh/h<r1b/g and r2x∗+c>mh holds. If r2c−mh/h<r1b/g, then E2 is unstable form Lemma 2 (a), since one of the eigenvalues of JE2 is expr1−y¯/b>1. Then, we can say that if E∗ exists, then E2 is unstable. Hence, there is less possibility for the occurrence of bistability for model (5).


## 4. Bifurcation Analysis

We have studied the stability properties ofE∗x∗,y∗ previously. In the subsequent section, we choose the Allee parameter m as a bifurcation parameter to analyze the bifurcation behaviour of system (5). We derive the conditions to obtain flip bifurcation (FB) and Neimark–Sacker bifurcation (NSB) at E∗x∗,y∗. Moreover, we utilize the center manifold theorem and the normal form theory [42, 43] to discuss the direction and stability property of the flip bifurcation and Neimark–Sacker bifurcation at E∗x∗,y∗.


### 4.1. Flip Bifurcation

System (5) undergoes flip bifurcation if one of the eigenvalues ofJE∗ must be −1 and the other should not be 1 or −1. Using this, we assume one of the eigenvalues of JE∗ as −1. Then, from (12), we have the following:(15)F−1:=4−2ax∗−2r2y∗2y∗+m2+2gx∗y∗x∗+b2+ar2x∗y∗2y∗+m2+ghx∗y∗2x∗+bx∗+c2−gr2x∗y∗3x∗+b2y∗+m2=0.Note that the left-hand side of expression (13) is in terms of parameter m in denominator with power. Hence, it is difficult to find the explicit value of m for which system (5) undergoes flip bifurcation. Hence, we denote m=mh as the critical parameter value that satisfies (13). Also, E∗x∗,y∗ is derived at the critical value m=mf.Next, ifa2>0, a3<0, and hy∗+m/r2>c holds, we define that the neighborhood Θf for system (5) undergoes flip bifurcation near the equilibrium point E∗ at some critical value m=mf as follows:(16)Θf=r1,r2,a,b,c,g,h,m:m=mf,r1,r2,a,b,c,g,h>0.Hence, we assume that system (5) undergoes flip bifurcation when m=mf changes in the neighborhood Θf.Now, we investigate the direction and stability of the possible occurrence of flip bifurcation for system (5) atE∗x∗,y∗ with the steps followed in [32, 42]. Since, r1,r2,a,b,c,g,h,m∈Θf on giving a perturbation m1≪1 of critical value m=mf for which system (5) undergoes flip bifurcation, then the perturbation system is given by the following:(17)xn+1=xnexpr1−axn−gynxn+b,yn+1=ynexpr2ynyn+m1+mf−hynxn+c.Next, we transform the equilibrium pointE∗x∗,y∗ into the origin by letting pn=xn−x∗ and qn=yn−y∗. Form (17), we obtain the following:(18)pn+1=pn+x∗expr1−apn+x∗−gqn+y∗pn+x∗+b−x∗,qn+1=qn+y∗expr2qn+y∗qn+y∗+m1+mf−hqn+y∗pn+x∗+c−y∗.By Taylor’s series expansion, (16) becomes,(19)pn+1=α1pn+α2qn+α3m1+α4pn2+α5qn2+α6m12+α7pnqn+α8pnm1+α9qnm1+opn+qn+m12,qn+1=β1pn+β2qn+β3m1+β4pn2+β5qn2+β6m12+β7pnqn+β8pnm1+β9qnm1+opn+qn+m12,where(20)α1=1+x∗M1,α2=−gx∗x∗+b′,α3=0,α4=−gx∗y∗x∗+b+M1x∗2M12,α5=g2x∗2x∗+b2,α6=0,α7−gbx∗+b2−gM1x∗x∗+b,α8=0,α9=0,β1=hy∗2x∗+c2,β2=1+y∗M2,β3=−r2y∗2y∗+mf2,β4=−hy∗2x∗+c3+h2y∗32x∗+c4,β5=−r2mfy∗y∗+mf3+M2+y∗M222,β6=r22y∗32y∗+mf4+r2y∗2y∗+mf3,β7=2hy∗x∗+c2+hy∗2M2x∗+c2,β8=−hr2y∗3x∗+c2y∗+b2,β9=−2mfr2y∗y∗+mf3−r2M2y∗2y∗+mf2,M1=−α+gy∗x∗+b2,M2=−r2y∗y∗+mf2.Consider(21)JE∗=1−ax∗+gx∗y∗x∗+b2−gx∗x∗+bhy∗2x∗+c21−r2y∗2y∗+mf2.Let us assume that the matrixJE∗ has eigenvalues.(22)λ1=−1,λ2=1−ax∗+gx∗y∗x∗+b2+ghx∗y∗2y∗+mf2x∗+bx∗+c22y∗+mf2−r2y∗,satisfiesλ1=1 and λ2≠1.Next, we construct the nonsingular matrixL as follows:(23)L=α2α2−1−α1λ2−α1.Using translationpnqn=LPnQn. Then, (17) can be written as(24)Pn+1=−Pn+G1pn,qn,m1+opn+qn+m12,Qn+1=λ2Qn+G2pn,qn,m1+opn+qn+m12,where(25)G1pn,qn,m1=ϕ1m1+ϕ2pn2+ϕ3qn2+ϕ4m12+ϕ5pnqn+ϕ6pnm1+ϕ7qnm1,G2pn,qn,m1=θ1m1+θ2pn2+θ3qn2+θ4m12+θ5pnqn+θ6pnm1+θ7qnm1,and(26)ϕ1=λ2−α1α3−α2β3α21+λ2,ϕ2=λ2−α1α4−α2β4α21+λ2,ϕ3=λ2−α1α5−α2β5α21+λ2,ϕ4=λ2−α1α6−α2β6α21+λ21+λ2,ϕ5=λ2−α1α7−α2β7α21+λ2,ϕ6=λ2−α1α8−α2β8α21+λ2,ϕ7=λ2−α1α9−α2β9α21+λ2,θ1=1+α1α3+α2β3α21+λ2,θ2=1+α1α4+α2β4α21+λ2,θ3=1+α1α5+α2β5α21+λ2,θ4=1+α1α6+α2β6α21+λ2,θ5=1+α1α7+α2β7α21+λ2,θ6=1+α1α8+α2β8α21+λ2,θ7=1+α1α9+α2β9α21+λ2.Now, by the center manifold theorem, we obtain the center manifoldGc0,0 of system (19) at 0,0 in a small neighborhood of m1=0 as follows:(27)Gc0,0=Pn,Qn:Qn=hPn,m1,=Pn,Qn:Qn=e1m1+e2Pn2+e3m1Pn+e4m12+oPn+m12.On applying map (19) on both sides of Qn=hPn,m1, we have the following:(28)λ2Qn+G2pn,qn,m1=e1m1+e2−Pn+G1pn,qn,m12+e3m1−Pn+G1pn,qn,m1+e4m12+oPn+m12,where(29)pn=α2Pn+Qn=α2Pn+hPn,m1,qn=−1−α1Pn+λ2−α1Qn=−1−α1Pn+λ2−α1hPn,m1,e1=θ11−λ2,e2=11−λ2θ2α22+θ31+α12−θ5α2−1−α1,e3=11+λ2−2e2ϕ1−θ6α2+θ7−1−α1−2e1θ2α22+2e1θ3−1−α1λ2−α1−e1θ5α2β2−α1,e4=11−λ2e12θ2α22+e12θ3λ2−α12+e12θ5α2λ2−α1+e1θ6α2+e1θ7λ2−α1−e2ϕ12+θ4−ϕ1e3.Therefore, on the center manifoldMc at origin, we have the following:(30)pn2=β22Pn2+2PnQn+Qn2,pnqn=−β21+β1Pn2+β2δ2−β1PnQn+β2λ2−β1Qn2,qn2=1+β12Pn2−21+β1λ2−β1PnQn+λ2−β12Qn2,where(31)PnQn=e1m1Pn+e2Pn3+e3m1Pn2+e4m12Pn+oPn+m13,Qn2=e12m12+2e1e2m1Pn2+2e1e3m12Pn+2e1e4m13+oPn+m13.Moreover, the mapG∗ restricted to CM Gc0,0 is as follows:(32)G∗Pn=−Pn+G1pn,qn,m1=−Pn+g1m1+g2Pn2+g3Pnm1+g4m12+g5Pn2m1+g6Pnm12+g7Pn3+g8m13+oPn+m13,where(33)g1=ϕ1,g2=ϕ2β22+ϕ31+β12−ϕ5β21+β1,g3=2e1ϕ2β22−2e1ϕ31+β1λ2−β1+e1ϕ5β2λ2−β1+ϕ6β2−ϕ71+β1,g4=e12ϕ2β22+e12ϕ3λ2−β12+ϕ4+e12ϕ5β2λ2−β1+e1ϕ6β2+e1ϕ7λ2−β1,g5=2e3ϕ2β22+2e1e2ϕ2β22−2e3ϕ31+β1λ2−β1+2e1e2ϕ3λ2−β12+e3ϕ5β2λ2−β1+2e1e2ϕ5β2λ2−β1+e2ϕ6β2+e2ϕ7λ2−β1,g6=2e4ϕ2β22+2e1e3ϕ2β22−2e4ϕ31+β1λ2−β1+2e1e3ϕ3λ2−β12+e4ϕ5β2λ2−β1+2e1e3ϕ5β2λ2−β1+e3ϕ6β2+e3ϕ7λ2−β1,g7=2e2ϕ2β22−2e2ϕ31+β1λ2−β1+e2ϕ5β2λ2−β1,g8=2e1e4ϕ2β22+2e1e4ϕ3λ2−β12+2e1e4ϕ5β2λ2−β1+e4ϕ6β2+e4ϕ7λ2−β1.From [43], we define ζ1 and ζ2 as follows:(34)ζ1=GPnm1∗+12Gm1∗GPnPn∗|Pn,m1=0,0=g3+g1g2,(35)ζ2=16GPnPnPn∗+12GPnPn∗2|Pn,m1=0,0=g7+g22.Therefore, we have got the following findings of flip bifurcation from the aforementioned study.Theorem 1.
Ifζ1≠0 and ζ2≠0, then system (5) exhibits a flip bifurcation at E∗x∗,y∗ when m1 changes in a small neighborhood of the origin. Moreover, if ζ2>0 (or ζ2<0), then the existing period-two orbits from E∗x∗,y∗ are stable (or unstable).


### 4.2. Neimark–Sacker Bifurcation

In this subsection, we find the condition for the existence of Neimark–Sacker (NS) bifurcation and also its direction and stability properties for system (5) nearE∗x∗,y∗ using the method followed in [32, 42]. For the existence of NS bifurcation of system (5), the complex conjugate eigenvalues of the characteristic equation at E∗x∗,y∗ should have absolute value one. For this, it is necessary to satisfy the following:(36)B1x∗,y∗2−4B2x∗,y∗<0andB2x∗,y∗=1,which provides(37)Am:=2−ax∗+gx∗y∗x∗+b2−r2y∗2y∗+m22−41−ax∗+gx∗y∗x∗+b2−r21−ax∗y∗2y∗+m2−gr2x∗y∗3x∗+b2y∗+m2+ghx∗y∗2x∗+bx∗+c2<0,Bm:=−ax∗+gx∗y∗x∗+b2−r21−ax∗y∗2y∗+m2−gr2x∗y∗3x∗+b2y∗+m2+ghx∗y∗2x∗+b2x∗+c2=0.Note that the left-hand side of (37) is in terms of m since m is taken as a bifurcation parameter. It is difficult to find the explicit expression for the critical value m to satisfy (24) for the occurrence of NS bifurcation of system (5). Therefore, we assume that the critical value is m=mh. Also, the interior equilibrium E∗x∗,y∗ is calculated at the critical value m=mh.Next, ifa2>0, a3<0, and hy∗+m/r2>c holds, we define the neighborhood as follows:(38)Θh=r1,r2,a,b,c,g,h,m:m=mh,r1,r2,a,b,c,g,h>0.Then, the equilibrium pointE∗x∗,y∗ can arise NS bifurcation at m=mh when it changes in the neighborhood of Θh. Now, we analyze the properties of possible NS bifurcation at E∗x∗,y∗ for system (5) if (24) holds for some mh. Given a perturbation m2≪1 of critical value mh, then the perturbation system is described as follows:(39)xn+1=xnexpr1−axn−gynxn+b,yn+1=ynexpr2ynyn+m2+mh−hynxn+c.Let us use the transformpn=xn−x∗, qn=yn−y∗ and shift E∗x∗,y∗ to 0,0. Then, system (28) takes the following form:(40)pn+1=pn+x∗expr1−apn+x∗−gqn+y∗E∗x∗,y∗+b−x∗,qn+1=qn+y∗expr2qn+y∗qn+y∗+m2+mh−hqn+y∗E∗x∗,y∗+c−y∗.Therefore, by Taylor expansion of (40),(41)pn+1=ρ1pn+ρ2qn+ρ3pn2+ρ4pnqn+ρ5qn2+ρ6pn3+ρ7pn2qn+ρ8pnvn2+ρ9vn3+opn+qn3,qn+1=η1pn+η2qn+η3pn2+η4pnqn+η5qn2+η6pn3+η7pn2qn+η8pnvn2+η9vn3+opn+qn3,where(42)ρ1=1+x∗N1,ρ2=−gx∗x∗+b,ρ3=−gx∗y∗x∗+b3+N1+x∗N122,ρ4=−gbx∗+b2−gN1x∗x∗+b,ρ5=g2x∗2x∗+b2,ρ6=−bgy∗x∗+b4−gN1x∗y∗x∗+b3+N122+x∗N136,ρ7=g2x∗y∗x∗+b4+gbx∗+b3−gbN1x∗+b2−gN12x∗2x∗+b,ρ8=−g2x∗x∗+b3+g21+x∗N12x∗+b2,ρ9=−g3x∗6x∗+b3,η1=hy∗2x∗+c2,η2=1+y∗N2,η3=−hy∗2x∗+c3+h2y∗32x∗+c4,η4=2hy∗x∗+c2+hy∗2N2x∗+c2,η5=−r2mhy∗y∗+mh3+N2+y∗N222,η6=hy∗2x∗+c4−h2y∗3x∗+c5+h3y∗46x∗+c6,η7=−hy∗2+y∗N2x∗+c3+h2y∗23+y∗N22x∗+c4,η8=h1+y∗2N222x∗+c2+2hy∗N2x∗+c2−r2mhy∗+mh2,η9=r2mhy∗y∗+mh4−r2mh1+N2y∗y∗+mh3+N222+N23y∗6,N1=−a+gy∗x∗+b2,N2=−r2y∗y∗+mh2.The characteristic polynomial equation associated with the linearized system (30) at the origin can be given as follows:(43)λ2+q1m2λ+q2m2=0,where(44)q1m2=−1+ax∗−gx∗y∗x∗+b2Ω1−1−r2y∗2y∗+mh+m22Ω2,q2m2=1−ax∗+gx∗y∗x∗+b21−r2y∗2y∗+mh+m22+ghx∗y∗2x∗+bx∗+c2Ω1Ω2,WithΩ1=expr1−ax∗−gy∗/x∗+b and Ω2=expr2y∗/y∗+mh+m2−hy∗/x∗+c. Now, the roots of (43) are the pair of complex conjugates.(45)λ1,2=12−q1m2±i4q2m2−q1m22.Sincer1,r2,a,b,c,g,h,m∈Θh, we have,λ1,2=q2m2and(46)dλ1,2dm2=12q201−ax∗+gx∗y∗x∗+b22r2y∗2y∗+mh3−r2y∗y∗+mh21−ax∗+gx∗y∗x∗+b2×1−r2y∗2y∗+mh2+ghx∗y∗2x∗+bx∗+c2≠0.Furthermore, it is requiredλ1k,λ2k≠1 for k=1,2,3,4, when m2=0, which implies q10=±2,0,−1. Therefore, q10=−2+ax∗−gx∗y∗/x∗+b2+r2y∗2/y∗+mh2≠±2. We only require q10≠0,1, i.e.,(47)ax∗−gx∗y∗x∗+b2+r2y∗2y∗+mh2≠2,3.Letm2=0, σ=−q10/2, and θ=4q20−q120/2, and construct the nonsingular matrix.(48)L=ρ20σ−ρ1θ,Use the translationpnqn=LPnQn. Thus, system (30) takes the following form:(49)Pn+1=σPn+θQn+R1Pn,Qn+opn+qn3,Qn+1=−θPn+σQn+R2Pn,Qn+opn+qn3,where(50)R1Pn,Qn=1ρ2ρ3ρ22+ρ4ρ2σ−ρ1+ρ5σ−ρ12Pn2+ρ4ρ2θ+2θρ5σ−ρ1PnQn+ρ5θ2Qn2+ρ6ρ23+ρ7ρ22σ−ρ1+ρ8ρ2σ−ρ12+ρ9σ−ρ13Pn3+ρ7ρ22+2θρ8ρ2σ−ρ1+3θρ9σ−ρ1σ−ρ12Pn2Qn+θ2ρ8ρ2+3θ2ρ9σ−ρ1PnQn2+θ3ρ9Qn3,R2Pn,Qn=1ρ2θρ22ρ3ρ1−σ+ρ2η3+ρ2σ−ρ1(ρ4ρ1−σ+ρ2η4)+σ−ρ12ρ5ρ1−σ+ρ2η5Pn2+θρ2(ρ4×ρ1−σ+ρ2η4)+2θσ−ρ1ρ5ρ1−σ+ρ2η5PnQn+θ2ρ5ρ1−σ+ρ2η5Qn2+ρ23ρρ1−σ+ρ2η6+ρ22σ−ρ1ρ7ρ1−σ+ρ2η7+ρ2σ−ρ12(ρ8ρ1−σ+ρ2η8)σ−ρ13ρ9ρ1−σ+ρ2η9Pn3+θρ22(ρ7ρ1−σ+ρ2η7)+2θρ2σ−ρ1ρ8ρ1−σ+ρ2η8+3θσ−ρ12×ρ9ρ1−σ+ρ2η9Pn2Qn+θ2ρ2ρ8ρ1−σ+ρ2η8+3θ2σ−ρ1ρ9ρ1−σ+ρ2η9PnQn2+θ3(ρ9ρ1−σ+ρ2η9)Qn3.Next, we require the nonzero quantityχ∗ to ensure that system (34) undergoes NS bifurcation.(51)χ∗=−Re1−2λλ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(52)ξ20=18R1PnPn−R1QnQn+2R2PnQn+iR2PnPn−R2QnQn−2R1PnQn,ξ11=14R1PnPn+R1QnQn+iR2PnPn+R2QnQn,ξ02=18R1PnPn−R1QnQn−2R2PnQn+iR2PnPn−R2QnQn+2R1PnQn,ξ21=116R1PnPnPn+R1PnQnQn+R2PnPnQn+R2QnQnQn+iR2PnPnPn+R2PnQnQn−R1PnPnQn−R1QnQnQn.Therefore, from [42], we can state the subsequent results.Theorem 2.
If (32) and (33) hold and also the quantity χ∗ is nonzero, then system (28) admits NS bifurcation at E∗x∗,y∗ when mh changes in the neighborhood of Θh. Additionally, if the quantity χ∗<0 (or χ∗>0), then the stable (or unstable) invariant closed curve starts to bifurcate from E∗x∗,y∗.


## 4.1. Flip Bifurcation

System (5) undergoes flip bifurcation if one of the eigenvalues ofJE∗ must be −1 and the other should not be 1 or −1. Using this, we assume one of the eigenvalues of JE∗ as −1. Then, from (12), we have the following:(15)F−1:=4−2ax∗−2r2y∗2y∗+m2+2gx∗y∗x∗+b2+ar2x∗y∗2y∗+m2+ghx∗y∗2x∗+bx∗+c2−gr2x∗y∗3x∗+b2y∗+m2=0.Note that the left-hand side of expression (13) is in terms of parameter m in denominator with power. Hence, it is difficult to find the explicit value of m for which system (5) undergoes flip bifurcation. Hence, we denote m=mh as the critical parameter value that satisfies (13). Also, E∗x∗,y∗ is derived at the critical value m=mf.Next, ifa2>0, a3<0, and hy∗+m/r2>c holds, we define that the neighborhood Θf for system (5) undergoes flip bifurcation near the equilibrium point E∗ at some critical value m=mf as follows:(16)Θf=r1,r2,a,b,c,g,h,m:m=mf,r1,r2,a,b,c,g,h>0.Hence, we assume that system (5) undergoes flip bifurcation when m=mf changes in the neighborhood Θf.Now, we investigate the direction and stability of the possible occurrence of flip bifurcation for system (5) atE∗x∗,y∗ with the steps followed in [32, 42]. Since, r1,r2,a,b,c,g,h,m∈Θf on giving a perturbation m1≪1 of critical value m=mf for which system (5) undergoes flip bifurcation, then the perturbation system is given by the following:(17)xn+1=xnexpr1−axn−gynxn+b,yn+1=ynexpr2ynyn+m1+mf−hynxn+c.Next, we transform the equilibrium pointE∗x∗,y∗ into the origin by letting pn=xn−x∗ and qn=yn−y∗. Form (17), we obtain the following:(18)pn+1=pn+x∗expr1−apn+x∗−gqn+y∗pn+x∗+b−x∗,qn+1=qn+y∗expr2qn+y∗qn+y∗+m1+mf−hqn+y∗pn+x∗+c−y∗.By Taylor’s series expansion, (16) becomes,(19)pn+1=α1pn+α2qn+α3m1+α4pn2+α5qn2+α6m12+α7pnqn+α8pnm1+α9qnm1+opn+qn+m12,qn+1=β1pn+β2qn+β3m1+β4pn2+β5qn2+β6m12+β7pnqn+β8pnm1+β9qnm1+opn+qn+m12,where(20)α1=1+x∗M1,α2=−gx∗x∗+b′,α3=0,α4=−gx∗y∗x∗+b+M1x∗2M12,α5=g2x∗2x∗+b2,α6=0,α7−gbx∗+b2−gM1x∗x∗+b,α8=0,α9=0,β1=hy∗2x∗+c2,β2=1+y∗M2,β3=−r2y∗2y∗+mf2,β4=−hy∗2x∗+c3+h2y∗32x∗+c4,β5=−r2mfy∗y∗+mf3+M2+y∗M222,β6=r22y∗32y∗+mf4+r2y∗2y∗+mf3,β7=2hy∗x∗+c2+hy∗2M2x∗+c2,β8=−hr2y∗3x∗+c2y∗+b2,β9=−2mfr2y∗y∗+mf3−r2M2y∗2y∗+mf2,M1=−α+gy∗x∗+b2,M2=−r2y∗y∗+mf2.Consider(21)JE∗=1−ax∗+gx∗y∗x∗+b2−gx∗x∗+bhy∗2x∗+c21−r2y∗2y∗+mf2.Let us assume that the matrixJE∗ has eigenvalues.(22)λ1=−1,λ2=1−ax∗+gx∗y∗x∗+b2+ghx∗y∗2y∗+mf2x∗+bx∗+c22y∗+mf2−r2y∗,satisfiesλ1=1 and λ2≠1.Next, we construct the nonsingular matrixL as follows:(23)L=α2α2−1−α1λ2−α1.Using translationpnqn=LPnQn. Then, (17) can be written as(24)Pn+1=−Pn+G1pn,qn,m1+opn+qn+m12,Qn+1=λ2Qn+G2pn,qn,m1+opn+qn+m12,where(25)G1pn,qn,m1=ϕ1m1+ϕ2pn2+ϕ3qn2+ϕ4m12+ϕ5pnqn+ϕ6pnm1+ϕ7qnm1,G2pn,qn,m1=θ1m1+θ2pn2+θ3qn2+θ4m12+θ5pnqn+θ6pnm1+θ7qnm1,and(26)ϕ1=λ2−α1α3−α2β3α21+λ2,ϕ2=λ2−α1α4−α2β4α21+λ2,ϕ3=λ2−α1α5−α2β5α21+λ2,ϕ4=λ2−α1α6−α2β6α21+λ21+λ2,ϕ5=λ2−α1α7−α2β7α21+λ2,ϕ6=λ2−α1α8−α2β8α21+λ2,ϕ7=λ2−α1α9−α2β9α21+λ2,θ1=1+α1α3+α2β3α21+λ2,θ2=1+α1α4+α2β4α21+λ2,θ3=1+α1α5+α2β5α21+λ2,θ4=1+α1α6+α2β6α21+λ2,θ5=1+α1α7+α2β7α21+λ2,θ6=1+α1α8+α2β8α21+λ2,θ7=1+α1α9+α2β9α21+λ2.Now, by the center manifold theorem, we obtain the center manifoldGc0,0 of system (19) at 0,0 in a small neighborhood of m1=0 as follows:(27)Gc0,0=Pn,Qn:Qn=hPn,m1,=Pn,Qn:Qn=e1m1+e2Pn2+e3m1Pn+e4m12+oPn+m12.On applying map (19) on both sides of Qn=hPn,m1, we have the following:(28)λ2Qn+G2pn,qn,m1=e1m1+e2−Pn+G1pn,qn,m12+e3m1−Pn+G1pn,qn,m1+e4m12+oPn+m12,where(29)pn=α2Pn+Qn=α2Pn+hPn,m1,qn=−1−α1Pn+λ2−α1Qn=−1−α1Pn+λ2−α1hPn,m1,e1=θ11−λ2,e2=11−λ2θ2α22+θ31+α12−θ5α2−1−α1,e3=11+λ2−2e2ϕ1−θ6α2+θ7−1−α1−2e1θ2α22+2e1θ3−1−α1λ2−α1−e1θ5α2β2−α1,e4=11−λ2e12θ2α22+e12θ3λ2−α12+e12θ5α2λ2−α1+e1θ6α2+e1θ7λ2−α1−e2ϕ12+θ4−ϕ1e3.Therefore, on the center manifoldMc at origin, we have the following:(30)pn2=β22Pn2+2PnQn+Qn2,pnqn=−β21+β1Pn2+β2δ2−β1PnQn+β2λ2−β1Qn2,qn2=1+β12Pn2−21+β1λ2−β1PnQn+λ2−β12Qn2,where(31)PnQn=e1m1Pn+e2Pn3+e3m1Pn2+e4m12Pn+oPn+m13,Qn2=e12m12+2e1e2m1Pn2+2e1e3m12Pn+2e1e4m13+oPn+m13.Moreover, the mapG∗ restricted to CM Gc0,0 is as follows:(32)G∗Pn=−Pn+G1pn,qn,m1=−Pn+g1m1+g2Pn2+g3Pnm1+g4m12+g5Pn2m1+g6Pnm12+g7Pn3+g8m13+oPn+m13,where(33)g1=ϕ1,g2=ϕ2β22+ϕ31+β12−ϕ5β21+β1,g3=2e1ϕ2β22−2e1ϕ31+β1λ2−β1+e1ϕ5β2λ2−β1+ϕ6β2−ϕ71+β1,g4=e12ϕ2β22+e12ϕ3λ2−β12+ϕ4+e12ϕ5β2λ2−β1+e1ϕ6β2+e1ϕ7λ2−β1,g5=2e3ϕ2β22+2e1e2ϕ2β22−2e3ϕ31+β1λ2−β1+2e1e2ϕ3λ2−β12+e3ϕ5β2λ2−β1+2e1e2ϕ5β2λ2−β1+e2ϕ6β2+e2ϕ7λ2−β1,g6=2e4ϕ2β22+2e1e3ϕ2β22−2e4ϕ31+β1λ2−β1+2e1e3ϕ3λ2−β12+e4ϕ5β2λ2−β1+2e1e3ϕ5β2λ2−β1+e3ϕ6β2+e3ϕ7λ2−β1,g7=2e2ϕ2β22−2e2ϕ31+β1λ2−β1+e2ϕ5β2λ2−β1,g8=2e1e4ϕ2β22+2e1e4ϕ3λ2−β12+2e1e4ϕ5β2λ2−β1+e4ϕ6β2+e4ϕ7λ2−β1.From [43], we define ζ1 and ζ2 as follows:(34)ζ1=GPnm1∗+12Gm1∗GPnPn∗|Pn,m1=0,0=g3+g1g2,(35)ζ2=16GPnPnPn∗+12GPnPn∗2|Pn,m1=0,0=g7+g22.Therefore, we have got the following findings of flip bifurcation from the aforementioned study.Theorem 1.
Ifζ1≠0 and ζ2≠0, then system (5) exhibits a flip bifurcation at E∗x∗,y∗ when m1 changes in a small neighborhood of the origin. Moreover, if ζ2>0 (or ζ2<0), then the existing period-two orbits from E∗x∗,y∗ are stable (or unstable).


## 4.2. Neimark–Sacker Bifurcation

In this subsection, we find the condition for the existence of Neimark–Sacker (NS) bifurcation and also its direction and stability properties for system (5) nearE∗x∗,y∗ using the method followed in [32, 42]. For the existence of NS bifurcation of system (5), the complex conjugate eigenvalues of the characteristic equation at E∗x∗,y∗ should have absolute value one. For this, it is necessary to satisfy the following:(36)B1x∗,y∗2−4B2x∗,y∗<0andB2x∗,y∗=1,which provides(37)Am:=2−ax∗+gx∗y∗x∗+b2−r2y∗2y∗+m22−41−ax∗+gx∗y∗x∗+b2−r21−ax∗y∗2y∗+m2−gr2x∗y∗3x∗+b2y∗+m2+ghx∗y∗2x∗+bx∗+c2<0,Bm:=−ax∗+gx∗y∗x∗+b2−r21−ax∗y∗2y∗+m2−gr2x∗y∗3x∗+b2y∗+m2+ghx∗y∗2x∗+b2x∗+c2=0.Note that the left-hand side of (37) is in terms of m since m is taken as a bifurcation parameter. It is difficult to find the explicit expression for the critical value m to satisfy (24) for the occurrence of NS bifurcation of system (5). Therefore, we assume that the critical value is m=mh. Also, the interior equilibrium E∗x∗,y∗ is calculated at the critical value m=mh.Next, ifa2>0, a3<0, and hy∗+m/r2>c holds, we define the neighborhood as follows:(38)Θh=r1,r2,a,b,c,g,h,m:m=mh,r1,r2,a,b,c,g,h>0.Then, the equilibrium pointE∗x∗,y∗ can arise NS bifurcation at m=mh when it changes in the neighborhood of Θh. Now, we analyze the properties of possible NS bifurcation at E∗x∗,y∗ for system (5) if (24) holds for some mh. Given a perturbation m2≪1 of critical value mh, then the perturbation system is described as follows:(39)xn+1=xnexpr1−axn−gynxn+b,yn+1=ynexpr2ynyn+m2+mh−hynxn+c.Let us use the transformpn=xn−x∗, qn=yn−y∗ and shift E∗x∗,y∗ to 0,0. Then, system (28) takes the following form:(40)pn+1=pn+x∗expr1−apn+x∗−gqn+y∗E∗x∗,y∗+b−x∗,qn+1=qn+y∗expr2qn+y∗qn+y∗+m2+mh−hqn+y∗E∗x∗,y∗+c−y∗.Therefore, by Taylor expansion of (40),(41)pn+1=ρ1pn+ρ2qn+ρ3pn2+ρ4pnqn+ρ5qn2+ρ6pn3+ρ7pn2qn+ρ8pnvn2+ρ9vn3+opn+qn3,qn+1=η1pn+η2qn+η3pn2+η4pnqn+η5qn2+η6pn3+η7pn2qn+η8pnvn2+η9vn3+opn+qn3,where(42)ρ1=1+x∗N1,ρ2=−gx∗x∗+b,ρ3=−gx∗y∗x∗+b3+N1+x∗N122,ρ4=−gbx∗+b2−gN1x∗x∗+b,ρ5=g2x∗2x∗+b2,ρ6=−bgy∗x∗+b4−gN1x∗y∗x∗+b3+N122+x∗N136,ρ7=g2x∗y∗x∗+b4+gbx∗+b3−gbN1x∗+b2−gN12x∗2x∗+b,ρ8=−g2x∗x∗+b3+g21+x∗N12x∗+b2,ρ9=−g3x∗6x∗+b3,η1=hy∗2x∗+c2,η2=1+y∗N2,η3=−hy∗2x∗+c3+h2y∗32x∗+c4,η4=2hy∗x∗+c2+hy∗2N2x∗+c2,η5=−r2mhy∗y∗+mh3+N2+y∗N222,η6=hy∗2x∗+c4−h2y∗3x∗+c5+h3y∗46x∗+c6,η7=−hy∗2+y∗N2x∗+c3+h2y∗23+y∗N22x∗+c4,η8=h1+y∗2N222x∗+c2+2hy∗N2x∗+c2−r2mhy∗+mh2,η9=r2mhy∗y∗+mh4−r2mh1+N2y∗y∗+mh3+N222+N23y∗6,N1=−a+gy∗x∗+b2,N2=−r2y∗y∗+mh2.The characteristic polynomial equation associated with the linearized system (30) at the origin can be given as follows:(43)λ2+q1m2λ+q2m2=0,where(44)q1m2=−1+ax∗−gx∗y∗x∗+b2Ω1−1−r2y∗2y∗+mh+m22Ω2,q2m2=1−ax∗+gx∗y∗x∗+b21−r2y∗2y∗+mh+m22+ghx∗y∗2x∗+bx∗+c2Ω1Ω2,WithΩ1=expr1−ax∗−gy∗/x∗+b and Ω2=expr2y∗/y∗+mh+m2−hy∗/x∗+c. Now, the roots of (43) are the pair of complex conjugates.(45)λ1,2=12−q1m2±i4q2m2−q1m22.Sincer1,r2,a,b,c,g,h,m∈Θh, we have,λ1,2=q2m2and(46)dλ1,2dm2=12q201−ax∗+gx∗y∗x∗+b22r2y∗2y∗+mh3−r2y∗y∗+mh21−ax∗+gx∗y∗x∗+b2×1−r2y∗2y∗+mh2+ghx∗y∗2x∗+bx∗+c2≠0.Furthermore, it is requiredλ1k,λ2k≠1 for k=1,2,3,4, when m2=0, which implies q10=±2,0,−1. Therefore, q10=−2+ax∗−gx∗y∗/x∗+b2+r2y∗2/y∗+mh2≠±2. We only require q10≠0,1, i.e.,(47)ax∗−gx∗y∗x∗+b2+r2y∗2y∗+mh2≠2,3.Letm2=0, σ=−q10/2, and θ=4q20−q120/2, and construct the nonsingular matrix.(48)L=ρ20σ−ρ1θ,Use the translationpnqn=LPnQn. Thus, system (30) takes the following form:(49)Pn+1=σPn+θQn+R1Pn,Qn+opn+qn3,Qn+1=−θPn+σQn+R2Pn,Qn+opn+qn3,where(50)R1Pn,Qn=1ρ2ρ3ρ22+ρ4ρ2σ−ρ1+ρ5σ−ρ12Pn2+ρ4ρ2θ+2θρ5σ−ρ1PnQn+ρ5θ2Qn2+ρ6ρ23+ρ7ρ22σ−ρ1+ρ8ρ2σ−ρ12+ρ9σ−ρ13Pn3+ρ7ρ22+2θρ8ρ2σ−ρ1+3θρ9σ−ρ1σ−ρ12Pn2Qn+θ2ρ8ρ2+3θ2ρ9σ−ρ1PnQn2+θ3ρ9Qn3,R2Pn,Qn=1ρ2θρ22ρ3ρ1−σ+ρ2η3+ρ2σ−ρ1(ρ4ρ1−σ+ρ2η4)+σ−ρ12ρ5ρ1−σ+ρ2η5Pn2+θρ2(ρ4×ρ1−σ+ρ2η4)+2θσ−ρ1ρ5ρ1−σ+ρ2η5PnQn+θ2ρ5ρ1−σ+ρ2η5Qn2+ρ23ρρ1−σ+ρ2η6+ρ22σ−ρ1ρ7ρ1−σ+ρ2η7+ρ2σ−ρ12(ρ8ρ1−σ+ρ2η8)σ−ρ13ρ9ρ1−σ+ρ2η9Pn3+θρ22(ρ7ρ1−σ+ρ2η7)+2θρ2σ−ρ1ρ8ρ1−σ+ρ2η8+3θσ−ρ12×ρ9ρ1−σ+ρ2η9Pn2Qn+θ2ρ2ρ8ρ1−σ+ρ2η8+3θ2σ−ρ1ρ9ρ1−σ+ρ2η9PnQn2+θ3(ρ9ρ1−σ+ρ2η9)Qn3.Next, we require the nonzero quantityχ∗ to ensure that system (34) undergoes NS bifurcation.(51)χ∗=−Re1−2λλ¯21−λξ11ξ20−12ξ112−ξ022+Reλ¯ξ21,where(52)ξ20=18R1PnPn−R1QnQn+2R2PnQn+iR2PnPn−R2QnQn−2R1PnQn,ξ11=14R1PnPn+R1QnQn+iR2PnPn+R2QnQn,ξ02=18R1PnPn−R1QnQn−2R2PnQn+iR2PnPn−R2QnQn+2R1PnQn,ξ21=116R1PnPnPn+R1PnQnQn+R2PnPnQn+R2QnQnQn+iR2PnPnPn+R2PnQnQn−R1PnPnQn−R1QnQnQn.Therefore, from [42], we can state the subsequent results.Theorem 2.
If (32) and (33) hold and also the quantity χ∗ is nonzero, then system (28) admits NS bifurcation at E∗x∗,y∗ when mh changes in the neighborhood of Θh. Additionally, if the quantity χ∗<0 (or χ∗>0), then the stable (or unstable) invariant closed curve starts to bifurcate from E∗x∗,y∗.


## 5. Chaos Control

In this section, we study the chaos control analyses for system (5). Firstly, we use the state feedback control method as in [27] to control the chaotic system. For system (5), we consider the following corresponding controlled system:(53)xn+1=xnexpr1−axn−gynxn+b−uxn,yn,yn+1=ynexpr2ynyn+m−hynxn+c,where uxn,yn=h1xn−x∗+h2yn−y∗ is the feedback controlling force with feedback gains h1 and h2, with equilibrium E∗x∗,y∗ of (53). For the controlled system (36), the Jacobian matrix at E∗x∗,y∗ is given as(54)Jx∗,y∗=κ11−h1κ12−h2κ21κ22,where(55)κ11=1−ax∗+gx∗y∗x∗+b2,κ12=−gx∗x∗+b,κ21=hy∗2x∗+b2,κ22=1−r2y∗2y∗+m2.Then, forJx∗,y∗, we have the characteristic equation as follows:(56)λ¯2−κ11+κ22−h1λ¯+κ22κ11−h1−κ21κ12−h2=0.Let λ¯1 and λ¯2 be the eigenvalues of (54), which yields the following:(57)λ¯1+λ¯2=κ11+κ22−h1,(58)λ¯1λ¯2=κ22κ11−h1−κ21κ12−h2.Using (57) and (58), λ¯1λ¯2=1, λ¯1=1, and λ¯1=−1. Also, ensure that |λ¯1,2|<1. Then, we derive the marginal stability lines as follows:(59)L¯1:κ11κ22−κ21κ12−1=h1κ22−h2κ21,L¯2:h11−κ22+h2κ21=κ11+κ22−1−κ11κ22+κ21κ12,L¯3:h11+κ22−h2κ21=κ11+κ22+1+κ11κ22−κ21κ12.Moreover, the triangular region enclosed by linesL¯1, L¯2, and L¯3 have stable eigenvalues for the Jacobian matrix (37).Next, the pole-placement control method, as in [31], is utilized to control the unstable dynamics of system (5). By taking m as control parameter, we rewrite system (5) as follows:(60)xn+1=xnexpr1−axn−gynxn+b=f1xn,yn,yn+1=ynexpr2ynyn+m−hynxn+c=f2xn,yn.Furthermore,m is needed to lie in some interval |m−m0|<δ with δ>0 and m0 denotes the nominal value, for which system (5) has unstable dynamics. Now, we utilize the state feedback control method to shift the trajectory to the expected state. Let the equilibrium point E∗x∗,y∗ be unstable for system (5) because of NS bifurcation. Then, system (5) can be approximated near E∗x∗,y∗ using the linear map, which is given by the following:(61)xn+1−x∗yn+1−y∗≈Axn−x∗yn−y∗+Bm−m0,where(62)A=∂f1x∗,y∗,m0∂xn∂f1x∗,y∗,m0∂yn∂f2x∗,y∗,m0∂xn∂f2x∗,y∗,m0∂yn,B=∂f1x∗,y∗,m0∂m∂f2x∗,y∗,m0∂m=0−r2y∗y∗+m2.It is clear to observe that system (44) is controllable, provided matrix C has rank 2, which is,(63)C=B:AB=0∂f1x∗,y∗,m0∂yn−r2y∗2y∗+m2−r2y∗y∗+m2∂f2x∗,y∗,m0∂yn−r2y∗y∗+m2.Furthermore,−r2y∗2/y∗+m2≠0, and assume that ∂f2x∗,y∗,m0/∂yn≠0. Then, system (44) is controllable. Next, we assume that m−m0=−Kxn−x∗yn−y∗, where K=s1s2. Then, system (44) can be written as follows:(64)xn+1−x∗yn+1−y∗≈A−BKxn−x∗yn−y∗.Then, the controller system is given by the following:(65)xn+1=xnexpr1−axn−gynxn+b=f1xn,yn,yn+1=ynexpr2ynm0−S+yn−hynxn+c=f2xn,yn.where S=s1xn−x∗+s2yn−y∗. Furthermore, if the matrix A−BK has eigenvalues λ¯1 and λ¯2 that lie in an open unit disk, then E∗x∗,y∗ is locally asymptotically stable. Thus, we have the following:(66)A−BK=κ11κ12κ21−θs1κ22−θs1,where θ=−r2y∗2/y∗+m2, and κ11, κ12, κ21, and κ22 are the same as in (54). Then, the characteristic polynomial of (66) can be written as follows:(67)λ¯2−κ11+κ22−θs2λ¯+κ11κ22−θs2+κ12θs1−κ21=0.The lines of marginal stability are obtained, which are as follows:(68)L¯1:κ11κ22−θs2+κ12θs1−κ21=1,L¯2:κ11+κ22=1+θs2+κ11κ22−θs2+κ12θs1−κ21,L¯3:θs2=κ11+κ22+1+κ11κ22−θs2+κ12θs1−κ21.Therefore, matrix (49) has stable eigenvalues enclosed by the straight lines L¯1, L¯2, and L¯3 in s1,s2- plane.Next, we apply the hybrid control feedback methodology [28, 29] for controlling the bifurcation behavior of the system near the equilibrium point E∗x∗,y∗. Then, the controlled system can be written as follows:(69)xn+1=ϵxnexpr1−axn−gynxn+b+1−ϵxn,yn+1=ϵynexpr2ynyn+m−hynxn+c+1−ϵyn,where 0<ϵ<1 is the controlled strategy of the combination of both feedback control and parameter perturbation. The Jacobian matrix evaluated for system (54) at E∗x∗,y∗ is given by the following:(70)1−ϵax∗+ϵgx∗y∗x∗+b2−ϵgx∗x∗+bϵhy∗2x∗+C21−ϵr2y∗2y∗+m2.Note that one can select the appropriate value forϵ to ensure that all eigenvalues of the above matrix satisfy |λ¯1,2|<1.Remark 1.
System (2) is the extension of a continuous-time system (1) with the Allee effect in the predator’s growth term. It is worth mentioning here that the discrete form of system (2) is not studied elsewhere in the literature. Hence, we obtained discrete-time system (5) from (2) using methods similar to those in [11, 12]. Moreover, the influence of the Allee effect is shown in terms of flip and Neimark–Sacker bifurcations for system (5). Using the results in [42, 43], the direction and stability properties of both bifurcations can be discussed with the help of Theorem 1 and 2. Also, the bifurcation and chaos control analyses are carried out by utilizing the methods in [27–31].


## 6. Numerical Simulations

In the subsequent section, we perform some simulations for system (5) at E∗x∗,y∗ to ensure our mathematical results obtained in the previous sections.Case (i):.
firstly, let us take the parameter values asr1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, and c=0.5 and varying m∈0,1 for system (5). Then, system (5) exhibits flip bifurcation when m reaches the critical value mf=0.28805, and the interior equilibrium point at mf is E∗x∗,y∗=0.76387948,2.11933949, which also satisfies (13). Also, the characteristic (12) is given by the following:(71)λ2+1.11345044λ+0.11345926=0,where λ1,2=−0.11346,−1 are the eigenvalues of the Jacobian matrix at E∗x∗,y∗|mf. Then, from (34) and (35), we obtain ζ1=1.3842 and ζ2=5.17693, and the properties of flip bifurcation are illustrated in Theorem 1.
The interior equilibrium point varies accordingly for various values ofm, which is given in Table 1. The nature of system (5) near E∗x∗,y∗ is shown for different values of m, for m=0.2 is stable, m=0.27 is period-2, m=0.4 is period-4, m=0.5 is period-8, m=0.53 is chaotic, and m=0.9 is chaotic, which are shown by the time series plots in Figure 1 and the phase portraits in Figure 2. It shows that the system becomes chaotic from stable via period-doubling cascade for larger values of m. Also, the chaotic nature of the system is confirmed by one parameter bifurcation diagram and the maximum Lyapunov exponent in Figure 3.Table 1
The equilibrium point values of Case (i).
mx∗y∗0.20.755 2532.190 9580.270.762 1142.134 0270.40.774 7892.028 1700.50.784 4811.946 6310.530.787 3791.922 1510.90.822 7601.619 543Figure 1
The time series plots for system (5) withr1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, and c=0.5 (a) locally stable for m=0.2, (b) period-2 for m=0.27, (c) period-4 for m=0.4, (d) period-8 for m=0.5, (e) 8-chaotic like attractor for m=0.53, and (f) chaotic attractor for m=0.9.
(a)(b)(c)(d)(e)(f)Figure 2
The phase portraits for system (5) with r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, and c=0.5 (a) locally stable for m=0.2, (b) period-2 for m=0.27, (c) period-4 for m=0.4, (d) period-8 for m=0.5, (e) 8-chaotic like attractor for m=0.53, and (f) chaotic attractor for m=0.9. The red points represent the interior equilibrium point E∗x∗,y∗, and its values are given in Table 1.
(a)(b)(c)(d)(e)(f)Figure 3
((a), (b)). The flip bifurcation diagrams. (c) The maximum Lyapunov exponents for system (5) with r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, and m∈0,1.
(a)(b)(c)Case (ii):.
next, choose the system parameter values asr1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, and varying m∈0,1. Then, system (5) exhibits NK bifurcation at the critical value m=mh=0.199023 and has the equilibrium point E∗x∗,y∗=0.69031513,2.24467052. Then, the characteristic (12) takes the following form:(72)λ2−1.53322562λ+1=0,Then, the above characteristic polynomial has the eigenvalues λ1,2=0.766613±i0.64211. From (51), we obtain χ∗=−0.0946108. Then, the properties of NK bifurcation are illustrated in Theorem 2.
The interior equilibrium point varies for various values ofm, which is given in Table 2. The nature of system (5) near E∗x∗,y∗ is shown for different values of m, for m=0.18 is asymptotically stable, m=0.2,0.5 are invariant circles, m=0.58 is period-8, m=0.88 is invariant circle, and m=0.95 is locally asymptotically stable, as shown by the time plots in Figure 4 and the phase portraits in Figure 5. It shows that system (5) undergoes NK bifurcation and again attains stability via NK bifurcation. Also, the one parameter bifurcation diagram and the maximum Lyapunov exponent are shown in Figure 6.Table 2
The equilibrium point values of Case (ii).
mx∗y∗0.180.655 4892.230 8390.20.692 0542.245 3340.51.099 5892.329 8010.581.184 1782.329 6020.881.457 3702.287 3300.91.513 8872.270 648Figure 4
The time series plots for system (5) with r1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, and c=1.9 (a) asymptotically stable for m=0.18, ((b), (c)) invariant circles for m=0.2,0.5, (d) period-8 for m=0.58, (e) invariant circle for m=0.88, and (e) asymptotically stable for m=0.95.
(a)(b)(c)(d)(e)(f)Figure 5
The phase portraits for system (5) with r1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9 (a) locally asymptotically stable for m=0.18; ((b), (c)) invariant circles for m=0.2,0.5, (d) period-8 for m=0.58, (e) invariant circle for m=0.88, and (f) asymptotically stable for m=0.95. The red points represent the interior equilibrium point E∗x∗,y∗, and its values are given in Table 2.
(a)(b)(c)(d)(e)(f)Figure 6
((a), (b)).The Neimark–Sacker bifurcation diagram. (c) The maximum Lyapunov exponents for system (5) with r1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, and m∈0,1.
(a)(b)(c)


### 6.1. Chaos Control

Next, we apply the state feedback control method for existing chaos in the flip bifurcation. The system (5) is in chaotic state for the parameter values in case (i),r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5 and m=0.53. Then, E∗x∗,y∗=0.787379,1.922151 is unstable, and its phase portrait is given in Figure 2(e). We need to shift the unstable equilibrium toward the stable state. For this, take m=0.53, and the corresponding controlled system is given by the following:(73)xn+1=xnexp3.3−3.3xn−1.2ynxn+2.5−uxn,yn,yn+1=ynexpynyn+0.53−0.525ynxn+0.5,where uxn,yn=h1xn−x∗+h2yn−y∗, and h1, h2 are feedback gains. Furthermore, system (56) has the Jacobian matrix as follows:(74)J=−1.4303−h1−0.287419−h21.170370.385558.Then, we obtain the characteristic equation as follows:(75)λ¯2+1.04474+h1λ¯−0.385558h1+1.17037h2−0.215078=0.Furthermore, system (56) has the marginal stability lines, which are given by the following:(76)L¯1:=0.38558h1−1.17037h2=−1.21508,L¯2:=0.614442h1+1.17037h2=−1.82966,L¯3:=1.38556h1−1.17037h2=−0.259816.Then, the controlled system (56) has eigenvalues that may then be shown to be located within the triangular region defined by the straight lines L¯1, L¯2, and L¯3 (see Figure 7(a)). On choosing k1=−1, E∗x∗,y∗=0.787379,1.922151 is locally stable if and only if h2∈−1,0.8. Take h1=−1 and h2∈0,1.4. Then, the bifurcation diagram for system (56) is plotted in Figures 7(b) and 7(c).Figure 7
For the controlled system (36) with r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, m=0.53, and E∗x∗,y∗=0.787379,1.922151. (a) The stability triangle and ((b), (c) the bifurcation diagrams with h1=−1 and h2=0,1.4.
(a)(b)(c)Secondly, we utilize the pole-placement control method by taking parameter values as in case (ii)r1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, and m=0.5. Then, system (5) is unstable near the equilibrium point E∗x∗,y∗=1.099589,2.329801. To shift the unstable trajectory to the desired stable state, take m=0.5, and system (48) is given by the following:(77)xn+1=xnexp4−1.15xn−2.7ynxn+1.2=f1xn,yn,yn+1=ynexpyn0.5−S+yn−1.06ynxn+1.9=f2xn,yn,where S=s1xn−x∗+s2yn−y∗, s1 and s2 are feedback gains. Then, we have the following:(78)A=1.04349−1.291050.639470.322162,B=0−0.677838,C=00.875124−0.677838−0.218374,and the above-controlled system has the Jacobian matrix, which is of the following form:(79)A−BK=1.04349−1.291050.63947+0.677838s10.322162+0.677838s2,and its characteristic polynomial is written as follows:(80)λ¯2−1.36565+0.677838s2λ¯+0.875124s1+0.707315s2+1.16176=0.Furthermore, the lines of marginal stability of (77) are computed as follows:(81)L1:=1.16176+0.875124s1−0.707315s2=1,L2:=2.16176+0.875124s1+0.0294771s2=1.36565,L3:=3.52741+0.875124s1+0.707315s2=−0.677838s2.The eigenvalues may then be shown to be located within the triangular region defined by the straight linesL¯1, L¯2, and L¯3 for system (57) (see Figure 8(a)). On choosing s2=−2, E∗x∗,y∗ is locally asymptotically stable if and only if s1∈−0.84,1.43. Take s1 and s2∈0,2. Then, the bifurcation diagrams for the controlled system (57) are depicted in figures 8(b) and 8(c).Figure 8
For the controlled system (48) withr1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, m=0.5, and E∗x∗,y∗=1.099589,2.329801. (a) The stability triangle and ((b), (c) depict the bifurcation diagrams with s2=−2 and s1=0,2.
(a)(b)(c)Finally, we choose the parameter values as in case (i)r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, and m=0.53 to explore the hybrid control strategy. Now, the controlled system (54) takes the following form:(82)xn+1=ϵxnexp3.3−3.3xn−1.2ynxn+2.5+1−ϵxn,yn+1=ϵynexpynyn+0.53−0.525ynxn+0.5+1−ϵyn.Moreover, we haveE∗x∗,y∗=0.787379,1.922151 for system (58), and its Jacobian matrix is written as follows:(83)1−2.4303ϵ−0.28419ϵ1.17037ϵ1−0.614442ϵ.If0<ε<0.900535, then the above matrix has stable eigenvalues. Then, system (58) is stable near E∗x∗,y∗ for 0<ϵ<0.900535, which is depicted in the bifurcation diagram in Figure 9.Figure 9
The bifurcation diagrams for the controlled system (54) withr1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, m=0.53, E∗x∗,y∗=0.787379,1.922151, and ϵ∈0,1.
(a)(b)Remark 2.
Allee effect could have a stabilizing or destabilizing effect or both, which depends on the parameters in the prey-predator model. For example, the growth of prey can cause chaotic dynamics in the discrete-time predator-prey model, and how the populations are changed from extinction to persistence because of the Allee effect in both prey and predator have been reported in [24]. The authors in [26] showed that without the Allee effect for a larger growth rate of prey, the system becomes chaotic near the boundary equilibrium point. In the presence of the Allee effect, the system becomes stable via reverse periodic doubling. Also, the considered system undergoes flip and Neimark–Sacker bifurcation near the interior equilibrium point at some critical Allee parameter value. System (5) considered in the present study has destabilizing effect for larger m, which is clearly shown in Figure 3. Also, with a different set of parameters, system (5) has both stabilizing and destabilizing effect because of m via Neimark–Sacker bifurcation, as shown in Figure 6.Remark 3.
Tassaddiq et al. [35] used a nonstandard finite difference scheme and obtained the discrete ratio-dependent prey-predator model. They showed that the considered system undergoes Neimark–Sacker bifurcation for larger values of the catchability coefficient. Moreover, they applied a pole-placement control strategy to obtain the controlled system and showed the controllable region in the feedback control space. Also, they applied the hybrid control strategy by choosing the catchability coefficient in the chaotic region and found the stability interval for the control parameter. Din [32] obtained the discrete system for (1) with a death rate of a predator by the method of piecewise constant arguments and showed that the considered system undergoes both types of bifurcations. Also, the author successfully implemented the state feedback control strategy to control FB and NSB and pole placement to control NSB. Hence, in the present study, we successfully implemented state feedback and hybrid control method to control FB, and the stable regions and bifurcation diagrams for feedback gains are shown in Figures 7 and 9 for system (5). Also, we implemented a pole-placement control method to control NSB, and it is depicted in Figure 8.


## 6.1. Chaos Control

Next, we apply the state feedback control method for existing chaos in the flip bifurcation. The system (5) is in chaotic state for the parameter values in case (i),r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5 and m=0.53. Then, E∗x∗,y∗=0.787379,1.922151 is unstable, and its phase portrait is given in Figure 2(e). We need to shift the unstable equilibrium toward the stable state. For this, take m=0.53, and the corresponding controlled system is given by the following:(73)xn+1=xnexp3.3−3.3xn−1.2ynxn+2.5−uxn,yn,yn+1=ynexpynyn+0.53−0.525ynxn+0.5,where uxn,yn=h1xn−x∗+h2yn−y∗, and h1, h2 are feedback gains. Furthermore, system (56) has the Jacobian matrix as follows:(74)J=−1.4303−h1−0.287419−h21.170370.385558.Then, we obtain the characteristic equation as follows:(75)λ¯2+1.04474+h1λ¯−0.385558h1+1.17037h2−0.215078=0.Furthermore, system (56) has the marginal stability lines, which are given by the following:(76)L¯1:=0.38558h1−1.17037h2=−1.21508,L¯2:=0.614442h1+1.17037h2=−1.82966,L¯3:=1.38556h1−1.17037h2=−0.259816.Then, the controlled system (56) has eigenvalues that may then be shown to be located within the triangular region defined by the straight lines L¯1, L¯2, and L¯3 (see Figure 7(a)). On choosing k1=−1, E∗x∗,y∗=0.787379,1.922151 is locally stable if and only if h2∈−1,0.8. Take h1=−1 and h2∈0,1.4. Then, the bifurcation diagram for system (56) is plotted in Figures 7(b) and 7(c).Figure 7
For the controlled system (36) with r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, m=0.53, and E∗x∗,y∗=0.787379,1.922151. (a) The stability triangle and ((b), (c) the bifurcation diagrams with h1=−1 and h2=0,1.4.
(a)(b)(c)Secondly, we utilize the pole-placement control method by taking parameter values as in case (ii)r1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, and m=0.5. Then, system (5) is unstable near the equilibrium point E∗x∗,y∗=1.099589,2.329801. To shift the unstable trajectory to the desired stable state, take m=0.5, and system (48) is given by the following:(77)xn+1=xnexp4−1.15xn−2.7ynxn+1.2=f1xn,yn,yn+1=ynexpyn0.5−S+yn−1.06ynxn+1.9=f2xn,yn,where S=s1xn−x∗+s2yn−y∗, s1 and s2 are feedback gains. Then, we have the following:(78)A=1.04349−1.291050.639470.322162,B=0−0.677838,C=00.875124−0.677838−0.218374,and the above-controlled system has the Jacobian matrix, which is of the following form:(79)A−BK=1.04349−1.291050.63947+0.677838s10.322162+0.677838s2,and its characteristic polynomial is written as follows:(80)λ¯2−1.36565+0.677838s2λ¯+0.875124s1+0.707315s2+1.16176=0.Furthermore, the lines of marginal stability of (77) are computed as follows:(81)L1:=1.16176+0.875124s1−0.707315s2=1,L2:=2.16176+0.875124s1+0.0294771s2=1.36565,L3:=3.52741+0.875124s1+0.707315s2=−0.677838s2.The eigenvalues may then be shown to be located within the triangular region defined by the straight linesL¯1, L¯2, and L¯3 for system (57) (see Figure 8(a)). On choosing s2=−2, E∗x∗,y∗ is locally asymptotically stable if and only if s1∈−0.84,1.43. Take s1 and s2∈0,2. Then, the bifurcation diagrams for the controlled system (57) are depicted in figures 8(b) and 8(c).Figure 8
For the controlled system (48) withr1=4, a=1.15, g=2.7, b=1.2, r2=1, h=1.06, c=1.9, m=0.5, and E∗x∗,y∗=1.099589,2.329801. (a) The stability triangle and ((b), (c) depict the bifurcation diagrams with s2=−2 and s1=0,2.
(a)(b)(c)Finally, we choose the parameter values as in case (i)r1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, and m=0.53 to explore the hybrid control strategy. Now, the controlled system (54) takes the following form:(82)xn+1=ϵxnexp3.3−3.3xn−1.2ynxn+2.5+1−ϵxn,yn+1=ϵynexpynyn+0.53−0.525ynxn+0.5+1−ϵyn.Moreover, we haveE∗x∗,y∗=0.787379,1.922151 for system (58), and its Jacobian matrix is written as follows:(83)1−2.4303ϵ−0.28419ϵ1.17037ϵ1−0.614442ϵ.If0<ε<0.900535, then the above matrix has stable eigenvalues. Then, system (58) is stable near E∗x∗,y∗ for 0<ϵ<0.900535, which is depicted in the bifurcation diagram in Figure 9.Figure 9
The bifurcation diagrams for the controlled system (54) withr1=3.3, a=3.3, g=1.2, b=2.5, r2=1, h=0.525, c=0.5, m=0.53, E∗x∗,y∗=0.787379,1.922151, and ϵ∈0,1.
(a)(b)Remark 2.
Allee effect could have a stabilizing or destabilizing effect or both, which depends on the parameters in the prey-predator model. For example, the growth of prey can cause chaotic dynamics in the discrete-time predator-prey model, and how the populations are changed from extinction to persistence because of the Allee effect in both prey and predator have been reported in [24]. The authors in [26] showed that without the Allee effect for a larger growth rate of prey, the system becomes chaotic near the boundary equilibrium point. In the presence of the Allee effect, the system becomes stable via reverse periodic doubling. Also, the considered system undergoes flip and Neimark–Sacker bifurcation near the interior equilibrium point at some critical Allee parameter value. System (5) considered in the present study has destabilizing effect for larger m, which is clearly shown in Figure 3. Also, with a different set of parameters, system (5) has both stabilizing and destabilizing effect because of m via Neimark–Sacker bifurcation, as shown in Figure 6.Remark 3.
Tassaddiq et al. [35] used a nonstandard finite difference scheme and obtained the discrete ratio-dependent prey-predator model. They showed that the considered system undergoes Neimark–Sacker bifurcation for larger values of the catchability coefficient. Moreover, they applied a pole-placement control strategy to obtain the controlled system and showed the controllable region in the feedback control space. Also, they applied the hybrid control strategy by choosing the catchability coefficient in the chaotic region and found the stability interval for the control parameter. Din [32] obtained the discrete system for (1) with a death rate of a predator by the method of piecewise constant arguments and showed that the considered system undergoes both types of bifurcations. Also, the author successfully implemented the state feedback control strategy to control FB and NSB and pole placement to control NSB. Hence, in the present study, we successfully implemented state feedback and hybrid control method to control FB, and the stable regions and bifurcation diagrams for feedback gains are shown in Figures 7 and 9 for system (5). Also, we implemented a pole-placement control method to control NSB, and it is depicted in Figure 8.


## 7. Conclusion

This article deals with a discrete-time modified Leslie–Gower system with the Allee effect in the predator population. The presence of the Allee effect in the prey-predator model can have a stabilizing or destabilizing effect. The existence criteria of biologically meaningful equilibrium points have been investigated, and their stability analysis has also been carried out. In system (1), even if the prey goes extinct, the predator can survive by changing its food habits and going for an alternative food source. Note that in system (5), the prey extinction equilibrium pointE1=0,r2c/h−m exists only if r2c/h>m, i.e., the Allee parameter m is certainly less than r2c/h. Otherwise, the predator also goes to extension in the absence of prey. In the absence of the Allee effect m=0, the prey extinction equilibrium point E1 always exists. The existence and local stability of an interior equilibrium E∗x∗,y∗ can be achieved from Lemma 1 and Lemma 3. Figures 1(a), 4(a), and 4(f) show how the impact of m ensures the long-term survival of both species. We derived the conditions for the occurrence of flip and Neimark–Sacker bifurcations for system (5). Furthermore, we have discussed the direction and stability of both bifurcations with the help of the center manifold theorem and normal form theory. We showed that the size of the prey population increases and the predator population decreases when the Allee parameter m increases in Case (i), which is given in Table 1. Also, the sizes of both prey and predator populations increase when the Allee parameter m increases in Case (ii), which is given in Table 2. We showed that system (5) becomes chaotic from stable via flip bifurcation and also changes from stable to unstable via Neimark–Sacker bifurcation in numerical simulations for different sets of model parameters, varying the Allee parameter, see Figures 3 and 6. We verified the chaotic nature of the system with the help of a one-parameter bifurcation diagram and a maximum Lyapunov exponent, which show the stabilizing and destabilizing effect of the Allee parameter m. We observed that the proposed system exhibits complex behavior and is sensitive to the choice of parameter values and initial conditions. The presence of complex and chaotic behaviors in system (5) by varying the Allee parameter m can cause both prey and predator populations to have a higher risk of extinction. The presence of chaotic behavior in system (5) affects the prediction of the sizes of both populations in future generations. Also, we showed that the state feedback, pole-placement, and hybrid control methods help to shift the unstable equilibrium to a stable state for a suitable range of control parameters in Figures 7, 8, and 9.In short, the presence of the Allee effect on the predator population has a high impact on the dynamics of the modified Leslie–Gower system with populations that have no overlap between generations, resulting in a new discrete system (5) with various dynamical behaviors. Thus, it could be interesting and meaningful to study the dynamics of the discrete predator-prey system (5) with the Allee effects in both prey and predator, also with other interaction functions. However, these terms will increase the complexity of system (5), and we will leave this as future research.



---

*Source: 6931354-2022-03-28.xml*