# Bifurcation and Chaos of a Nonlinear Discrete-Time Predator-Prey Model Involving the Nonlinear Allee Effect

Authors: Ning Song
Journal: Discrete Dynamics in Nature and Society
 (2023)

---

## Abstract

The complex dynamics of a nonlinear discretized predator-prey model with the nonlinear Allee effect in prey and both populations are investigated. First, the rigorous results are derived from the existence and stability of the fixed points of the model. Second, we establish a model with the Allee effect in prey undergoing codimension-one bifurcations (flip bifurcation and Neimark–Sacker bifurcation) and codimension-two bifurcation associated with 1 : 2 strong resonance by using center manifold theorem and bifurcation theory, and the direction of bifurcations is also evaluated. In particular, chaos in the sense of Marotto is proved at some certain conditions. Third, numerical simulations are performed to illustrate the effectiveness of the theoretical results and other complex dynamical behaviors, such as the period-3, 4, 6, 8, 9, 30, and 43 orbits, attracting invariant cycles, coexisting chaotic sets, and so forth. Of most interest is the finding of coexisting attractors and multistability. Moreover, a moderate Allee effect in predators can stabilize the dynamical behavior. Finally, the hybrid feedback control strategy is implemented to stabilize chaotic orbits existing in the model.


---

## Body

## 1. Introduction

The studies on biological mathematics indicate that discrete-time models described by different equations are more pertinent mathematical tools to formulate and simulate the dynamical properties than continuous ones, particularly when populations have nonoverlapping generations [1]. Moreover, using the discrete case is also crucial because it may reveal richer dynamical behaviors and provide more efficient computations. For instance, chaos can be observed in single-species discrete models, whereas for continuous-time models, minimum three-species models are needed [1]. The models of discrete-time take the following form:(1)Qn+1=fQn.As early as the 1970s, May [2, 3] had clearly illustrated the complex dynamics possible in simple discrete-time models, whether it is from the standpoint of ecology or a mathematical point of view, to investigate the complex dynamics of discrete models are very meaningful [4–9].The most important and fascinating topic of current research in biological mathematics is the impact of the Allee effect on population dynamics. Ecologist Allee [10] first described that the Allee effect shows a positive relationship between population density and the growth rate per capita, and the existence of it indicates that there is a minimum population size necessary for a population to sustain itself in the natural world [11, 12]. There are various kinds of biological and physical factors behind this Allee mechanism, such as complications in finding mates, lessened defenses against predators, reproductive facilitation, inbreeding depressions, social thermoregulation, environment conditioning, and so on [13, 14]. This effect has already been discovered in numerous organisms in the real world, including insects [15], birds and mammals [16], plants [17], bacteria [18], marine invertebrates [19], and plant-herbivore interactions [20]. However, the Allee effects of different species, because of their different mechanisms [21], may affect the dynamics of populations differently. It is now well understood that the Allee effect plays an important role in the conservation of endangered and exploited species. Therefore, the study of the Allee effect is important to conservation biology.In recent decades, the Allee effect in the discrete-time predator-prey model [22–28] has developed quickly. Celik and Duman [29] studied the stability of the positive equilibrium point for the discrete predator-prey model without the Allee effect:(2)xn+1=xn+rxn1−xn−axnyn,yn+1=yn+aynxn−yn,and the model with the Allee effect in the prey population:(3)xn+1=xn+rxn1−xnxnm+xn−axnyn,yn+1=yn+aynxn−yn,where xn and yn represent the densities of prey and predator populations at time n, respectively. r is the intrinsic growth rate parameter of prey and a is the conversion efficiency parameter of prey to predator. r and a are positive constants. xn/m+xn in model (3) denotes the Allee effect function [12, 13, 30], where m is the Allee constant which means the population size at which fitness is half its maximum value [12, 14] and 0<m<r/a. Thereafter, Yuan and Yang [31] derived the existence and bifurcation directions of flip bifurcation, Neimark–Sacker bifurcation, and codimension-two bifurcations for model (2) with 1 : 2 resonance by using center manifold theorem, bifurcation theory, and the normal form method. Chen et al. [32] found some similar dynamic characteristics for models (2) and (3) and the important role played by the Allee constant in model (3). On the basis of model (2), Wu and Zhao [33] took xn−d0<d<1 as the Allee effect function on prey population and analyzed the detailed bifurcation behaviors of codimension-one and codimension-two with 1 : 2 resonance of the positive fixed point. Wang et al. [34] took 1−e−εx as the mate-finding Allee effect (i.e., it was hard to find a compatible and receptive mate at low population size or density [14, 19], ε is the Allee effect constant) on prey population and considered the following model:(4)xn+1=xn+rxn1−xn1−e−εxn−axnyn,yn+1=yn+aynxn−yn.Further, they considered the Allee effects on both prey and predator populations as the model(5)xn+1=xn+rxn1−xn1−e−εxn−axnyn,yn+1=yn+aynxn−ynynu+yn,where yn/u+yn incorporates the Allee effect on predator population and u is the Allee effect constant. Bigger u reflects that the stronger the Allee effect and the predator density with respect to the prey population will be decreased. For these two models, they obtained the asymptotically stable conditions of the equilibrium points and showed the bifurcation or chaotic phenomena by numerical simulations but did not theoretically prove the existence of them. Enlightened by the abovementioned works, in this paper, we are committed to analyze the complex dynamics (such as codimension-one bifurcations, codimension-two bifurcation, and Marotto’s chaos) of models (4) and (5) in detail by both analytical and numerical methods and obtain some new points (such as the coexisting attractors and multistability). The results obtained demonstrate that the Allee effect may promote persistence for both prey and predator and remain its stable dynamical behavior.The paper is organized as follows. In Section2, the rigorous results are derived from the existence of various fixed points of models (4) and (5) and their local dynamics. In Section 3, the details about the bifurcation analysis of codimension-one and codimension-two are given. In Section 4, chaos in the sense of Marotto is proved when certain conditions are satisfied. In Section 5, the extensive results of numerical simulations are used to verify the analytical results and to display rich dynamical behaviors. In Section 6, the hybrid feedback control strategy is utilized to control the chaos caused by the bifurcations of the model (4). In Section 7, our findings are summarized.


## 2. Analysis of Fixed Points

### 2.1. Existence and Stability of Fixed Points of Model (4)

Fixed points of the model (4) are determined by solving the following nonlinear equations:(6)x=x+rx1−x1−e−εx−axy,y=y+ayx−y.For all parameter valuesr,ε, and a, the trivial fixed point Z10,0 and the axial fixed point Z21,0 always exist, they are two boundary fixed points.In the following theorem, we give the existence and the number of the positive (coexistence) fixed point inR+2.Theorem 1.
Assumer,ε, andaare positive real numbers; ifa<rε, the model (4) has a unique positive (coexistence) fixed pointZ3x∗,y∗, wherex∗=y∗and(7)r1−x∗1−e−εx∗−ax∗=0.
Ifa≥rε, the model (4) has no positive fixed point.Proof.
Let(8)Fx=r1−x1−e−εx−ax.
For convenience, we assume(9)F1x=r1−x1−e−εx,F2x=ax.
It is easy to know thatF10=F11=0, and(10)F1′x=−r1−e−εx+r1−xεe−εx,F1″x=−2rεe−εx−r1−xε2e−εx<0,for x∈0,1, so F1x is a concave function in (0, 1). F2x is the straight line crossing Z10,0 with slope a. Thus, there will be only two cases for the roots of Fx=0 as follows:(i)
IfF1′0=rε≤a, then there is no root in (0, 1) of Fx=0, in other words, no positive (coexistence) fixed point exists, as shown in Figure 1(a).(ii)
IfF1′0=rε>a, then there exist exactly one root in (0, 1) of Fx=0, in other words, exactly one positive (coexistence) fixed point Z3x∗,y∗ exists, where x∗=y∗ and r1−x∗1−e−εx∗−ax∗=0, as shown in Figure 1(b).
The Jacobian matrixJ of the model (4) evaluated at any fixed point x,y is as follows:(11)J=∑i,j=12∂i,j2,2Jij,where(12)J11=1+r−2rx1−e−εx+rx1−xεe−εx−ay,J12=−ax,J21=ay,J22=1+ax−2ay.
In equation (12), the symbol∂i,jm,nM denotes a partitioned matrix with m row blocks and n column blocks, and the (i, j)-th block is M, a smaller matrix, all other blocks are zero matrices (see more details in [35]).
The characteristic equation of the Jacobian matrixJ is as follows:(13)Fλ=λ2+pλ+q=0,where(14)p=−trJ=−2−r+rxxε−ε−2+re−εx−ax+3ay+2rx,q=detJ=1+r−2rx1−e−εx+rx1−xεe−εx−ay∙1+ax−2ay+a2xy.
In order to obtain a clear insight into the local stability analysis of the fixed pointZx,y, we shall introduce the following definition and lemma from references [8, 36].Figure 1
The positive fixed point is given by the intersection point of the two curves (C1:Fx=r1−x1−e−εx, C2:Fx=ax). (a) No fixed point: rε≤a and (b) unique fixed point: rε>a.
(a)(b)Definition 1.
The fixed pointZx,yis called(i)
sink ifλ1<1 and λ2<1, and it is locally asymptotic stable;(ii)
source ifλ1>1 and λ2>1, and it is locally unstable;(iii)
saddle ifλ1<1 and λ2>1 or vice versa (λ1>1 and λ2<1), and it is locally unstable;(iv)
nonhyperbolic if eitherλ1=1 or λ2=1.Lemma 1.
LetFλ=λ2+pλ+q, suppose thatF1>0, λ1andλ2are two roots of equationFλ=0, then(i)
λ1<1andλ2<1if and only ifF−1>0andq<1;(ii)
λ1>1andλ2>1if and only ifF−1>0andq>1;(iii)
λ1<1andλ2>1(orλ1>1andλ2<1) if and only ifF−1<0;(iv)
λ1=−1andλ2≠1if and only ifF−1=0andq≠0,2;(v)
λ1andλ2are complex conjugates withλ1=λ2=1if and only ifp2−4q<0andq=1.Using Definition1 and Lemma 1, we have the following results.Theorem 2.
Model (4) admits:(1)
The trivial fixed pointZ10,0 is always nonhyperbolic.(2)
The axial fixed pointZ21,0 is(i)
Source if and only ife−ε<1−2/rr>2;(ii)
Saddle if and only ife−ε>1−2/r;(iii)
Nonhyperbolic if and only ife−ε=1−2/rr>2.Proof.
The Jacobian matrixJ of model (4) evaluated at Z10,0 and Z21,0 is given by equation (15):(15)JZ1=I2,JZ2=∂1,12,21−r1−e−ε+∂1,22,2−a+∂2,22,21+a,respectively. Obviously, the trivial fixed point Z10,0 is always nonhyperbolic. It can be easily observed that the eigenvalues of JZ2 are 1+a>1 and 1−r1−e−ε, thus Z21,0 is saddle point if 1−r1−e−ε<1, i.e., e−ε>1−2/r; Z21,0 is source if 1−r1−e−ε>1, i.e., e−ε<1−2/r and 1−2/r>0r>2; Z21,0 is nonhyperbolic if 1−r1−e−ε=1, i.e., e−ε=1−2/r and 1−2/r>0r>2.Theorem 3.
Suppose thatF1>0, then the unique positive fixed pointZ3x∗,y∗is(1)
sink if and only if(16)max−ax∗,ax∗−2J11∗−2<ax∗ax∗−J11∗<1−J11∗,(2)
source if and only if(17)ax∗ax∗−J11∗>max−ax∗,ax∗−2J11∗−2,1−J11∗,(3)
saddle if and only if(18)−ax∗<ax∗ax∗−J11∗<ax∗−2J11∗−2,whereJ11∗=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ax∗.Proof.
We first obtain thatF1>0 holds if ax∗−J11∗>−1. (1)
F−1>0 holds if and only if ax∗ax∗−J11∗>ax∗−2J11∗−2. On the other hand, detJ11∗<1 holds if and only if ax∗ax∗−J11∗<1−J11∗. Hence, Z3x∗,y∗ is sink if and only if(19)max−ax∗,ax∗−2J11∗−2<ax∗ax∗−J11∗<1−J11∗.(2)
The conclusionF−1>0 is the same as the case (1). We obtain detJ11∗>1 holds if and only if ax∗ax∗−J11∗>1−J11∗. Hence, Z3x∗,y∗ is source if and only if(20)ax∗ax∗−J11∗>max−ax∗,ax∗−2J11∗−2,1−J11∗.(3)
We obtainF−1<0 holds if and only if ax∗ax∗−J11∗<ax∗−2J11∗−2. Hence, Z3x∗,y∗ is saddle if and only if(21)−ax∗<ax∗ax∗−J11∗<ax∗−2J11∗−2.It is well known that a sink which is locally asymptotic stable; thus, the condition (16) is also the locally asymptotically stable condition for the positive fixed point Z3x∗,y∗.


### 2.2. Existence and Stability of Fixed Points of Model (5)

Fixed points of model (5) are determined by solving the following nonlinear equations.(22)x=x+rx1−x1−e−εx−axy,y=y+ayx−yyu+y.Using the same method in the previous subsection, the existence result about fixed points of the model (5) in R+2 is summarized as follows.Proposition 1.
Supposer,a,εanduare positive parameter values, model (5) has two boundary fixed pointsZ40,0,Z51,0and unique positive (coexistence) fixed pointZ6x˜∗,y˜∗providing thata<rε, wherex˜∗=y˜∗satisfys(23)r1−x˜∗1−e−εx˜∗−ax˜∗=0.The Jacobian matrixJ˜ of model (5) evaluated at any fixed point x,y is as follows:(24)J˜=∑i,j=12∂i,j2,2J˜ij,where(25)J˜11=1+r1−2x1−e−εx+rx1−xεe−εx−ay,J˜12=−ax,J˜21=ay2u+y,J˜22=1+2ayx−y−ay2u+y−ay2x−yu+y2.The characteristic equation of the Jacobian matrixJ˜ is as follows:(26)F˜λ=λ2+p˜λ+q˜=0,where(27)p˜=−trJ˜=−2−r1−2x1−e−εx−rx1−xεe−εx+ay+ay2−2ayx−yu+y+ay2x−yu+y2,q˜=detJ˜=1+r−2rx1−e−εx+rx1−xεe−εx−ay∙1+2ayx−y−ay2u+y−ay2x−yu+y2+a2xy2u+y.It follows from Lemma1 that the boundary fixed point Z40,0 and Z51,0 are always nonhyperbolic. Using the same method of Theorem 2, some simple manipulation yields the asymptotically stable condition of the positive fixed point Z6x˜∗,y˜∗ in the following theorem.Theorem 4.
The unique positive fixed pointZ6x˜∗,y˜∗is locally asymptotically stable if and only if(28)maxJ˜11∗−ax˜∗2u+x˜∗,−2−J˜11∗+ax˜∗2u+x˜∗<J˜11∗1−ax˜∗2u+x˜∗+a2x˜∗3u+x˜∗<1,where J˜11∗=1+r1−2x∗1−e−εx∗+rx∗1−x∗εe−εx∗−ax∗.


## 2.1. Existence and Stability of Fixed Points of Model (4)

Fixed points of the model (4) are determined by solving the following nonlinear equations:(6)x=x+rx1−x1−e−εx−axy,y=y+ayx−y.For all parameter valuesr,ε, and a, the trivial fixed point Z10,0 and the axial fixed point Z21,0 always exist, they are two boundary fixed points.In the following theorem, we give the existence and the number of the positive (coexistence) fixed point inR+2.Theorem 1.
Assumer,ε, andaare positive real numbers; ifa<rε, the model (4) has a unique positive (coexistence) fixed pointZ3x∗,y∗, wherex∗=y∗and(7)r1−x∗1−e−εx∗−ax∗=0.
Ifa≥rε, the model (4) has no positive fixed point.Proof.
Let(8)Fx=r1−x1−e−εx−ax.
For convenience, we assume(9)F1x=r1−x1−e−εx,F2x=ax.
It is easy to know thatF10=F11=0, and(10)F1′x=−r1−e−εx+r1−xεe−εx,F1″x=−2rεe−εx−r1−xε2e−εx<0,for x∈0,1, so F1x is a concave function in (0, 1). F2x is the straight line crossing Z10,0 with slope a. Thus, there will be only two cases for the roots of Fx=0 as follows:(i)
IfF1′0=rε≤a, then there is no root in (0, 1) of Fx=0, in other words, no positive (coexistence) fixed point exists, as shown in Figure 1(a).(ii)
IfF1′0=rε>a, then there exist exactly one root in (0, 1) of Fx=0, in other words, exactly one positive (coexistence) fixed point Z3x∗,y∗ exists, where x∗=y∗ and r1−x∗1−e−εx∗−ax∗=0, as shown in Figure 1(b).
The Jacobian matrixJ of the model (4) evaluated at any fixed point x,y is as follows:(11)J=∑i,j=12∂i,j2,2Jij,where(12)J11=1+r−2rx1−e−εx+rx1−xεe−εx−ay,J12=−ax,J21=ay,J22=1+ax−2ay.
In equation (12), the symbol∂i,jm,nM denotes a partitioned matrix with m row blocks and n column blocks, and the (i, j)-th block is M, a smaller matrix, all other blocks are zero matrices (see more details in [35]).
The characteristic equation of the Jacobian matrixJ is as follows:(13)Fλ=λ2+pλ+q=0,where(14)p=−trJ=−2−r+rxxε−ε−2+re−εx−ax+3ay+2rx,q=detJ=1+r−2rx1−e−εx+rx1−xεe−εx−ay∙1+ax−2ay+a2xy.
In order to obtain a clear insight into the local stability analysis of the fixed pointZx,y, we shall introduce the following definition and lemma from references [8, 36].Figure 1
The positive fixed point is given by the intersection point of the two curves (C1:Fx=r1−x1−e−εx, C2:Fx=ax). (a) No fixed point: rε≤a and (b) unique fixed point: rε>a.
(a)(b)Definition 1.
The fixed pointZx,yis called(i)
sink ifλ1<1 and λ2<1, and it is locally asymptotic stable;(ii)
source ifλ1>1 and λ2>1, and it is locally unstable;(iii)
saddle ifλ1<1 and λ2>1 or vice versa (λ1>1 and λ2<1), and it is locally unstable;(iv)
nonhyperbolic if eitherλ1=1 or λ2=1.Lemma 1.
LetFλ=λ2+pλ+q, suppose thatF1>0, λ1andλ2are two roots of equationFλ=0, then(i)
λ1<1andλ2<1if and only ifF−1>0andq<1;(ii)
λ1>1andλ2>1if and only ifF−1>0andq>1;(iii)
λ1<1andλ2>1(orλ1>1andλ2<1) if and only ifF−1<0;(iv)
λ1=−1andλ2≠1if and only ifF−1=0andq≠0,2;(v)
λ1andλ2are complex conjugates withλ1=λ2=1if and only ifp2−4q<0andq=1.Using Definition1 and Lemma 1, we have the following results.Theorem 2.
Model (4) admits:(1)
The trivial fixed pointZ10,0 is always nonhyperbolic.(2)
The axial fixed pointZ21,0 is(i)
Source if and only ife−ε<1−2/rr>2;(ii)
Saddle if and only ife−ε>1−2/r;(iii)
Nonhyperbolic if and only ife−ε=1−2/rr>2.Proof.
The Jacobian matrixJ of model (4) evaluated at Z10,0 and Z21,0 is given by equation (15):(15)JZ1=I2,JZ2=∂1,12,21−r1−e−ε+∂1,22,2−a+∂2,22,21+a,respectively. Obviously, the trivial fixed point Z10,0 is always nonhyperbolic. It can be easily observed that the eigenvalues of JZ2 are 1+a>1 and 1−r1−e−ε, thus Z21,0 is saddle point if 1−r1−e−ε<1, i.e., e−ε>1−2/r; Z21,0 is source if 1−r1−e−ε>1, i.e., e−ε<1−2/r and 1−2/r>0r>2; Z21,0 is nonhyperbolic if 1−r1−e−ε=1, i.e., e−ε=1−2/r and 1−2/r>0r>2.Theorem 3.
Suppose thatF1>0, then the unique positive fixed pointZ3x∗,y∗is(1)
sink if and only if(16)max−ax∗,ax∗−2J11∗−2<ax∗ax∗−J11∗<1−J11∗,(2)
source if and only if(17)ax∗ax∗−J11∗>max−ax∗,ax∗−2J11∗−2,1−J11∗,(3)
saddle if and only if(18)−ax∗<ax∗ax∗−J11∗<ax∗−2J11∗−2,whereJ11∗=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ax∗.Proof.
We first obtain thatF1>0 holds if ax∗−J11∗>−1. (1)
F−1>0 holds if and only if ax∗ax∗−J11∗>ax∗−2J11∗−2. On the other hand, detJ11∗<1 holds if and only if ax∗ax∗−J11∗<1−J11∗. Hence, Z3x∗,y∗ is sink if and only if(19)max−ax∗,ax∗−2J11∗−2<ax∗ax∗−J11∗<1−J11∗.(2)
The conclusionF−1>0 is the same as the case (1). We obtain detJ11∗>1 holds if and only if ax∗ax∗−J11∗>1−J11∗. Hence, Z3x∗,y∗ is source if and only if(20)ax∗ax∗−J11∗>max−ax∗,ax∗−2J11∗−2,1−J11∗.(3)
We obtainF−1<0 holds if and only if ax∗ax∗−J11∗<ax∗−2J11∗−2. Hence, Z3x∗,y∗ is saddle if and only if(21)−ax∗<ax∗ax∗−J11∗<ax∗−2J11∗−2.It is well known that a sink which is locally asymptotic stable; thus, the condition (16) is also the locally asymptotically stable condition for the positive fixed point Z3x∗,y∗.


## 2.2. Existence and Stability of Fixed Points of Model (5)

Fixed points of model (5) are determined by solving the following nonlinear equations.(22)x=x+rx1−x1−e−εx−axy,y=y+ayx−yyu+y.Using the same method in the previous subsection, the existence result about fixed points of the model (5) in R+2 is summarized as follows.Proposition 1.
Supposer,a,εanduare positive parameter values, model (5) has two boundary fixed pointsZ40,0,Z51,0and unique positive (coexistence) fixed pointZ6x˜∗,y˜∗providing thata<rε, wherex˜∗=y˜∗satisfys(23)r1−x˜∗1−e−εx˜∗−ax˜∗=0.The Jacobian matrixJ˜ of model (5) evaluated at any fixed point x,y is as follows:(24)J˜=∑i,j=12∂i,j2,2J˜ij,where(25)J˜11=1+r1−2x1−e−εx+rx1−xεe−εx−ay,J˜12=−ax,J˜21=ay2u+y,J˜22=1+2ayx−y−ay2u+y−ay2x−yu+y2.The characteristic equation of the Jacobian matrixJ˜ is as follows:(26)F˜λ=λ2+p˜λ+q˜=0,where(27)p˜=−trJ˜=−2−r1−2x1−e−εx−rx1−xεe−εx+ay+ay2−2ayx−yu+y+ay2x−yu+y2,q˜=detJ˜=1+r−2rx1−e−εx+rx1−xεe−εx−ay∙1+2ayx−y−ay2u+y−ay2x−yu+y2+a2xy2u+y.It follows from Lemma1 that the boundary fixed point Z40,0 and Z51,0 are always nonhyperbolic. Using the same method of Theorem 2, some simple manipulation yields the asymptotically stable condition of the positive fixed point Z6x˜∗,y˜∗ in the following theorem.Theorem 4.
The unique positive fixed pointZ6x˜∗,y˜∗is locally asymptotically stable if and only if(28)maxJ˜11∗−ax˜∗2u+x˜∗,−2−J˜11∗+ax˜∗2u+x˜∗<J˜11∗1−ax˜∗2u+x˜∗+a2x˜∗3u+x˜∗<1,where J˜11∗=1+r1−2x∗1−e−εx∗+rx∗1−x∗εe−εx∗−ax∗.


## 3. Analysis of Bifurcations

In this section, the qualitative behavior of model (4) is analyzed by studying codimension-one bifurcations (flip bifurcation and Neimark–Sacker bifurcation) and codimension-two bifurcation at the unique positive (coexistence) fixed point Z3x∗,y∗. The phenomena of bifurcations are verified by using the center manifold theorem and bifurcation theory [37–39].


### 3.1. Flip Bifurcation

In this subsection, we show that there exist some values of parameters such that model (4) undergoes flip bifurcation by choosing a as a bifurcation parameter. The characteristic equation of model (4) at a positive fixed point Z3x∗,y∗ is as follows:(29)λ2+paλ+qa=0,where(30)pa=re−εx∗1+x∗2ε−x∗ε+2−ax∗+3ay∗+2rx∗−2−r,qa=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗1+ax∗−2ay∗+a2x∗y∗.Let the critical parameter value isa0, if a0 satisfies(31)1+εx∗2−ε+2x∗ra0x∗−2e−εx∗+2+a0x∗−22x∗−1r+2a0x∗−12=0,then, the roots of characteristic (4) at the fixed point Z3x∗a0,y∗a0 are as follows:(32)λ1=−1,λ2=1−pa0.The conditionλ2≠1 implies that pa0≠0,2, i.e.,(33)1+εx∗2−ε+2x∗re−εx∗+2x∗−1r+2a0x∗≠2,4.We define the following set:(34)SFB=a0,r,ε∈R+3:3.2,3.3arefulfilled,a<rε,a,r,ε>0.The variation of the parameters in the small neighborhood of the setSFB gives emergence of flip bifurcation of model (4). Suppose that a0,r,ε∈SFB in the following analysis.Letw=x−x∗,v=y−y∗ and a¯=a−a0. We transform the fixed point Z3x∗,y∗ of model (4) into the origin and expand the Taylor series around w,v,a¯=0,0,0 to the third order, take the a¯ as a new dependent variable, we have the following:(35)wn+1=a1wn+a2vn+a3a¯+a11wn2+a12wnvn+a13wna¯+a23vna¯+a111wn3+a123wnvna¯,vn+1=b1wm+b2vn+b12wnvn+b13wna¯+b22vn2+b23vna¯+b223vn2a¯+b123wnvna¯,where the coefficients of model (35) are given in Appendix A.We construct an invertible matrix(36)T=∂1,12,2a2+∂1,22,2a2+∂2,12,2−a1−1+∂2,22,2−a1+λ2,and use the translation(37)wnvn=TXnYn,then, model (35) can be changed into(38)Xn+1=−Xn+FXn,Yn,a¯,Yn+1=λ2Yn+GXn,Yn,a¯,where(39)FXn,Yn,a¯=1a2λ2+1c3a¯+c11Xn2+c12XnYn+c13Xna¯+c22Yn2+c23Yna¯+c111Xn3+c112Xn2Yn+c113Xn2a¯+c122XnYn2+c123XnYna¯+c222Yn3+c223Yn2a¯+OXn+Yn+a¯4,GXn,Yn,a¯=1a2λ2+1d3a¯+d11Xn2+d12XnYn+d13Xna¯+d22Yn2+d23Yna¯+d111Xn3+d112Xn2Yn+d113Xn2a¯+d122XnYn2+d123XnYna¯+d222Yn3+d223Yn2a¯+OXn+Yn+a¯4,where the coefficients of FXn,Yn,a¯,GXn,Yn,a¯ are given in Appendix A.Now, we use the center manifold theory to determine the center manifoldWc0,0 of model (38) at the fixed point (0, 0) in a small neighborhood of a¯=0. The center manifold Wc0,0 can be presented as follows:(40)Wc0,0=Xn,Yn∈R2Yn=H∗Xn,a¯,H∗0,0=0,DH∗0,0=0.Assume that the center manifold takes the following form:(41)H∗Xn,a¯=ξ0a¯+ξ1Xn2+ξ2Xna¯+ξ3a¯2+OXn+a¯3,then, apply the model (38) on two sides of Yn=H∗Xn,a¯ simultaneously, we have(42)H∗−Xn+FXn,H∗Xn,a¯,a¯,a¯−λ2H∗Xn,a¯−GXn,H∗Xn,a¯,a¯=0.By comparing coefficients of (42), we obtain(43)ξ0=d3a21−λ22,ξ1=d11a21−λ22,ξ2=−2c3ξ1+d13a21+λ22,ξ3=c3a2ξ2λ2+a2ξ2+c3ξ1a221+λ22λ2−1.Therefore, model (38) restricted to the center manifold as follows:(44)ΘXn=−Xn+η1a¯+η2Xn2+η3Xna¯+η4Xn3+η5Xn2a¯+η6Xna¯2+η7a¯2+η8a¯3+OXn+a¯4,where(45)η1=c3a21+λ2,η2=c11a21+λ2,η3=c12ξ0+c13a21+λ2,η4=c12ξ1+c111a21+λ2,η5=c12ξ2+c23ξ1+c112ξ0+2c22ξ0ξ1+c113a21+λ2,η6=c12ξ3+c23ξ2+c123ξ0+c122ξ02+2c22ξ0ξ2a21+λ2,η7=c22ξ02+c23ξ0a21+λ2,η8=2c22ξ0ξ3+c23ξ3+c223ξ02+c222ξ03a21+λ2.In order for map (44) to undergo a flip bifurcation, we require that the following two discriminatory quantities α1 and α2 to be nonzero, where(46)α1=∂2Θ∂Xn∂a¯+12∂Θ∂Xn∂2Θ∂Xn20,0=η3+η1η2≠0,and(47)α2=16∂3Θ∂Xn3+12∂2Θ∂Xn220,0=η4+η22≠0.From the previously mentioned analysis, we obtain the following theorem on flip bifurcation.Theorem 5.
If the bifurcation parameteravaries in the small neighborhood ofa0, and(48)α1≠0,α2≠0,then, model (4) undergoes a flip bifurcation at positive (coexistence) fixed pointZ3x∗,y∗. What’s more, ifα2>0(α2<0) in equation (47), the period-2 orbits that bifurcate from this fixed pointZ3x∗,y∗are stable (respectively unstable).


### 3.2. Neimark–Sacker Bifurcation

Next, we give the existence criteria and the direction for Neimark–Sacker bifurcation by choosinga as a bifurcation parameter. The characteristic equation of the model (4) at positive fixed point Z3x∗a,y∗a is as follows:(49)λ2+paλ+qa=0,and the roots of characteristic (49) have the following form:(50)λ1,2a=−pa±p2a−4qa2,where(51)pa=−J11+1+ax∗−2ay∗,qa=J111+ax∗−2ay∗+a2x∗y∗,J11=1+r1−2x∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗.The eigenvaluesλ1,2a will be a pair of complex conjugate for p2a−4qa<0, that is,(52)22x∗−1r+ax∗−y∗r1+εx∗2−ε+2x∗e−εx∗+r21+εx∗2−ε+2x∗2e−2εx∗+2x∗−12r2+22x∗−1ax∗−y∗r+a2x∗2−6x∗y∗+y∗2<0.Let the critical parameter value bea^, suppose a^ satisfies(53)r1−1a^1−e−ε/a^−1=0.We can obtain thatqa^=1 and λ1,2a^=−pa^/2±i4qa^−p2a^/2=ρ±iυ.We define the following set(54)SNSB=a^,r,ε∈R+3:3.13,3.14aresatisfied,a<rε,a,r,ε>0.The variation of the parameters in the small neighborhood of the setSNSB gives emergence of the Neimark–Sacker bifurcation of model (4). Suppose that a^,r,ε∈SNSB in the following analysis.Under the conditions (52) and (53), we have(55)λ1,2a=qa,(56)d=dλ1,2adaa=a^=12a^3ra^2−ε+2a^+εe−ε/a^−a^r−2a^−2r≠0.In addition, ifpa^≠0,1, which implies(57)1a^2a^2+−ε−2a^+εre−ε/a^−a^r−2a^−2r≠2,3.Correspondingly, we haveλ1,2na^≠1n=1,2,3,4.Ata=a^, we shift the fixed point Z3x∗,y∗ of model (4) into the origin by letting w=x−x∗,v=y−y∗. After the substitution, we obtain(58)wn+1=a1′wn+a2′vn+a11′wn2+a12′wnvn+a111′wn3,vn+1=b1′wm+b2′vn+b12′wnvn+b22′vn2,where a1′,a2′,a11′,a12′,a111′,b1′,b2′,b12′,b22′ are given in Appendix B.Let the invertible matrix(59)T=∂1,12,2a2′+∂1,22,20+∂2,12,2ρ−a1′+∂2,22,2−υ,and under the transformation(60)wnvn=TXnYn,model (58) becomes(61)Xn+1Yn+1=ρ−υυρXn+1Yn+1+FXn,YnGXn,Yn,where(62)FXn,Yn=ζ11Xn2+ζ12XnYn+ζ111Xn3+OXn+Yn4,GXn,Yn=σ11Xn2+σ12XnYn+σ22Yn2+σ111Y3+OXn+Yn4,where the coefficients ζ11,ζ12,ζ111,σ11,σ12,σ22,σ111 are given in Appendix B.In order to guarantee the existence of Neimark–Sacker bifurcation of model (4), the following discriminatory quantity Π will be used to determine the direction of the bifurcation should not be zero:(63)Π=−R1−2λ1λ221−λ1τ1τ2−12τ22−τ32+Rλ2τ4,where(64)τ1=18FXX−FYY+2GXY+iGXX−GYY−2FXY,τ2=14FXX+FYY+iGXX+GYY,τ3=18FXX−FYY−2GXY+iGXX−GYY+2FXY,τ4=116FXXX+FXYY+GXXY+GYYY+iGXXX+GXYY−FXXY−FYYY.In which,(65)FXX=ζ11,FXY=ζ12,FXXX=ζ111,GXX=σ11,GXY=σ12,GYY=σ22,GXXX=σ111,FYY=FXXY=FXYY=FYYY=GXXY=GXYY=GYYY=0.After some calculation, we obtain(66)Π=−132ρ+12−υ2ACE−ADF−BCF−BDE−132FXX2+GXX+GYY2−116−ρFXXX−υGXXX−164FXX−2GXY2+2FXY+GXX−GYY2,where(67)A=2ρ4−2υ4−3ρ3+ρυ2+ρ2−υ2,B=−4υ3ρ+5ρ2υ−4υρ3−2ρυ+υ3,C=FXX+2GXY,D=−2FXY+GXX−GYY,E=FXX,F=GXX+GYY.From the aforementioned discussions, we have the following theorem related to Neimark–Sacker bifurcation.Theorem 6.
If the bifurcation parameteravaries in the small neighborhood ofa^satisfying the conditions (56), (57) and(68)Π≠0,then, model (4) undergoes a Neimark–Sacker bifurcation at positive (coexistence) fixed pointZ3x∗,y∗. More concretely, ifΠ<0(respectively,Π>0), then Neimark–Sacker bifurcation is supercritical (respectively, subcritical) and an attracting (respectively, repelling) invariant closed curve bifurcates from the fixed pointZ3x∗,y∗fora>a^(respectively,a<a^).


### 3.3. Codimension-Two Bifurcation

Now, we discuss the conditions for the occurrence of codimension-two bifurcation associated with 1 : 2 strong resonance of the model (4) by selecting a and r as bifurcation parameters. The characteristic equation of model (4) at a positive fixed point Z3x∗,y∗ is as follows:(69)Fλ=λ2+pa,rλ+qa,r=0,where(70)pa,r=re−εx∗1+x∗2ε−x∗ε+2−ax∗+3ay∗+2rx∗−2−r,qa,r=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗∙1+ax∗−2ay∗+a2x∗y∗.By the relation between roots and coefficients,Fλ has a root −1 with multiplicity 2 if and only if F−1=0 and F0=1. We define the following set:(71)SCT=a,r,ε∈R+3:a=1x∗,r=2εx∗2−εx∗−2x∗+1e−εx∗+2x∗−1,a<rε,a,r,ε>0.It is clear see that ifa,r,ε∈SCT, then the two eigenvalues of JZ3 are λ1,2=−1. Thus, model (4) undergoes codimension-two bifurcation with 1 : 2 strong resonance at Z3x∗,y∗.Taking parametersa1,r1,ε arbitrarily from SCT, we consider model (4) with a=a1+a∗,r=r1+r∗ as follows(72)xn+1=xn+r1+r∗xn1−xn1−e−εxn−a1+a∗xnyn,yn+1=yn+a1+a∗ynxn−yn,where a∗,r∗≪1 are small perturbations of parameters a,r.Letw=x−x∗,v=y−y∗, then we can transform the fixed point Z3x∗,y∗ to the origin and expand the Taylor series at w,v,a∗,r∗=0,0,0,0 to the third-order then, model (72) takes the following form:(73)wn+1=1+m1wn+m2vn+m20wn2+m11wnvn+m30wn3,vn+1=n1wm+1−n2vn+n11wnvn+n02vn2,where the coefficients m1,m2,m20,m11,m30 and n1,n2,n11,n02 are given in Appendix C.Make an invertible matrix(74)T=∂1,12,2m22+m1+∂1,11,2m2−2−m1+∂1,12,1−1+∂2,22,20,and use the transformation(75)wnvn=Tx˜ny˜n,then, model (73) can be rewritten as follows:(76)x˜n+1y˜n+1=−1+θ1a,r1+θ2a,rθ3a,r−1+θ4a,rx˜ny˜n+f˜x˜n,y˜n,a,rg˜x˜n,y˜n,a,r,where(77)f˜x˜n,y˜n,a,r=∑2≤j+k≤3f˜jka,rx˜njy˜nk,g˜x˜n,y˜n,a,r=∑2≤j+k≤3g˜jka,rx˜njy˜nk,with θ1a,r=2−n2−m2n1/2+m1,θ2a,r=n1m2/2+m1−1, θ3a,r=2−n2−m2n1/2+m1, and θ4a,r=2+m1+m2n1/2+m1; the coefficients of f˜x˜n,y˜n,a,r,g˜x˜n,y˜n,a,r are given in Appendix C.We denote the invertible matrix(78)K=∂1,12,21+θ2a,r+∂1,22,20+∂2,12,2−θ1a,r+∂2,22,21,and take the nonsingular linear coordinate transformation(79)x˜ny˜n=Kxnyn,then, model (76) becomes(80)xn+1yn+1=−11θ5a,r−1+θ6a,rx˜ny˜n+f^xn,yn,a,rg^xn,yn,a,r,where θ5a,r=1+θ2a,rθ3a,r−θ1a,rθ4a,r, θ6a,r=θ1a,r+θ4a,r, and(81)f^xn,yn,a,r=11+θ2a,rf˜x˜n,y˜n,a,r=∑2≤j+k≤3f^jka,rxnjynk,g^xn,yn,a,r=θ1a,r1+θ2a,rf˜x˜n,y˜n,a,r+g˜x˜n,y˜n,a,r=∑2≤j+k≤3g^jka,rxnjynk,the coefficients of f^xn,yn,a,r,g^xn,yn,a,r are given in Appendix C.Lastly, we take the following transformation(82)xn=ξn+∑2≤j+k≤3ψjkξnjηnk,yn=ηn+∑2≤j+k≤3ϕjkξnjηnk,where coefficients ψjk and ϕjk2≤j+k≤3 will be confirmed in the following discussion. The inverse transformation of the variable from transformation (82) is(83)ξn=xn+∑2≤j+k≤3ψ¯jkxnjynk+Oxn+yn4,ηn=yn+∑2≤j+k≤3ϕ¯jkxnjynk+Oxn+yn4,in which(84)ψ¯30=ψ30−2ψ202−ψ11ϕ20,ψ¯21=−ψ11ϕ11+ψ21−2ψ02ϕ20−3ψ11ψ20,ψ¯03=ψ03−2ψ02ϕ02−ψ11ψ02,ψ¯12=−2ψ02ϕ11+ψ12−ψ112−ψ11ϕ02−2ψ20ψ02,ϕ¯30=ϕ30−2ϕ20ψ20−ϕ11ϕ20,ϕ¯21=−ϕ11ψ20+ϕ21−2ψ11ϕ20−ϕ112−2ϕ02ϕ20,ϕ¯03=ϕ03−2ϕ022−ψ02ϕ11,ϕ¯12=−ψ11ϕ11+ϕ12−3ϕ11ϕ02−2ψ02ϕ20.Using transformation (82) and (83) for model (80), we can obtain as follows:(85)ξn+1ηn+1=−11θ5a,r−1+θ6a,rξnηn+Ψξn,ηnΦξn,ηn,where(86)Ψξn,ηn=∑2≤j+k≤3Ψjkξnjηnk+Oξn+ηn4,Φξn,ηn=∑2≤j+k≤3Φjkξnjηnk+Oξn+ηn4,and the coefficients of Ψξn,ηn,Φξn,ηn are given in Appendix C.In order to eliminate all quadratic terms, we take(87)Φ20=Φ11=Φ02=Ψ20=Ψ11=Ψ02,and can obtain ψjk and ϕjk for all nonnegative integers j,k and j+k=2. Further, in order to remove all cubic but those resonant terms, we take(88)Φ30=Φ21=Φ12=Φ03=Ψ12=Ψ03,and can obtain ψjk and ϕjk for all nonnegative integers j,k and j+k=3. After the previously mentioned transformations, model (85) can finally be reduced into the following form(89)ξn+1=−ξn+ηn+Oξn+ηn4,ηn+1=θ5a,rξn+−1+θ6a,rηn+Ca,rξn3+Da,rξn2ηn+Oξn+ηn4,where Ca,r=Ψ30 and Da,r=Ψ21. When a,r=a1,r1, we have(90)Ca1,r1=g^30a1,r1+f^20a1,r1g^20a1,r1+12g^202a1,r1+12g^20a1,r1g^11a1,r1,Da1,r1=g^21a1,r1+3f^30a1,r1+12f^20a1,r1g^11a1,r1+54g^20a1,r1g^11a1,r1+g^20a1,r1g^02a1,r1+3f^202a1,r1+52f^20a1,r1g^20a1,r1+52f^11a1,r1g^20a1,r1+g^202a1,r1+12g^112a1,r1.By making use of the results on the existence of 1 : 2 strong resonance bifurcation theorems in [37, 39] and the analysis performed, we have the following theorem.Theorem 7.
If the bifurcation parametersaandrvary in the small neighborhood ofSCT, and(91)Ca1,r1≠0,Da1,r1≠0,then, model (4) undergoes a codimension-two bifurcation with 1 : 2 strong resonance at positive (coexistence) fixed point Z3x∗,y∗.


## 3.1. Flip Bifurcation

In this subsection, we show that there exist some values of parameters such that model (4) undergoes flip bifurcation by choosing a as a bifurcation parameter. The characteristic equation of model (4) at a positive fixed point Z3x∗,y∗ is as follows:(29)λ2+paλ+qa=0,where(30)pa=re−εx∗1+x∗2ε−x∗ε+2−ax∗+3ay∗+2rx∗−2−r,qa=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗1+ax∗−2ay∗+a2x∗y∗.Let the critical parameter value isa0, if a0 satisfies(31)1+εx∗2−ε+2x∗ra0x∗−2e−εx∗+2+a0x∗−22x∗−1r+2a0x∗−12=0,then, the roots of characteristic (4) at the fixed point Z3x∗a0,y∗a0 are as follows:(32)λ1=−1,λ2=1−pa0.The conditionλ2≠1 implies that pa0≠0,2, i.e.,(33)1+εx∗2−ε+2x∗re−εx∗+2x∗−1r+2a0x∗≠2,4.We define the following set:(34)SFB=a0,r,ε∈R+3:3.2,3.3arefulfilled,a<rε,a,r,ε>0.The variation of the parameters in the small neighborhood of the setSFB gives emergence of flip bifurcation of model (4). Suppose that a0,r,ε∈SFB in the following analysis.Letw=x−x∗,v=y−y∗ and a¯=a−a0. We transform the fixed point Z3x∗,y∗ of model (4) into the origin and expand the Taylor series around w,v,a¯=0,0,0 to the third order, take the a¯ as a new dependent variable, we have the following:(35)wn+1=a1wn+a2vn+a3a¯+a11wn2+a12wnvn+a13wna¯+a23vna¯+a111wn3+a123wnvna¯,vn+1=b1wm+b2vn+b12wnvn+b13wna¯+b22vn2+b23vna¯+b223vn2a¯+b123wnvna¯,where the coefficients of model (35) are given in Appendix A.We construct an invertible matrix(36)T=∂1,12,2a2+∂1,22,2a2+∂2,12,2−a1−1+∂2,22,2−a1+λ2,and use the translation(37)wnvn=TXnYn,then, model (35) can be changed into(38)Xn+1=−Xn+FXn,Yn,a¯,Yn+1=λ2Yn+GXn,Yn,a¯,where(39)FXn,Yn,a¯=1a2λ2+1c3a¯+c11Xn2+c12XnYn+c13Xna¯+c22Yn2+c23Yna¯+c111Xn3+c112Xn2Yn+c113Xn2a¯+c122XnYn2+c123XnYna¯+c222Yn3+c223Yn2a¯+OXn+Yn+a¯4,GXn,Yn,a¯=1a2λ2+1d3a¯+d11Xn2+d12XnYn+d13Xna¯+d22Yn2+d23Yna¯+d111Xn3+d112Xn2Yn+d113Xn2a¯+d122XnYn2+d123XnYna¯+d222Yn3+d223Yn2a¯+OXn+Yn+a¯4,where the coefficients of FXn,Yn,a¯,GXn,Yn,a¯ are given in Appendix A.Now, we use the center manifold theory to determine the center manifoldWc0,0 of model (38) at the fixed point (0, 0) in a small neighborhood of a¯=0. The center manifold Wc0,0 can be presented as follows:(40)Wc0,0=Xn,Yn∈R2Yn=H∗Xn,a¯,H∗0,0=0,DH∗0,0=0.Assume that the center manifold takes the following form:(41)H∗Xn,a¯=ξ0a¯+ξ1Xn2+ξ2Xna¯+ξ3a¯2+OXn+a¯3,then, apply the model (38) on two sides of Yn=H∗Xn,a¯ simultaneously, we have(42)H∗−Xn+FXn,H∗Xn,a¯,a¯,a¯−λ2H∗Xn,a¯−GXn,H∗Xn,a¯,a¯=0.By comparing coefficients of (42), we obtain(43)ξ0=d3a21−λ22,ξ1=d11a21−λ22,ξ2=−2c3ξ1+d13a21+λ22,ξ3=c3a2ξ2λ2+a2ξ2+c3ξ1a221+λ22λ2−1.Therefore, model (38) restricted to the center manifold as follows:(44)ΘXn=−Xn+η1a¯+η2Xn2+η3Xna¯+η4Xn3+η5Xn2a¯+η6Xna¯2+η7a¯2+η8a¯3+OXn+a¯4,where(45)η1=c3a21+λ2,η2=c11a21+λ2,η3=c12ξ0+c13a21+λ2,η4=c12ξ1+c111a21+λ2,η5=c12ξ2+c23ξ1+c112ξ0+2c22ξ0ξ1+c113a21+λ2,η6=c12ξ3+c23ξ2+c123ξ0+c122ξ02+2c22ξ0ξ2a21+λ2,η7=c22ξ02+c23ξ0a21+λ2,η8=2c22ξ0ξ3+c23ξ3+c223ξ02+c222ξ03a21+λ2.In order for map (44) to undergo a flip bifurcation, we require that the following two discriminatory quantities α1 and α2 to be nonzero, where(46)α1=∂2Θ∂Xn∂a¯+12∂Θ∂Xn∂2Θ∂Xn20,0=η3+η1η2≠0,and(47)α2=16∂3Θ∂Xn3+12∂2Θ∂Xn220,0=η4+η22≠0.From the previously mentioned analysis, we obtain the following theorem on flip bifurcation.Theorem 5.
If the bifurcation parameteravaries in the small neighborhood ofa0, and(48)α1≠0,α2≠0,then, model (4) undergoes a flip bifurcation at positive (coexistence) fixed pointZ3x∗,y∗. What’s more, ifα2>0(α2<0) in equation (47), the period-2 orbits that bifurcate from this fixed pointZ3x∗,y∗are stable (respectively unstable).


## 3.2. Neimark–Sacker Bifurcation

Next, we give the existence criteria and the direction for Neimark–Sacker bifurcation by choosinga as a bifurcation parameter. The characteristic equation of the model (4) at positive fixed point Z3x∗a,y∗a is as follows:(49)λ2+paλ+qa=0,and the roots of characteristic (49) have the following form:(50)λ1,2a=−pa±p2a−4qa2,where(51)pa=−J11+1+ax∗−2ay∗,qa=J111+ax∗−2ay∗+a2x∗y∗,J11=1+r1−2x∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗.The eigenvaluesλ1,2a will be a pair of complex conjugate for p2a−4qa<0, that is,(52)22x∗−1r+ax∗−y∗r1+εx∗2−ε+2x∗e−εx∗+r21+εx∗2−ε+2x∗2e−2εx∗+2x∗−12r2+22x∗−1ax∗−y∗r+a2x∗2−6x∗y∗+y∗2<0.Let the critical parameter value bea^, suppose a^ satisfies(53)r1−1a^1−e−ε/a^−1=0.We can obtain thatqa^=1 and λ1,2a^=−pa^/2±i4qa^−p2a^/2=ρ±iυ.We define the following set(54)SNSB=a^,r,ε∈R+3:3.13,3.14aresatisfied,a<rε,a,r,ε>0.The variation of the parameters in the small neighborhood of the setSNSB gives emergence of the Neimark–Sacker bifurcation of model (4). Suppose that a^,r,ε∈SNSB in the following analysis.Under the conditions (52) and (53), we have(55)λ1,2a=qa,(56)d=dλ1,2adaa=a^=12a^3ra^2−ε+2a^+εe−ε/a^−a^r−2a^−2r≠0.In addition, ifpa^≠0,1, which implies(57)1a^2a^2+−ε−2a^+εre−ε/a^−a^r−2a^−2r≠2,3.Correspondingly, we haveλ1,2na^≠1n=1,2,3,4.Ata=a^, we shift the fixed point Z3x∗,y∗ of model (4) into the origin by letting w=x−x∗,v=y−y∗. After the substitution, we obtain(58)wn+1=a1′wn+a2′vn+a11′wn2+a12′wnvn+a111′wn3,vn+1=b1′wm+b2′vn+b12′wnvn+b22′vn2,where a1′,a2′,a11′,a12′,a111′,b1′,b2′,b12′,b22′ are given in Appendix B.Let the invertible matrix(59)T=∂1,12,2a2′+∂1,22,20+∂2,12,2ρ−a1′+∂2,22,2−υ,and under the transformation(60)wnvn=TXnYn,model (58) becomes(61)Xn+1Yn+1=ρ−υυρXn+1Yn+1+FXn,YnGXn,Yn,where(62)FXn,Yn=ζ11Xn2+ζ12XnYn+ζ111Xn3+OXn+Yn4,GXn,Yn=σ11Xn2+σ12XnYn+σ22Yn2+σ111Y3+OXn+Yn4,where the coefficients ζ11,ζ12,ζ111,σ11,σ12,σ22,σ111 are given in Appendix B.In order to guarantee the existence of Neimark–Sacker bifurcation of model (4), the following discriminatory quantity Π will be used to determine the direction of the bifurcation should not be zero:(63)Π=−R1−2λ1λ221−λ1τ1τ2−12τ22−τ32+Rλ2τ4,where(64)τ1=18FXX−FYY+2GXY+iGXX−GYY−2FXY,τ2=14FXX+FYY+iGXX+GYY,τ3=18FXX−FYY−2GXY+iGXX−GYY+2FXY,τ4=116FXXX+FXYY+GXXY+GYYY+iGXXX+GXYY−FXXY−FYYY.In which,(65)FXX=ζ11,FXY=ζ12,FXXX=ζ111,GXX=σ11,GXY=σ12,GYY=σ22,GXXX=σ111,FYY=FXXY=FXYY=FYYY=GXXY=GXYY=GYYY=0.After some calculation, we obtain(66)Π=−132ρ+12−υ2ACE−ADF−BCF−BDE−132FXX2+GXX+GYY2−116−ρFXXX−υGXXX−164FXX−2GXY2+2FXY+GXX−GYY2,where(67)A=2ρ4−2υ4−3ρ3+ρυ2+ρ2−υ2,B=−4υ3ρ+5ρ2υ−4υρ3−2ρυ+υ3,C=FXX+2GXY,D=−2FXY+GXX−GYY,E=FXX,F=GXX+GYY.From the aforementioned discussions, we have the following theorem related to Neimark–Sacker bifurcation.Theorem 6.
If the bifurcation parameteravaries in the small neighborhood ofa^satisfying the conditions (56), (57) and(68)Π≠0,then, model (4) undergoes a Neimark–Sacker bifurcation at positive (coexistence) fixed pointZ3x∗,y∗. More concretely, ifΠ<0(respectively,Π>0), then Neimark–Sacker bifurcation is supercritical (respectively, subcritical) and an attracting (respectively, repelling) invariant closed curve bifurcates from the fixed pointZ3x∗,y∗fora>a^(respectively,a<a^).


## 3.3. Codimension-Two Bifurcation

Now, we discuss the conditions for the occurrence of codimension-two bifurcation associated with 1 : 2 strong resonance of the model (4) by selecting a and r as bifurcation parameters. The characteristic equation of model (4) at a positive fixed point Z3x∗,y∗ is as follows:(69)Fλ=λ2+pa,rλ+qa,r=0,where(70)pa,r=re−εx∗1+x∗2ε−x∗ε+2−ax∗+3ay∗+2rx∗−2−r,qa,r=1+r−2rx∗1−e−εx∗+rx∗1−x∗εe−εx∗−ay∗∙1+ax∗−2ay∗+a2x∗y∗.By the relation between roots and coefficients,Fλ has a root −1 with multiplicity 2 if and only if F−1=0 and F0=1. We define the following set:(71)SCT=a,r,ε∈R+3:a=1x∗,r=2εx∗2−εx∗−2x∗+1e−εx∗+2x∗−1,a<rε,a,r,ε>0.It is clear see that ifa,r,ε∈SCT, then the two eigenvalues of JZ3 are λ1,2=−1. Thus, model (4) undergoes codimension-two bifurcation with 1 : 2 strong resonance at Z3x∗,y∗.Taking parametersa1,r1,ε arbitrarily from SCT, we consider model (4) with a=a1+a∗,r=r1+r∗ as follows(72)xn+1=xn+r1+r∗xn1−xn1−e−εxn−a1+a∗xnyn,yn+1=yn+a1+a∗ynxn−yn,where a∗,r∗≪1 are small perturbations of parameters a,r.Letw=x−x∗,v=y−y∗, then we can transform the fixed point Z3x∗,y∗ to the origin and expand the Taylor series at w,v,a∗,r∗=0,0,0,0 to the third-order then, model (72) takes the following form:(73)wn+1=1+m1wn+m2vn+m20wn2+m11wnvn+m30wn3,vn+1=n1wm+1−n2vn+n11wnvn+n02vn2,where the coefficients m1,m2,m20,m11,m30 and n1,n2,n11,n02 are given in Appendix C.Make an invertible matrix(74)T=∂1,12,2m22+m1+∂1,11,2m2−2−m1+∂1,12,1−1+∂2,22,20,and use the transformation(75)wnvn=Tx˜ny˜n,then, model (73) can be rewritten as follows:(76)x˜n+1y˜n+1=−1+θ1a,r1+θ2a,rθ3a,r−1+θ4a,rx˜ny˜n+f˜x˜n,y˜n,a,rg˜x˜n,y˜n,a,r,where(77)f˜x˜n,y˜n,a,r=∑2≤j+k≤3f˜jka,rx˜njy˜nk,g˜x˜n,y˜n,a,r=∑2≤j+k≤3g˜jka,rx˜njy˜nk,with θ1a,r=2−n2−m2n1/2+m1,θ2a,r=n1m2/2+m1−1, θ3a,r=2−n2−m2n1/2+m1, and θ4a,r=2+m1+m2n1/2+m1; the coefficients of f˜x˜n,y˜n,a,r,g˜x˜n,y˜n,a,r are given in Appendix C.We denote the invertible matrix(78)K=∂1,12,21+θ2a,r+∂1,22,20+∂2,12,2−θ1a,r+∂2,22,21,and take the nonsingular linear coordinate transformation(79)x˜ny˜n=Kxnyn,then, model (76) becomes(80)xn+1yn+1=−11θ5a,r−1+θ6a,rx˜ny˜n+f^xn,yn,a,rg^xn,yn,a,r,where θ5a,r=1+θ2a,rθ3a,r−θ1a,rθ4a,r, θ6a,r=θ1a,r+θ4a,r, and(81)f^xn,yn,a,r=11+θ2a,rf˜x˜n,y˜n,a,r=∑2≤j+k≤3f^jka,rxnjynk,g^xn,yn,a,r=θ1a,r1+θ2a,rf˜x˜n,y˜n,a,r+g˜x˜n,y˜n,a,r=∑2≤j+k≤3g^jka,rxnjynk,the coefficients of f^xn,yn,a,r,g^xn,yn,a,r are given in Appendix C.Lastly, we take the following transformation(82)xn=ξn+∑2≤j+k≤3ψjkξnjηnk,yn=ηn+∑2≤j+k≤3ϕjkξnjηnk,where coefficients ψjk and ϕjk2≤j+k≤3 will be confirmed in the following discussion. The inverse transformation of the variable from transformation (82) is(83)ξn=xn+∑2≤j+k≤3ψ¯jkxnjynk+Oxn+yn4,ηn=yn+∑2≤j+k≤3ϕ¯jkxnjynk+Oxn+yn4,in which(84)ψ¯30=ψ30−2ψ202−ψ11ϕ20,ψ¯21=−ψ11ϕ11+ψ21−2ψ02ϕ20−3ψ11ψ20,ψ¯03=ψ03−2ψ02ϕ02−ψ11ψ02,ψ¯12=−2ψ02ϕ11+ψ12−ψ112−ψ11ϕ02−2ψ20ψ02,ϕ¯30=ϕ30−2ϕ20ψ20−ϕ11ϕ20,ϕ¯21=−ϕ11ψ20+ϕ21−2ψ11ϕ20−ϕ112−2ϕ02ϕ20,ϕ¯03=ϕ03−2ϕ022−ψ02ϕ11,ϕ¯12=−ψ11ϕ11+ϕ12−3ϕ11ϕ02−2ψ02ϕ20.Using transformation (82) and (83) for model (80), we can obtain as follows:(85)ξn+1ηn+1=−11θ5a,r−1+θ6a,rξnηn+Ψξn,ηnΦξn,ηn,where(86)Ψξn,ηn=∑2≤j+k≤3Ψjkξnjηnk+Oξn+ηn4,Φξn,ηn=∑2≤j+k≤3Φjkξnjηnk+Oξn+ηn4,and the coefficients of Ψξn,ηn,Φξn,ηn are given in Appendix C.In order to eliminate all quadratic terms, we take(87)Φ20=Φ11=Φ02=Ψ20=Ψ11=Ψ02,and can obtain ψjk and ϕjk for all nonnegative integers j,k and j+k=2. Further, in order to remove all cubic but those resonant terms, we take(88)Φ30=Φ21=Φ12=Φ03=Ψ12=Ψ03,and can obtain ψjk and ϕjk for all nonnegative integers j,k and j+k=3. After the previously mentioned transformations, model (85) can finally be reduced into the following form(89)ξn+1=−ξn+ηn+Oξn+ηn4,ηn+1=θ5a,rξn+−1+θ6a,rηn+Ca,rξn3+Da,rξn2ηn+Oξn+ηn4,where Ca,r=Ψ30 and Da,r=Ψ21. When a,r=a1,r1, we have(90)Ca1,r1=g^30a1,r1+f^20a1,r1g^20a1,r1+12g^202a1,r1+12g^20a1,r1g^11a1,r1,Da1,r1=g^21a1,r1+3f^30a1,r1+12f^20a1,r1g^11a1,r1+54g^20a1,r1g^11a1,r1+g^20a1,r1g^02a1,r1+3f^202a1,r1+52f^20a1,r1g^20a1,r1+52f^11a1,r1g^20a1,r1+g^202a1,r1+12g^112a1,r1.By making use of the results on the existence of 1 : 2 strong resonance bifurcation theorems in [37, 39] and the analysis performed, we have the following theorem.Theorem 7.
If the bifurcation parametersaandrvary in the small neighborhood ofSCT, and(91)Ca1,r1≠0,Da1,r1≠0,then, model (4) undergoes a codimension-two bifurcation with 1 : 2 strong resonance at positive (coexistence) fixed point Z3x∗,y∗.


## 4. Analysis of Marotto’s Chaos

In the present section, Marotto’s theorem in references [40, 41] is used to prove that the unique positive (coexistence) fixed point Z3x∗,y∗ of model (4) is a snap-back repeller and model (4) possesses a chaotic behavior in the sense of Marotto.Before giving the conditions for the existence of chaotic phenomena, we first present the definitions and theorem of Marotto’s chaos given in references [32, 33]. For any map F:Rn⟶Rn, and any positive integer k, let Fk represent the composition of F with itself k times. For a differentiable function F, let DFZ denote the Jacobian matrix of F evaluated at the point Z∈Rn, and DFZ its determinant. Let UrZ denote the closed ball in Rn of radius r centered at the point Z and Ur0Z its interior. Also, let Z be the usual Euclidean norm of Z in Rn.Definition 2.
LetFbe differentiable inUrZ0. The pointZ0∈Rnis an expanding fixed point ofFinUrZ0, ifFZ0=Z0and all eigenvalues ofDFZexceed 1 in norm for allZ∈UrZ0.Definition 3.
Suppose thatZ0is an expanding fixed point of a given mapFinUrZ0withr>0, thenZ0is called a snap-back repeller ofFif there exists a pointZ∗∈UrZ0satisfiesZ∗≠Z0, FMZ∗=Z0, andDFMZ∗≠0for some positive integerM.Theorem 8.
IfFhas a snap-back repeller, then the mapFis chaotic.The characteristic equationFλ of the Jacobian matrix J of model (4) has a pair of conjugate complex roots λ1 and λ2, and λ1=λ2>1 if and only if(92)p2x,y−4qx,y<0,qx,y>1,where px,y=−trJ,qx,y=detJ.In order to give the conditions under which the fixed pointZ3x∗,y∗ is a snap-back repeller, a neighborhood UrZ3 of Z3x∗,y∗ is found to satisfy (92) for all Z∈UrZ3. For convenience, the following equations are proposed:(93)S1x,y=p2x,y−4qx,y,S2x,y=qx,y−1.By some calculations, we have(94)S1x,y=a2y2−2aβ1+2axy+β12,and(95)∆1=2aβ1+2ax2−4a2β12,where(96)β1=rxx−1εe−εx+r1−2xe−εx+ax+2rx−r.Then, the equationS1x,y=0 has two real roots(97)y¯1=β1+2ax−4a∆11a,y¯2=β1+2ax+4a∆11a,only if ∆11=xβ1+ax>0. Since ∆11x=0=0,∆11x=1=2a+r1−e−ε>0, we denote D1=x∆11>0. Thus, S1<0 if x∈D1 and y∈D3=y¯1,y¯2.On the other hand,(98)S2x,y=2a2y2+a2β2−1y+r−2rx−aβ2x+β3,with(99)∆2=a22β2−12−8a2r−2rx−aβ2x+β3,where(100)β2=rεx2−εx−2x+1e−εx+2rx−r−1,β3=−e−εxr1+εx2−εx−2x.Then, the equationS2x,y=0 has two real roots(101)y¯1′=1−2β2−∆224a,y¯2′=1−2β2+∆224a,when(102)∆22=Q2+2Q22x−1r+2ax−1−2ax+12+22x−1r+2ax−12+9≥0,where Q=21+εx2−ε+2xre−εx. Denote D2=x∆22≥0,D2′=x∆22<0, thus, S2x,y>0 if one of the following conditions is satisfied.(1)
x∈D2 and y∈D3′=−∞,y¯1′∪y¯2′,+∞;(2)
x∈D2′ and y∈R.According to the previous analysis, we obtain the following lemma.Lemma 2.
Leta>0,r>0, andε>0, denote the following conditions.(103)UP1:x∈D1∩D2andy∈D3∩D3′;UP2:x∈D1∩D2′andy∈D3.If the fixed pointZ3x∗,y∗ of model (4) satisfies(104)Z3x∗,y∗∈UrZ3=UP1,UP2,then, Z3x∗,y∗ is called an expanding fixed point of F¯ in UrZ3, F¯ is defined by the map (4).By Definition3 of a snap-back repeller, we need to find one point Z^x1,y1∈UrZ3 which satisfies the conditions Z^≠Z3,F¯MZ^=Z3 and DF¯MZ^≠0, where F¯M is the M times iterations, M is an integer and M>1.Note that(105)x1+rx11−x11−e−εx1−ax1y1=x2,y1+ay1x1−y1=y2,and(106)x2+rx21−x21−e−εx2−ax2y2=x∗,y2+ay2x2−y2=y∗.If equations (105) and (106) have solutions which are all different from Z3, then a map F¯2 is constructed to map the point Z^x1,y1 to the fixed point Z3x∗,y∗ after two iterations.The solutionsx2,y2 of (104) satisfy the following equations:(107)y2=1ax2x2−x∗+rx21−x21−e−εx2,ay∗x2=1x2−rx2x2−1e−εx2+a+rx22−rx2+x∗rx2x2−1e−εx2−rx22+r+1x2−x∗.Puttingx2 and y2 in (105), and solving x1,y1, we obtain(108)y1=1ax1x1−x2+rx11−x11−e−εx1,ay2x1=1x1−rx1x1−1e−εx1+a+rx12−rx1+x2rx1x1−1e−εx1−rx12+r+1x1−x2.By a simple calculation, we obtain(109)DF¯2Z^=−r1+x22ε−ε+2x2φ1e−ε2−φ12x2−1r+φ2−r1+x12ε−ε+2x1ψ1e−εx1−ψ12x1−1r+ψ2,where(110)φ1=1+x2−2y2a,φ2=1+2a2y22+x2−3y2a,ψ1=1+x1−2y1a,ψ2=1+2a2y12+x1−3y1a.According to the previous analysis, we have the following theorem.Theorem 9.
Assume that the conditions in Lemma 2hold. If the solutions of equations (107) and (108) satisfy that(111)xi,yi≠x∗,y∗,i=1,2,x1,y1∈UZ3,x1,y1≠0,0andDF¯2Z^≠0,thenZ3x∗,y∗is a snap-back repeller of model (4) inUZ3, and in this case the model (4) is chaotic in the sense of Marotto.Next, we give an example to illustrate the existence of conditions in Theorem9.Example 1.
Fora=2.1,r=2.55, and ε=4.5, the map (4) has a positive fixed point Z3=0.5236307748,0.5236307748 and the eigenvalues associated with Z3 is λ1,2=−0.1854542600±1.096631426i. According to the conditions of Lemma 2 and Theorem 9, we find that a region of Z3 (see Figure 2) is(112)U=x,y0.4767978417<x<0.5652456196,y¯1<y<y¯2⊂UP2,U=x,y0.3140639424<x<0.4767978417,y∈D3∩D3′⊂UP1,U=x,yx>0.5652456196,y∈D3∩D3′⊂UP1.and there exists a point Z^=0.5433963492,0.5396159152 satisfying that F¯2Z^=Z3 and(113)DF¯2Z^=1.544079215≠0,where(114)y¯1,2=5.42857143x+e−4.5x5.464285715x2−7.89285714x+1.214285714∓0.755928946e−4.5x38.25x3−55.25x2+8.5x+x31x−8.5−1.214285714,y¯1,2′=−1.214285714x+e−4.5x−2.732142858x2+3.94642857x−0.607142857∓0.566893424K1e−9x+K2e−4.5x+K3+0.9642857142,where K1=2322.75802x4−6710.18985x3+5878.585125x2−1491.1533x+114.7041, K2=2914.8336x3−5445.0711x2+2431.2771x−274.3902, and K3=199.3761+x836.6652x−811.7928. We can easily observe that Z3,Z^∈U. Therefore, Z3 is a snap-back repeller and the model (4) is chaotic in the sense of Marotto.Figure 2
The different regions ofZ3 with a=2.1,r=2.55, and ε=4.5.


## 5. Numerical Simulations

Diagrams for bifurcation, phase-plane, maximum Lyapunov exponents, the basin of attractor, two-parameter dynamical map, and time sequence of the model (4) and model (5) are drawn to validate the theoretical results in this section. Furthermore, we find some new, complex, and surprising dynamical behaviors as the parameters vary.


### 5.1. Numerical Simulations for Bifurcations and Marotto’s Chaos atZ3

#### 5.1.1. Flip Bifurcation

The case of flip bifurcation is given in Figure3. Suppose that parameters 0.1≤a≤1.1, r=2.9 and ε=4.5, based on Theorem 1, we know model (4) has only one positive fixed point Z3. Figure 3(a) shows the bifurcation diagram of the model with respect to the bifurcating parameter a in a,x plane with initial values (0.9, 0.9), the flip bifurcation (labeled “FB”) takes place for the fixed point Z3=0.8036104706,0.8036104706 with a=0.6896609431,α1=1.664046302 and α2=2.266430550>0. The chaotic attractor for a=0.2 is given in Figure 3(b) and the corresponding time sequence diagram is given in Figure 3(c) which shows the chaotic nature. Figure 3 confirms the correctness of the theoretical results in Theorem 5.Figure 3
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=0.2 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


#### 5.1.2. Neimark–Sacker Bifurcation

The case of Neimark–Sacker bifurcation is given in Figure4. Suppose that parameters 0.8≤a≤1.8,r=3 and ε=3.4, based on Theorem 1, we know model (4) has only one positive fixed point. Figure 4(a) shows that the bifurcation diagram of model (4) with respect to the bifurcating parameter a in a,x plane with initial values (0.5, 0.5), the Neimark–Sacker bifurcation (labeled “NSB”) takes place for fixed point Z3=0.6206980308,0.6206980308 and a=1.611089371 with d=0.7278891375>0 and Π=−1.722080647<0, an attracting invariant cycles appears. The chaotic attractor for a=1.778 is given in Figure 4(b) and the corresponding time sequence diagram is given in Figure 4(c) which shows the chaotic nature. Figure 4 confirms the correctness of the theoretical results in Theorem 6.Figure 4
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=1.778 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


#### 5.1.3. Codimension-Two Bifurcation

The case of codimension-two bifurcation is given in Figure5. Suppose that parameters a=1.307020867, r=4.733967524 and ε=3, based on Theorem 1, we know the model (4) has only one positive fixed point. After calculating, we obtain the positive fixed point Z3=0.7592142764,0.7592142764, and the eigenvalues of the Jacobian matrix JZ3 are λ1,2=−1. Furthermore, we have a=a1=1.307020867 and r=r1=4.733967524, thus a,r,ε∈CT, Ca1,r1=12.01168114≠0 and Da1,r1=107.5939727≠0. Figure 5(a) shows that the three-dimensional bifurcation of the model (4) in a,r,x plane for ε=3 and a,r vary in a neighborhood of a1,r1=1.307020867,4.733967524 with initial values (0.6, 0.6). Therefore, it follows from Theorem 7 that a fixed point Z3 is a 1 : 2 strong resonance bifurcation point. Figure 5(b) shows the two-dimensional bifurcation diagram in a,x plane when r=4.733967524. Figure 5 confirms the correctness of the theoretical results in 3.3.Figure 5
(a) Bifurcation diagram of model (4) in a,r,x plane. (b) Bifurcation diagram of model (4) in a,x plane with r=4.733967524.
(a)(b)


#### 5.1.4. Marotto’s Chaos

Using Example1, the case of Marotto’s chaos is given in Figure 6. Suppose that parameters 2.54≤r≤2.57, a=2.1 and ε=4.5, based on Theorem 1, we know the model (4) has only one positive fixed point Z3. Figure 6(a) shows the bifurcation diagram in r,x plane with initial values (0.3, 0.3). The maximum Lyapunov exponents corresponding to (a) are calculated and plotted in Figure 6(b). It can be seen from Figure 6(b) that all the maximal Lyapunov exponents are larger than 0, so the model (4) with the given initial values are all chaotic. From Example 1 and Theorem 9, one can know that the chaotic attractor is located in the chaotic region r∈2.54,2.57, which is in line with the maximal Lyapunov exponents in Figure 6(b). The Marotto’s chaotic attractor is exhibited in Figure 6(c) for r=2.55. Figure 6 verifies the correctness of the theoretical results in Theorem 9.Figure 6
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Maximum Lyapunov exponents corresponding to (a). (c) Chaotic attractor for r=2.55 in (a).
(a)(b)(c)


### 5.2. Numerical Simulations for Model (4) in Different Bifurcation Parameters

In this part, we investigate the influence of the parameters on the dynamics of the system (4), and some new and more complex dynamic behaviors are presented as the parameters vary. Bifurcation parameters are considered in the following four cases, respectively:(1)
Varya in range 0≤a≤1.9, and fix r=3,ε=4;(2)
Varya in range 0≤a≤1.9, and fix r=3,ε=5.61;(3)
Varyε in range 1≤ε≤3.4, and fix a=1.8,r=3;(4)
Varyr in range 4≤r≤6.5, and fix a=1.8,ε=1.


#### 5.2.1. Case (1)

Taking parametersr=3,ε=4 and a covering [0, 1.9] with initial values (0.66, 0.66), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 7(a) and the local amplification is given in Figure 7(b) for a∈0.3,0.8. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 7(c) which confirms the existence of the chaos and period orbits in the parametric space. In Figure 7(a), there is a stable fixed point for more large regions a∈0.714,1.57, the stable fixed point loses its stability at a∼1.57 and closed invariant circle appears when the parameter a exceeds 1.57 and eventually leads to chaos at a∼1.845. There are two different routes (invariant circle and inverse period-doubling bifurcation) to the chaos which occurs simultaneously. The inverse period-doubling bifurcation phenomena occurs at a∼0.714, and there are some cascades of period-doubling bifurcations lead to chaos. The chaotic behavior suddenly appears at a∼0. Figure 7(b) shows that there is an unending sequence of bifurcations until chaos. The phase-plane diagrams corresponding to Figure 7(a) are shown in Figures 7(d)–7(f) for the two-coexisting chaotic sets at a=0.315, four-coexisting chaotic sets at a=0.353, and the invariant circle at a=1.57799, respectively.Figure 7
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈0.3,1.2 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 6(a): (d) a=0.315, (e) a=0.353, and (f) a=1.57799.
(a)(b)(c)(d)(e)(f)


#### 5.2.2. Case (2)

Taking parametersr=3,ε=5.61 and a covering [0, 1.9] with initial values (0.08, 0.08), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 8(a) and the local amplifications are plotted in Figures 8(b) and 8(c) for a∈1.25,1.4 and a∈1.4,1.5, respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 8(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 8(a) shows that there is a stable fixed point for a∈0.8075,1.173 and the fixed point loses its stability as a increases. The period-3 orbits are observed for ∈1.173,1.297, and there happen two different routes to chaos at the same time; they are inverse period-doubling bifurcation and period-doubling bifurcation. When a increases at a=1.36, a period-6 orbit appears, and some cascades of period-doubling bifurcations lead to chaos. The chaotic behaviors suddenly appear at a∼1.65. One can also observe that there is an inverse period-doubling bifurcation which occurs at a∼0.8075 changes into chaos with complex period windows and the chaotic behavior suddenly appears at a∼0. Figure 8(c) depicts the alternates of model (4) between periodic and chaotic motions. The phase-plane diagrams corresponding to Figure 8(a) are shown in Figures 8(e)–8(i) for the two-coexisting chaotic sets at a=0.42, period-3 orbits at a=1.27, period-6 orbits at a=1.36, three-coexisting chaotic sets at a=1.571, and chaotic attractor is at a=1.65, respectively.Figure 8
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈1.25,1.4 corresponding to (a). (c) Local amplification for a∈1.4,1.5 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 7(a): (e) a=0.42, (f) a=1.27, (g) a=1.36, (h) a=1.571, and (i) a=1.65.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


#### 5.2.3. Case (3)

Taking parametersa=1.8,r=3 and ε covering [1, 3.4] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in ε,x plane is given in Figure 9(a). Figure 9(b) is the local amplifications for 2.93<ε<3.1. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 9(c) which confirms the existence of the chaos and period orbits in the parametric space. Figure 9(a) shows that the stability of fixed point happens for ε∈1,2.495, loses its stability as ε increases and through a Neimark–Sacker bifurcation for ε∼2.495. For ε>2.495, the fixed point becomes unstable and a closed invariant circle appears, the radius of the circle becomes bigger and bigger, and eventually leads to chaos. In particular, for ε>2.88, the invariant circle disappears, then it becomes period-4 orbits at ε∼2.976 suddenly. Model (4) undergoes period-doubling phenomena, which leads to chaos with complex period windows. The phase-plane diagrams for various values of ε are shown in Figures 9(d)–9(l), which clearly illustrate the act of how an invariant circle bifurcates from the stable fixed point and increases its radius, the invariant circle leads to chaotic attractors and the chaotic behavior suddenly disappears at ε∼3.41. Figure 9 shows that there are period-4, period-8, period-30, four-coexisting chaotic sets, and chaotic attractor.Figure 9
(a) Bifurcation diagram of model (4) in ε,x plane. (b) Local amplification for ε∈2.93,3.1 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of ε corresponding to Figure 8(a): (d) ε=1.9, (e) ε=2.5, (f) ε=2.875, (g) ε=2.88, (h) ε=2.976, (i) ε=3.048, (j) ε=3.101, (k) ε=3.131, and (l) ε=3.321.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


#### 5.2.4. Case (4)

Taking parametersa=1.8,ε=1 and r covering [4, 6.5] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in r,x plane is given in Figure 10(a) and the local amplifications for r∈5.82,5.96 and r∈5.96,6.3 are given in Figures 10(b) and 10(c), respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 10(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 10(a) shows that the stability areas of a fixed point is r∈4,5.279, the fixed point becomes unstable as r increases. This would suggest that a population increases the carrying capacity at a certain growth rate and then lose its stability. The phase-plane diagrams in Figures 10(e)–10(m) clearly depicts the process of an invariant curve becoming increasingly complex as parameter r growth and chaotic attractors appear. When r increases at certain values, for example, at r>5.279, a closed curve appears, indicating that model (4) undergoes a Neimark–Sacker bifurcation, and the quasiperiodic orbit arise. At r=5.85, the circle disappears and a period-9 orbit appears, then evolves into chaos. Figure 10 shows the rich types of dynamics for the model (4); there are period-9, period-43, quasiperiodic orbit, nonattracting chaotic set, and chaotic attractor.Figure 10
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Local amplification for r∈5.82,5.96 corresponding to (a). (c) Local amplification for r∈5.96,6.3 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of r corresponding to Figure 9(a): (e) r=4.9, (f) r=5.291, (g) r=5.5711, (h) r=5.582, (i) r=5.59, (j) r=5.85, (k) r=6.055, (l) r=6.08, and (m) r=6.16.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)


### 5.3. Coexisting Attractors and Multistability of Model (4)

As mentioned in Subsection5.2, model (4) can exhibit a variety of rich dynamical behaviors when varying the parameters a, r and ε, respectively. In this subsection, we will present some coexisting attractors and multistability to obtain more insightful information on the complex dynamics of the model (4).With a view of exploring more detailed dynamics of model (4), we first plot a two-parameter dynamical diagram [42] as illustrated in Figure 11 to present the impact of the parameters a, r, and ε. It can be clearly seen from Figure 11 that there are two color areas in which the yellow area denotes the periodic state and the dark-blue region represents the chaotic state. Figures 11(a) and 11(b) show the dynamics regions of the model (4) in the a,r plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Figures 11(c) and 11(d) show the dynamics regions of the model (4) in the a,ε plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Observed from Figures 11(a) and 11(b), there are some finely banded yellow regions occurring in the dark-blue area, indicating that the orbits of the model (4) alter many times in these regions between period and chaos. This also shows that model (4) is extremely sensitive in these regions, meaning the emergence of the complicated state transition behavior. Furthermore, when we vary both the parameters a and r within a certain range, one can observe that the color area is always relatively concentrated and complete in the yellow region (chaos) or dark-blue region (period), which means model (4) has strong structural stability. Especially, viewing the color of the dynamical map horizontally along the direction of r=3, the dynamical behaviors are in good agreement with the bifurcation diagram of Figure 12(a). Figures 11(c) and 11(d) show that when the parameter a is small, model (4) is in a chaotic state. Then it converts into the periodic state, next drop into chaos again, and finally go out of chaos into periodic orbits with a gradually increasing. We also see that period is mainly concentrated in the left half of the graph while chaos has mainly appeared in the right half. Remark that seeing the color of the dynamical map horizontally along the direction of ε=4.8, the dynamical behaviors also coincide well with the bifurcation diagram of Figure 12(a). A comparison of Figures 11(a) and 11(b) (or Figures 11(c) and 11(d)) shows that different initial conditions may be able to have different dynamical behaviors when taking the same parameters a and r. In conclusion, the two-parameter dynamical map, as a useful tool, can offer a better parameter selection for population evolution in the real world.Figure 11
(a, b) Dynamical map with respect to the parametersa and r with ε=4.8 and initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. (c, d) Dynamical map with respect to the parameters a and ε with r=3 and initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively.
(a)(b)(c)(d)Figure 12
Coexisting dynamics with respect to the parametera covering [0,1] with initial conditions (0.08, 0.08) (cyan) and initial conditions (0.7, 0.7) (blue). (a) Bifurcation diagram of the model (4) in a,x plane. (b) Maximum Lyapunov exponents corresponding to (a).
(a)(b)Based on the previously discussed analysis, our attention is further focused on the study of coexisting attractors and multistability. The coexisting attractors is an important nonlinear dynamics which means two or more attractors generate simultaneously under different initial conditions but with the same parameter set [43]. Particularly, owing to only one positive stable fixed point, the resulting attractors are all hidden [43, 44]. The trajectories of the model with multistability can asymptotically approach different states when varying initial conditions and fixing parametric values. In order to observe this interesting phenomenon in the model (4), suppose that r=3,ε=4.8, and select a as the adjustable parameter. The coexisting bifurcation diagram with varying parameter 0≤a≤1 is displayed in Figure 12(a) and the corresponding maximum Lyapunov exponents to (a) are plotted in Figure 12(b). Note that the cyan trajectory and blue trajectory in Figure 12 begin with the initial conditions (0.08, 0.08) and (0.7, 0.7), respectively. The parametric region in which Figure 12(a) presents different behaviors is called the coexisting zone. Obviously, from Figure 12(a), one can observe the coexistence of period and chaos in an extremely narrow area between the two vertical red dashed lines with 0.191≤a≤0.209. With the parameter varying in the coexisting zone, there are chaotic attractors coexist with the periodic attractors (e.g., period-3 orbits, period-6 orbits, and period-12 orbits). When setting the parameter a=0.206 and keeping other parameters unchanged, the phase portrait of the coexistence of the period-3 attractor and chaotic attractor are illustrated in Figure 13. It is worth noting that the bifurcation diagram matches well with its corresponding maximum Lyapunov exponents.Figure 13
Phase portraits for period-3 and chaos orbits.It is well known that the coexisting attractors can be characterized by the basins of attraction, and the multistability depends on the different basins of attraction [42]. An attractor’s attraction basin is the region of the initial condition plane, which tend to the attractor as the time approaches infinity. That is when we start with points that are within the basin of attraction, they will finally be iterated into the attractor [45, 46]. Fix the parameters a=0.206,r=3 and ε=4.8. The basin of attraction of the coexisting period-3 attractor and chaotic attractor of Figure 13 is shown in Figure 14. Different colors represent different attraction regions, where the orange-red and the royal blue indicate period-3 attractor and chaotic attractor, respectively. Figure 14 depicts such model is crucially dependent on the initial states, and selecting the initial value in different colors, the model (4) will generate chaotic or periodic attractors. What is more, it’s easy to observe that the model is delicate with respect to the initial conditions. Especially, the orange-red region and royal blue region are interweaved with each other in the bottom half of the graph, implying that a tiny perturbation to the initial state is likely to lead to a state transition behavior.Figure 14
The basins of attraction in thex−y initial plane.


### 5.4. Impact of the Allee Effect on the Stability of Model (5)

For better visualization of the impact of the Allee parameteru on the stability of model (5), we draw the bifurcation diagrams of model (5) in a,x plane, ε,x plane and r,x plane, respectively. Bifurcation parameters are considered in four cases which is the same as subsection 5.2.In Figure15, the four different values of uu=0.02,0.05,0.08,0.11 are selected. Figures 15(a)–15(d) show that if we gradually increase the value of u for the model (5) with the Allee effect on both populations, then model (5) becomes stable, i.e., the Allee effect in predators has a stabilizing effect on the model (5) with the set of parameter values, the chaotic dynamics is observed at large values of a than the model without Allee constant u (see Figures 7(a) and 8(a)). However, changing the values of u has no significant effect on the inverse period-doubling bifurcation (see Figures 15(a)–15(h)). In fact, we can use the feedback control technique to change the chaotic dynamics of inverse period-doubling bifurcation (see Section 6). A similar behavior is observed in the remaining subfigures ((d)-(h)), ((i)-(l)), and ((m)-(p)) also. The results we obtained mean that the Allee constant u plays a very important role in managing the population.Figure 15
Bifurcation diagrams of model (5) with different values of u. The parameter values for a−d corresponding to case (1). The parameter values for e−h corresponding to case (2). The parameter values for i−l corresponding to case (3). The parameter values for m−p corresponding to case (4). (a) u = 0.02, (b) u = 0.05, (c) u = 0.08, (d) u = 0.011, (e) u = 0.02, (f) u = 0.05, (g) u = 0.08, (h) u = 0.011, (i) u = 0.02, (j) u = 0.05, (k) u = 0.08, (l) u = 0.011, (m) u = 0.02, (n) u = 0.05, (o) u = 0.08, and (p) u = 0.011.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)


## 5.1. Numerical Simulations for Bifurcations and Marotto’s Chaos atZ3

### 5.1.1. Flip Bifurcation

The case of flip bifurcation is given in Figure3. Suppose that parameters 0.1≤a≤1.1, r=2.9 and ε=4.5, based on Theorem 1, we know model (4) has only one positive fixed point Z3. Figure 3(a) shows the bifurcation diagram of the model with respect to the bifurcating parameter a in a,x plane with initial values (0.9, 0.9), the flip bifurcation (labeled “FB”) takes place for the fixed point Z3=0.8036104706,0.8036104706 with a=0.6896609431,α1=1.664046302 and α2=2.266430550>0. The chaotic attractor for a=0.2 is given in Figure 3(b) and the corresponding time sequence diagram is given in Figure 3(c) which shows the chaotic nature. Figure 3 confirms the correctness of the theoretical results in Theorem 5.Figure 3
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=0.2 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


### 5.1.2. Neimark–Sacker Bifurcation

The case of Neimark–Sacker bifurcation is given in Figure4. Suppose that parameters 0.8≤a≤1.8,r=3 and ε=3.4, based on Theorem 1, we know model (4) has only one positive fixed point. Figure 4(a) shows that the bifurcation diagram of model (4) with respect to the bifurcating parameter a in a,x plane with initial values (0.5, 0.5), the Neimark–Sacker bifurcation (labeled “NSB”) takes place for fixed point Z3=0.6206980308,0.6206980308 and a=1.611089371 with d=0.7278891375>0 and Π=−1.722080647<0, an attracting invariant cycles appears. The chaotic attractor for a=1.778 is given in Figure 4(b) and the corresponding time sequence diagram is given in Figure 4(c) which shows the chaotic nature. Figure 4 confirms the correctness of the theoretical results in Theorem 6.Figure 4
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=1.778 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


### 5.1.3. Codimension-Two Bifurcation

The case of codimension-two bifurcation is given in Figure5. Suppose that parameters a=1.307020867, r=4.733967524 and ε=3, based on Theorem 1, we know the model (4) has only one positive fixed point. After calculating, we obtain the positive fixed point Z3=0.7592142764,0.7592142764, and the eigenvalues of the Jacobian matrix JZ3 are λ1,2=−1. Furthermore, we have a=a1=1.307020867 and r=r1=4.733967524, thus a,r,ε∈CT, Ca1,r1=12.01168114≠0 and Da1,r1=107.5939727≠0. Figure 5(a) shows that the three-dimensional bifurcation of the model (4) in a,r,x plane for ε=3 and a,r vary in a neighborhood of a1,r1=1.307020867,4.733967524 with initial values (0.6, 0.6). Therefore, it follows from Theorem 7 that a fixed point Z3 is a 1 : 2 strong resonance bifurcation point. Figure 5(b) shows the two-dimensional bifurcation diagram in a,x plane when r=4.733967524. Figure 5 confirms the correctness of the theoretical results in 3.3.Figure 5
(a) Bifurcation diagram of model (4) in a,r,x plane. (b) Bifurcation diagram of model (4) in a,x plane with r=4.733967524.
(a)(b)


### 5.1.4. Marotto’s Chaos

Using Example1, the case of Marotto’s chaos is given in Figure 6. Suppose that parameters 2.54≤r≤2.57, a=2.1 and ε=4.5, based on Theorem 1, we know the model (4) has only one positive fixed point Z3. Figure 6(a) shows the bifurcation diagram in r,x plane with initial values (0.3, 0.3). The maximum Lyapunov exponents corresponding to (a) are calculated and plotted in Figure 6(b). It can be seen from Figure 6(b) that all the maximal Lyapunov exponents are larger than 0, so the model (4) with the given initial values are all chaotic. From Example 1 and Theorem 9, one can know that the chaotic attractor is located in the chaotic region r∈2.54,2.57, which is in line with the maximal Lyapunov exponents in Figure 6(b). The Marotto’s chaotic attractor is exhibited in Figure 6(c) for r=2.55. Figure 6 verifies the correctness of the theoretical results in Theorem 9.Figure 6
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Maximum Lyapunov exponents corresponding to (a). (c) Chaotic attractor for r=2.55 in (a).
(a)(b)(c)


## 5.1.1. Flip Bifurcation

The case of flip bifurcation is given in Figure3. Suppose that parameters 0.1≤a≤1.1, r=2.9 and ε=4.5, based on Theorem 1, we know model (4) has only one positive fixed point Z3. Figure 3(a) shows the bifurcation diagram of the model with respect to the bifurcating parameter a in a,x plane with initial values (0.9, 0.9), the flip bifurcation (labeled “FB”) takes place for the fixed point Z3=0.8036104706,0.8036104706 with a=0.6896609431,α1=1.664046302 and α2=2.266430550>0. The chaotic attractor for a=0.2 is given in Figure 3(b) and the corresponding time sequence diagram is given in Figure 3(c) which shows the chaotic nature. Figure 3 confirms the correctness of the theoretical results in Theorem 5.Figure 3
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=0.2 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


## 5.1.2. Neimark–Sacker Bifurcation

The case of Neimark–Sacker bifurcation is given in Figure4. Suppose that parameters 0.8≤a≤1.8,r=3 and ε=3.4, based on Theorem 1, we know model (4) has only one positive fixed point. Figure 4(a) shows that the bifurcation diagram of model (4) with respect to the bifurcating parameter a in a,x plane with initial values (0.5, 0.5), the Neimark–Sacker bifurcation (labeled “NSB”) takes place for fixed point Z3=0.6206980308,0.6206980308 and a=1.611089371 with d=0.7278891375>0 and Π=−1.722080647<0, an attracting invariant cycles appears. The chaotic attractor for a=1.778 is given in Figure 4(b) and the corresponding time sequence diagram is given in Figure 4(c) which shows the chaotic nature. Figure 4 confirms the correctness of the theoretical results in Theorem 6.Figure 4
(a) Bifurcation diagram of model (4) in a,x plane. (b) Chaotic attractor for a=1.778 in (a). (c) Time sequence diagram of prey densities in (b).
(a)(b)(c)


## 5.1.3. Codimension-Two Bifurcation

The case of codimension-two bifurcation is given in Figure5. Suppose that parameters a=1.307020867, r=4.733967524 and ε=3, based on Theorem 1, we know the model (4) has only one positive fixed point. After calculating, we obtain the positive fixed point Z3=0.7592142764,0.7592142764, and the eigenvalues of the Jacobian matrix JZ3 are λ1,2=−1. Furthermore, we have a=a1=1.307020867 and r=r1=4.733967524, thus a,r,ε∈CT, Ca1,r1=12.01168114≠0 and Da1,r1=107.5939727≠0. Figure 5(a) shows that the three-dimensional bifurcation of the model (4) in a,r,x plane for ε=3 and a,r vary in a neighborhood of a1,r1=1.307020867,4.733967524 with initial values (0.6, 0.6). Therefore, it follows from Theorem 7 that a fixed point Z3 is a 1 : 2 strong resonance bifurcation point. Figure 5(b) shows the two-dimensional bifurcation diagram in a,x plane when r=4.733967524. Figure 5 confirms the correctness of the theoretical results in 3.3.Figure 5
(a) Bifurcation diagram of model (4) in a,r,x plane. (b) Bifurcation diagram of model (4) in a,x plane with r=4.733967524.
(a)(b)


## 5.1.4. Marotto’s Chaos

Using Example1, the case of Marotto’s chaos is given in Figure 6. Suppose that parameters 2.54≤r≤2.57, a=2.1 and ε=4.5, based on Theorem 1, we know the model (4) has only one positive fixed point Z3. Figure 6(a) shows the bifurcation diagram in r,x plane with initial values (0.3, 0.3). The maximum Lyapunov exponents corresponding to (a) are calculated and plotted in Figure 6(b). It can be seen from Figure 6(b) that all the maximal Lyapunov exponents are larger than 0, so the model (4) with the given initial values are all chaotic. From Example 1 and Theorem 9, one can know that the chaotic attractor is located in the chaotic region r∈2.54,2.57, which is in line with the maximal Lyapunov exponents in Figure 6(b). The Marotto’s chaotic attractor is exhibited in Figure 6(c) for r=2.55. Figure 6 verifies the correctness of the theoretical results in Theorem 9.Figure 6
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Maximum Lyapunov exponents corresponding to (a). (c) Chaotic attractor for r=2.55 in (a).
(a)(b)(c)


## 5.2. Numerical Simulations for Model (4) in Different Bifurcation Parameters

In this part, we investigate the influence of the parameters on the dynamics of the system (4), and some new and more complex dynamic behaviors are presented as the parameters vary. Bifurcation parameters are considered in the following four cases, respectively:(1)
Varya in range 0≤a≤1.9, and fix r=3,ε=4;(2)
Varya in range 0≤a≤1.9, and fix r=3,ε=5.61;(3)
Varyε in range 1≤ε≤3.4, and fix a=1.8,r=3;(4)
Varyr in range 4≤r≤6.5, and fix a=1.8,ε=1.


### 5.2.1. Case (1)

Taking parametersr=3,ε=4 and a covering [0, 1.9] with initial values (0.66, 0.66), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 7(a) and the local amplification is given in Figure 7(b) for a∈0.3,0.8. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 7(c) which confirms the existence of the chaos and period orbits in the parametric space. In Figure 7(a), there is a stable fixed point for more large regions a∈0.714,1.57, the stable fixed point loses its stability at a∼1.57 and closed invariant circle appears when the parameter a exceeds 1.57 and eventually leads to chaos at a∼1.845. There are two different routes (invariant circle and inverse period-doubling bifurcation) to the chaos which occurs simultaneously. The inverse period-doubling bifurcation phenomena occurs at a∼0.714, and there are some cascades of period-doubling bifurcations lead to chaos. The chaotic behavior suddenly appears at a∼0. Figure 7(b) shows that there is an unending sequence of bifurcations until chaos. The phase-plane diagrams corresponding to Figure 7(a) are shown in Figures 7(d)–7(f) for the two-coexisting chaotic sets at a=0.315, four-coexisting chaotic sets at a=0.353, and the invariant circle at a=1.57799, respectively.Figure 7
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈0.3,1.2 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 6(a): (d) a=0.315, (e) a=0.353, and (f) a=1.57799.
(a)(b)(c)(d)(e)(f)


### 5.2.2. Case (2)

Taking parametersr=3,ε=5.61 and a covering [0, 1.9] with initial values (0.08, 0.08), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 8(a) and the local amplifications are plotted in Figures 8(b) and 8(c) for a∈1.25,1.4 and a∈1.4,1.5, respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 8(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 8(a) shows that there is a stable fixed point for a∈0.8075,1.173 and the fixed point loses its stability as a increases. The period-3 orbits are observed for ∈1.173,1.297, and there happen two different routes to chaos at the same time; they are inverse period-doubling bifurcation and period-doubling bifurcation. When a increases at a=1.36, a period-6 orbit appears, and some cascades of period-doubling bifurcations lead to chaos. The chaotic behaviors suddenly appear at a∼1.65. One can also observe that there is an inverse period-doubling bifurcation which occurs at a∼0.8075 changes into chaos with complex period windows and the chaotic behavior suddenly appears at a∼0. Figure 8(c) depicts the alternates of model (4) between periodic and chaotic motions. The phase-plane diagrams corresponding to Figure 8(a) are shown in Figures 8(e)–8(i) for the two-coexisting chaotic sets at a=0.42, period-3 orbits at a=1.27, period-6 orbits at a=1.36, three-coexisting chaotic sets at a=1.571, and chaotic attractor is at a=1.65, respectively.Figure 8
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈1.25,1.4 corresponding to (a). (c) Local amplification for a∈1.4,1.5 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 7(a): (e) a=0.42, (f) a=1.27, (g) a=1.36, (h) a=1.571, and (i) a=1.65.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


### 5.2.3. Case (3)

Taking parametersa=1.8,r=3 and ε covering [1, 3.4] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in ε,x plane is given in Figure 9(a). Figure 9(b) is the local amplifications for 2.93<ε<3.1. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 9(c) which confirms the existence of the chaos and period orbits in the parametric space. Figure 9(a) shows that the stability of fixed point happens for ε∈1,2.495, loses its stability as ε increases and through a Neimark–Sacker bifurcation for ε∼2.495. For ε>2.495, the fixed point becomes unstable and a closed invariant circle appears, the radius of the circle becomes bigger and bigger, and eventually leads to chaos. In particular, for ε>2.88, the invariant circle disappears, then it becomes period-4 orbits at ε∼2.976 suddenly. Model (4) undergoes period-doubling phenomena, which leads to chaos with complex period windows. The phase-plane diagrams for various values of ε are shown in Figures 9(d)–9(l), which clearly illustrate the act of how an invariant circle bifurcates from the stable fixed point and increases its radius, the invariant circle leads to chaotic attractors and the chaotic behavior suddenly disappears at ε∼3.41. Figure 9 shows that there are period-4, period-8, period-30, four-coexisting chaotic sets, and chaotic attractor.Figure 9
(a) Bifurcation diagram of model (4) in ε,x plane. (b) Local amplification for ε∈2.93,3.1 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of ε corresponding to Figure 8(a): (d) ε=1.9, (e) ε=2.5, (f) ε=2.875, (g) ε=2.88, (h) ε=2.976, (i) ε=3.048, (j) ε=3.101, (k) ε=3.131, and (l) ε=3.321.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


### 5.2.4. Case (4)

Taking parametersa=1.8,ε=1 and r covering [4, 6.5] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in r,x plane is given in Figure 10(a) and the local amplifications for r∈5.82,5.96 and r∈5.96,6.3 are given in Figures 10(b) and 10(c), respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 10(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 10(a) shows that the stability areas of a fixed point is r∈4,5.279, the fixed point becomes unstable as r increases. This would suggest that a population increases the carrying capacity at a certain growth rate and then lose its stability. The phase-plane diagrams in Figures 10(e)–10(m) clearly depicts the process of an invariant curve becoming increasingly complex as parameter r growth and chaotic attractors appear. When r increases at certain values, for example, at r>5.279, a closed curve appears, indicating that model (4) undergoes a Neimark–Sacker bifurcation, and the quasiperiodic orbit arise. At r=5.85, the circle disappears and a period-9 orbit appears, then evolves into chaos. Figure 10 shows the rich types of dynamics for the model (4); there are period-9, period-43, quasiperiodic orbit, nonattracting chaotic set, and chaotic attractor.Figure 10
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Local amplification for r∈5.82,5.96 corresponding to (a). (c) Local amplification for r∈5.96,6.3 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of r corresponding to Figure 9(a): (e) r=4.9, (f) r=5.291, (g) r=5.5711, (h) r=5.582, (i) r=5.59, (j) r=5.85, (k) r=6.055, (l) r=6.08, and (m) r=6.16.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)


## 5.2.1. Case (1)

Taking parametersr=3,ε=4 and a covering [0, 1.9] with initial values (0.66, 0.66), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 7(a) and the local amplification is given in Figure 7(b) for a∈0.3,0.8. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 7(c) which confirms the existence of the chaos and period orbits in the parametric space. In Figure 7(a), there is a stable fixed point for more large regions a∈0.714,1.57, the stable fixed point loses its stability at a∼1.57 and closed invariant circle appears when the parameter a exceeds 1.57 and eventually leads to chaos at a∼1.845. There are two different routes (invariant circle and inverse period-doubling bifurcation) to the chaos which occurs simultaneously. The inverse period-doubling bifurcation phenomena occurs at a∼0.714, and there are some cascades of period-doubling bifurcations lead to chaos. The chaotic behavior suddenly appears at a∼0. Figure 7(b) shows that there is an unending sequence of bifurcations until chaos. The phase-plane diagrams corresponding to Figure 7(a) are shown in Figures 7(d)–7(f) for the two-coexisting chaotic sets at a=0.315, four-coexisting chaotic sets at a=0.353, and the invariant circle at a=1.57799, respectively.Figure 7
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈0.3,1.2 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 6(a): (d) a=0.315, (e) a=0.353, and (f) a=1.57799.
(a)(b)(c)(d)(e)(f)


## 5.2.2. Case (2)

Taking parametersr=3,ε=5.61 and a covering [0, 1.9] with initial values (0.08, 0.08), the bifurcation diagram of model (4) with the Allee effect on prey population in a,x plane is given in Figure 8(a) and the local amplifications are plotted in Figures 8(b) and 8(c) for a∈1.25,1.4 and a∈1.4,1.5, respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 8(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 8(a) shows that there is a stable fixed point for a∈0.8075,1.173 and the fixed point loses its stability as a increases. The period-3 orbits are observed for ∈1.173,1.297, and there happen two different routes to chaos at the same time; they are inverse period-doubling bifurcation and period-doubling bifurcation. When a increases at a=1.36, a period-6 orbit appears, and some cascades of period-doubling bifurcations lead to chaos. The chaotic behaviors suddenly appear at a∼1.65. One can also observe that there is an inverse period-doubling bifurcation which occurs at a∼0.8075 changes into chaos with complex period windows and the chaotic behavior suddenly appears at a∼0. Figure 8(c) depicts the alternates of model (4) between periodic and chaotic motions. The phase-plane diagrams corresponding to Figure 8(a) are shown in Figures 8(e)–8(i) for the two-coexisting chaotic sets at a=0.42, period-3 orbits at a=1.27, period-6 orbits at a=1.36, three-coexisting chaotic sets at a=1.571, and chaotic attractor is at a=1.65, respectively.Figure 8
(a) Bifurcation diagram of the model (4) in a,x plane. (b) Local amplification for a∈1.25,1.4 corresponding to (a). (c) Local amplification for a∈1.4,1.5 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of a corresponding to Figure 7(a): (e) a=0.42, (f) a=1.27, (g) a=1.36, (h) a=1.571, and (i) a=1.65.
(a)(b)(c)(d)(e)(f)(g)(h)(i)


## 5.2.3. Case (3)

Taking parametersa=1.8,r=3 and ε covering [1, 3.4] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in ε,x plane is given in Figure 9(a). Figure 9(b) is the local amplifications for 2.93<ε<3.1. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 9(c) which confirms the existence of the chaos and period orbits in the parametric space. Figure 9(a) shows that the stability of fixed point happens for ε∈1,2.495, loses its stability as ε increases and through a Neimark–Sacker bifurcation for ε∼2.495. For ε>2.495, the fixed point becomes unstable and a closed invariant circle appears, the radius of the circle becomes bigger and bigger, and eventually leads to chaos. In particular, for ε>2.88, the invariant circle disappears, then it becomes period-4 orbits at ε∼2.976 suddenly. Model (4) undergoes period-doubling phenomena, which leads to chaos with complex period windows. The phase-plane diagrams for various values of ε are shown in Figures 9(d)–9(l), which clearly illustrate the act of how an invariant circle bifurcates from the stable fixed point and increases its radius, the invariant circle leads to chaotic attractors and the chaotic behavior suddenly disappears at ε∼3.41. Figure 9 shows that there are period-4, period-8, period-30, four-coexisting chaotic sets, and chaotic attractor.Figure 9
(a) Bifurcation diagram of model (4) in ε,x plane. (b) Local amplification for ε∈2.93,3.1 corresponding to (a). (c) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of ε corresponding to Figure 8(a): (d) ε=1.9, (e) ε=2.5, (f) ε=2.875, (g) ε=2.88, (h) ε=2.976, (i) ε=3.048, (j) ε=3.101, (k) ε=3.131, and (l) ε=3.321.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)


## 5.2.4. Case (4)

Taking parametersa=1.8,ε=1 and r covering [4, 6.5] with initial values (0.4, 0.4), the bifurcation diagram of the model (4) with the Allee effect on prey population in r,x plane is given in Figure 10(a) and the local amplifications for r∈5.82,5.96 and r∈5.96,6.3 are given in Figures 10(b) and 10(c), respectively. The maximum Lyapunov exponents corresponding to (a) are plotted in Figure 10(d) which confirms the existence of the chaos and period orbits in the parametric space. Figure 10(a) shows that the stability areas of a fixed point is r∈4,5.279, the fixed point becomes unstable as r increases. This would suggest that a population increases the carrying capacity at a certain growth rate and then lose its stability. The phase-plane diagrams in Figures 10(e)–10(m) clearly depicts the process of an invariant curve becoming increasingly complex as parameter r growth and chaotic attractors appear. When r increases at certain values, for example, at r>5.279, a closed curve appears, indicating that model (4) undergoes a Neimark–Sacker bifurcation, and the quasiperiodic orbit arise. At r=5.85, the circle disappears and a period-9 orbit appears, then evolves into chaos. Figure 10 shows the rich types of dynamics for the model (4); there are period-9, period-43, quasiperiodic orbit, nonattracting chaotic set, and chaotic attractor.Figure 10
(a) Bifurcation diagram of the model (4) in r,x plane. (b) Local amplification for r∈5.82,5.96 corresponding to (a). (c) Local amplification for r∈5.96,6.3 corresponding to (a). (d) Maximum Lyapunov exponents corresponding to (a). Phase-plane diagrams for various values of r corresponding to Figure 9(a): (e) r=4.9, (f) r=5.291, (g) r=5.5711, (h) r=5.582, (i) r=5.59, (j) r=5.85, (k) r=6.055, (l) r=6.08, and (m) r=6.16.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)


## 5.3. Coexisting Attractors and Multistability of Model (4)

As mentioned in Subsection5.2, model (4) can exhibit a variety of rich dynamical behaviors when varying the parameters a, r and ε, respectively. In this subsection, we will present some coexisting attractors and multistability to obtain more insightful information on the complex dynamics of the model (4).With a view of exploring more detailed dynamics of model (4), we first plot a two-parameter dynamical diagram [42] as illustrated in Figure 11 to present the impact of the parameters a, r, and ε. It can be clearly seen from Figure 11 that there are two color areas in which the yellow area denotes the periodic state and the dark-blue region represents the chaotic state. Figures 11(a) and 11(b) show the dynamics regions of the model (4) in the a,r plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Figures 11(c) and 11(d) show the dynamics regions of the model (4) in the a,ε plane with the initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. Observed from Figures 11(a) and 11(b), there are some finely banded yellow regions occurring in the dark-blue area, indicating that the orbits of the model (4) alter many times in these regions between period and chaos. This also shows that model (4) is extremely sensitive in these regions, meaning the emergence of the complicated state transition behavior. Furthermore, when we vary both the parameters a and r within a certain range, one can observe that the color area is always relatively concentrated and complete in the yellow region (chaos) or dark-blue region (period), which means model (4) has strong structural stability. Especially, viewing the color of the dynamical map horizontally along the direction of r=3, the dynamical behaviors are in good agreement with the bifurcation diagram of Figure 12(a). Figures 11(c) and 11(d) show that when the parameter a is small, model (4) is in a chaotic state. Then it converts into the periodic state, next drop into chaos again, and finally go out of chaos into periodic orbits with a gradually increasing. We also see that period is mainly concentrated in the left half of the graph while chaos has mainly appeared in the right half. Remark that seeing the color of the dynamical map horizontally along the direction of ε=4.8, the dynamical behaviors also coincide well with the bifurcation diagram of Figure 12(a). A comparison of Figures 11(a) and 11(b) (or Figures 11(c) and 11(d)) shows that different initial conditions may be able to have different dynamical behaviors when taking the same parameters a and r. In conclusion, the two-parameter dynamical map, as a useful tool, can offer a better parameter selection for population evolution in the real world.Figure 11
(a, b) Dynamical map with respect to the parametersa and r with ε=4.8 and initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively. (c, d) Dynamical map with respect to the parameters a and ε with r=3 and initial conditions (0.08, 0.08) and initial conditions (0.7, 0.7), respectively.
(a)(b)(c)(d)Figure 12
Coexisting dynamics with respect to the parametera covering [0,1] with initial conditions (0.08, 0.08) (cyan) and initial conditions (0.7, 0.7) (blue). (a) Bifurcation diagram of the model (4) in a,x plane. (b) Maximum Lyapunov exponents corresponding to (a).
(a)(b)Based on the previously discussed analysis, our attention is further focused on the study of coexisting attractors and multistability. The coexisting attractors is an important nonlinear dynamics which means two or more attractors generate simultaneously under different initial conditions but with the same parameter set [43]. Particularly, owing to only one positive stable fixed point, the resulting attractors are all hidden [43, 44]. The trajectories of the model with multistability can asymptotically approach different states when varying initial conditions and fixing parametric values. In order to observe this interesting phenomenon in the model (4), suppose that r=3,ε=4.8, and select a as the adjustable parameter. The coexisting bifurcation diagram with varying parameter 0≤a≤1 is displayed in Figure 12(a) and the corresponding maximum Lyapunov exponents to (a) are plotted in Figure 12(b). Note that the cyan trajectory and blue trajectory in Figure 12 begin with the initial conditions (0.08, 0.08) and (0.7, 0.7), respectively. The parametric region in which Figure 12(a) presents different behaviors is called the coexisting zone. Obviously, from Figure 12(a), one can observe the coexistence of period and chaos in an extremely narrow area between the two vertical red dashed lines with 0.191≤a≤0.209. With the parameter varying in the coexisting zone, there are chaotic attractors coexist with the periodic attractors (e.g., period-3 orbits, period-6 orbits, and period-12 orbits). When setting the parameter a=0.206 and keeping other parameters unchanged, the phase portrait of the coexistence of the period-3 attractor and chaotic attractor are illustrated in Figure 13. It is worth noting that the bifurcation diagram matches well with its corresponding maximum Lyapunov exponents.Figure 13
Phase portraits for period-3 and chaos orbits.It is well known that the coexisting attractors can be characterized by the basins of attraction, and the multistability depends on the different basins of attraction [42]. An attractor’s attraction basin is the region of the initial condition plane, which tend to the attractor as the time approaches infinity. That is when we start with points that are within the basin of attraction, they will finally be iterated into the attractor [45, 46]. Fix the parameters a=0.206,r=3 and ε=4.8. The basin of attraction of the coexisting period-3 attractor and chaotic attractor of Figure 13 is shown in Figure 14. Different colors represent different attraction regions, where the orange-red and the royal blue indicate period-3 attractor and chaotic attractor, respectively. Figure 14 depicts such model is crucially dependent on the initial states, and selecting the initial value in different colors, the model (4) will generate chaotic or periodic attractors. What is more, it’s easy to observe that the model is delicate with respect to the initial conditions. Especially, the orange-red region and royal blue region are interweaved with each other in the bottom half of the graph, implying that a tiny perturbation to the initial state is likely to lead to a state transition behavior.Figure 14
The basins of attraction in thex−y initial plane.


## 5.4. Impact of the Allee Effect on the Stability of Model (5)

For better visualization of the impact of the Allee parameteru on the stability of model (5), we draw the bifurcation diagrams of model (5) in a,x plane, ε,x plane and r,x plane, respectively. Bifurcation parameters are considered in four cases which is the same as subsection 5.2.In Figure15, the four different values of uu=0.02,0.05,0.08,0.11 are selected. Figures 15(a)–15(d) show that if we gradually increase the value of u for the model (5) with the Allee effect on both populations, then model (5) becomes stable, i.e., the Allee effect in predators has a stabilizing effect on the model (5) with the set of parameter values, the chaotic dynamics is observed at large values of a than the model without Allee constant u (see Figures 7(a) and 8(a)). However, changing the values of u has no significant effect on the inverse period-doubling bifurcation (see Figures 15(a)–15(h)). In fact, we can use the feedback control technique to change the chaotic dynamics of inverse period-doubling bifurcation (see Section 6). A similar behavior is observed in the remaining subfigures ((d)-(h)), ((i)-(l)), and ((m)-(p)) also. The results we obtained mean that the Allee constant u plays a very important role in managing the population.Figure 15
Bifurcation diagrams of model (5) with different values of u. The parameter values for a−d corresponding to case (1). The parameter values for e−h corresponding to case (2). The parameter values for i−l corresponding to case (3). The parameter values for m−p corresponding to case (4). (a) u = 0.02, (b) u = 0.05, (c) u = 0.08, (d) u = 0.011, (e) u = 0.02, (f) u = 0.05, (g) u = 0.08, (h) u = 0.011, (i) u = 0.02, (j) u = 0.05, (k) u = 0.08, (l) u = 0.011, (m) u = 0.02, (n) u = 0.05, (o) u = 0.08, and (p) u = 0.011.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)


## 6. Chaos Control

Controlling chaos is one of the most vital and developed areas for population models in the current research. Chaos due to the emergence of bifurcations is considered to be unfavorable phenomenon for the survival of populations. Hence, it is necessary to execute chaos control techniques to avoid the population from unpredictable situations. In this section, the hybrid feedback control strategy (state feedback and parameter perturbation) in references [31, 47–49] is applied to control the chaotic behavior which develops because of the appearance of bifurcation of model (4). This strategy has significant characteristics in population models, especially in discrete-time population models associated with biological species, which can also make the chaotic behavior predictable and stable.Assuming that model (4) substantiates bifurcation at the unique positive (coexistence) fixed point Z3x∗,y∗, the controlled form of model (4) is given as follows:(115)xn+1=sxn+rxn1−xn1−e−εx−axnyn+1−sxn,yn+1=syn+aynxn−yn+1−syn,where s∈0,1 is the control parameter.The Jacobian matrix of controllable model (115) estimated at Z3x∗,y∗ is given as follows:(116)Ω=∑i,j=12∂i,j2,2Ωij,then, its characteristic equation is as follows:(117)λ2−Ω11+Ω22λ+Γ=0,where(118)Ω11=−1+εx∗2−ε+2x∗rse−εx∗+−2x∗+1r−ay∗s+1,Ω12=−sax∗,Ω21=say∗,Ω22=1+sx∗−2y∗a,Γ=Ω11Ω22−Ω12Ω21.Next, we give the following theorem on the stability of the controlled model (113).Theorem 10.
The unique positive (coexistence) fixed pointZ3x∗,y∗of the controlled model (115) is locally asymptotically stable when the following condition holds:(119)Ω11+Ω22<Ω+1<2.Proof.
According to Lemma1 (i), the point is locally asymptotically stable if and only if 1−Ω11+Ω22+Γ>0,1+Ω11+Ω22+Γ>0 and Γ−1<0. Thus, Ω11+Ω22<Ω+1<2 holds.The following are two examples to demonstrate how this strategy controls the unstable fixed point.Example 2.
As seen in subsection5.1, the model (4) undergoes flip bifurcation when r=2.9,ε=4.5,a=0.6897. We apply the hybrid feedback control strategy to control chaos which appears because of the flip bifurcation. The controlled model (115) takes the following form:(120)xn+1=sxn+2.9xn1−xn1−e−4.5x−0.6897xnyn+1−sxn,yn+1=syn+0.6897ynxn−yn+1−syn.
Table1 shows that the control interval of stability is 0.1≤a≤1.1. In addition, Figure 16(a) presents the bifurcation diagram for the controlled model (120), which also indicates that the flip bifurcation of the fixed point Z3 of the controlled model (120) can be advanced.Table 1
Controllable interval for model (120) with various values of a.
Values ofaStability intervals0.110< s < 0.72842170120.210< s < 0.75980765360.310< s < 0.79500953610.410< s < 0.83555191020.510< s < 0.88315591770.610< s < 0.94141204950.630< s < 0.9548600187Figure 16
(a) Bifurcation diagram for the controlled model (120) when s=0.83. (b) Bifurcation diagram for controlled model (121) when s=0.94.
(a)(b)Example 3.
As seen in subsection5.1, the model (4) undergoes Neimark–Sacker bifurcation when r=3,ε=3.4,a=1.611. We apply the hybrid feedback control strategy to control the chaos which appears because of the Neimark–Sacker bifurcation. For this, the controlled model (113) takes the following form:(121)xn+1=sxn+3xn1−xn1−e−3.4x−1.611xnyn+1−sxn,yn+1=syn+1.611ynxn−yn+1−syn.
Table2 shows that there exists a control interval of stability for 0.8≤a≤1.8. Moreover, Figure 16(b) illustrates the bifurcation diagram for the controlled model (121), which also indicates that the Neimark–Sacker bifurcation of the fixed point Z3 of the controlled system (121) can be delayed.
By choosing the appropriate value of the control parameters, the chaotic behavior is changed into different periodic orbits. For example, when we choose =0.2, which is drawn in Figure 17 that chaotic trajectories of flip bifurcation and Neimark–Sacker bifurcation are controlled in period-1 orbit. Figures 17(a) and 17(b) show a sharp contrast to Figures 3(c) and 4(c), respectively.Table 2
Controllable interval for model (115) with various values of a.
Values ofaStability intervals1.620< s < 0.99684665391.650< s < 0.98650729901.680< s < 0.97657831981.710< s < 0.96703868301.740< s < 0.95786880721.770< s < 0.94905044711.790< s < 0.9433583728Figure 17
(a) Time sequence diagram for the statex of the controlled model (120) for r=2.9,ε=4.5,a=0.2 and s=0.2. (b) Time sequence diagram for the state x of the controlled model (121) for r=3,ε=3.4,a=1.778 and s=0.2.
(a)(b)


## 7. Summary

For many species, population growth is seasonal rather than a continuous process. Meanwhile, studying population dynamics involving Allee effects becomes an important subject of contemporary research.In this study, a nonlinear discrete-time predator-prey model with a nonlinear Allee effect has been established and analyzed in detail. By combing mathematical analysis and numerical simulations, we have shown that the model with the Allee effect in prey can undergo flip bifurcation, Neimark–Sacker bifurcation, codimension-two bifurcation associated with 1 : 2 strong resonance, and Marotto’s chaos under certain parametric conditions. With the help of a bifurcation diagram, phase portraits, time series, maximum Lyapunov exponents, dual-parameter dynamical map, the basin of attraction, the model reveals copious dynamics behaviors, including orbits of period-3, 4, 6, 8, 9, 30, and 43, invariant cycles, quasiperiodic orbits, types of two, three and four coexisting chaotic sets, and nonattracting chaotic set. More interestingly, diverse coexisting attractors and multistability are also observed when selecting appropriate parameter sets and different initial conditions. Particularly, we can discover that the Allee effect has a stabilizing effect. More precisely, an appropriate increase in the Allee effect strength on predators can stabilize population dynamics, thus reducing the probability of extinction. These results could be useful for understanding the dynamics of the discrete-time predator-prey model. Furthermore, these highly complicated dynamics of the model may provide the possibility of real-world applications [43] for the model under study and provide reasonable guidance in population dynamics.From an ecological perspective, attracting invariant closed curves, the cascade of period-doubling, and the chaotic sets imply that the prey and predator will coexist within a certain time. However, the predator-prey model is unstable if chaotic behavior happens. Chaos can lead the population to run a higher risk of extinction; thus they are harmful to the breeding of biological populations [48]. Specifically, the hybrid feedback control strategy is applied to stabilize the chaotic orbits existing in the model (4), which will make predator and prey maintain stable coexistence. Remarkably, it should be stressed that such model is often affected by external disturbances or noise [45], so there is a richness of complicated dynamical properties of this model still to be uncovered. Thus, we will pay attention to that questions in our future research work.



---

*Source: 5475999-2023-01-17.xml*