# Bifurcation Analysis and Chaos Control in a Discrete-Time Parasite-Host Model

Authors: Xueli Chen; Lishun Ren
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

A discrete-time parasite-host system with bifurcation is investigated in detail in this paper. The existence and stability of nonnegative fixed points are explored and the conditions for the existence of flip bifurcation and Neimark-Sacker bifurcation are derived by using the center manifold theorem and bifurcation theory. And we also prove the chaos in the sense of Marotto. The numerical simulations not only illustrate the consistence with the theoretical analysis, but also exhibit other complex dynamical behaviors, such as bifurcation diagrams, Maximum Lyapunov exponents, and phase portraits. More specifically, when the integral step size is chosen as a bifurcation parameter, this paper presents the finding of period orbits, attracting invariant cycles and chaotic attractors of the discrete-time parasite-host system. Specifically, we have stabilized the chaotic orbits at an unstable fixed point by using the feedback control method.


---

## Body

## 1. Introduction

In the theory of ecology, population dynamics are generally governed by continuous-time and discrete-time systems. In recent years, more and more attention is being paid to discrete-time population models [1–20]. We can get more accurate numerical simulation results and richer dynamical behaviors from discrete-time models. In addition, many scholars are interested in how to discretize continuous population dynamical models and study their dynamical properties [21–24]. The recent years, a number of articles which investigated the flip bifurcation, fold bifurcation, pitchfork bifurcation and in the sense of Marotto’s chaos of the discrete system were present in [25–29].We investigate the discrete-time parasite-host system in this paper. In some cases, parasites can reduce host density and induce host population extinction. Ebert et al. [30] formulated the following epidemiological microparasite model:(1)x˙t=ax+θy1-x+yK-dx-βxy,y˙t=βxy-d+αy,to understand how six microparasites regulate Daphnia populations and drive the populations to extinction, where x(t) and y(t) denote uninfected (susceptible) and infected (infected) hosts densities at time t, respectively; a is the maximum per capita birth rate of uninfected hosts; θ(0⩽θ⩽1) is the relative fecundity of an infected host; K is the carrying capacity of for the host population; d is the parasite-independent host background mortality; β is the constant infection rate and α is the parasite-induced excess death rate; all parameters are positive.In this model, the microparasite transmission is assumed via a mass action process, the fecundity of uninfected host is density-dependent, and the fecundity of infected host may be reduced due to being infected compared with that of uninfected host. Here,0≤θ≤1 means that the fecundity of infected hosts is reduced, but they still have certain fecundity, θ=0 means that infected hosts completely lose fecundity; θ=1 means that the fecundity of infected hosts is not affected by parasite infection. For model (1), when a>d, there is always equilibrium (0, 0) which is a saddle. This implies that extinction of host is impossible when a>d, that is to say, when a>d host always persists. But, in this model, the simulation for the stochastic model indicates that extinction of host likely occurs in certain parameter regions.Rescaling system (1) by the following as [21, 31],(2)x¯=xK,y¯=yK,t¯=atand removing the bars, then system (1) becomes in the following form:(3)x˙t=x+θy1-x-y-δx-sxy,y˙t=sxy-δ+γy,where s=βK/a, δ=d/a, and γ=α/a are positive constants. Applying the forward Euler discrete scheme to system (3), we obtain the following discrete-time parasite-host system as follows:(4)F:xy⟶x+hx+θy1-x-y-δx-sxyy+hsxy-δ+γy,where h is the integral step size. The fixed points of system (4) satisfy the following equations:(5)x+hx+θy1-x-y-δx-sxy=x,y+hsxy-δ+γy=y,and then(6)x+θy1-x-y-δx-sxy=0,sxy-δ+γy=0.In our paper, we apply the forward Euler Scheme to discrete the parasite-past model and mainly focus on the existence and stability of nonnegative fixed points and flip bifurcation, Neimark-Sacker bifurcation, and possible chaos in the sense of Marotto’s definition [32] in the discrete-time parasite-host system by using the center manifold theorem [33] and the bifurcation theory [17–20, 33–36]. When the integral step size h is taken as a bifurcation parameter, the detailed existence conditions of flip bifurcation and Neimark-Sacker bifurcation are given in a very strict mathematical way. Numerical simulations are shown, including maximum Lyapunov exponents, bifurcation diagrams, and phase portraits, to verify theoretical analyses, and display some new and interesting nonlinear dynamical behaviors of the parasite-host system. In particular, one of our investigations demonstrates that the integral step size makes difference with respect to nonlinear dynamical behaviors of the discrete parasite-host system when the integral step size is taken into account as a bifurcation parameter.The organization of this paper is as follows. In Section2, the theorem on the existence and stability of fixed points for system (4) is proved. In Section 3, we give the sufficient conditions of existence for flip bifurcation and Neimark-Sacker bifurcation by using the center manifold theorem and bifurcation theory. In Section 4, we first rigorously prove the existence of chaos in the sense of Marotto’s definition. Numerical simulations are presented not only to verify the theoretical analysis but also to exhibit other complex dynamics in Section 5. In Section 6, chaos is controlled to an unstable fixed point using the feedback control method. Finally, we conclude this paper with comments and discuss the future work.


## 2. The Existence and Stability of Fixed Points

In this section, we discuss the existence and stability of the fixed points. It is clear that the fixed points of system (4) satisfy the following equations:(7)x=x+hx+θy1-x-y-δx-sxy,y=y+hsxy-δ+γy.By calculation of the above system, we get the following results: system (4) has three fixed points (0,0), (1-δ,0) and the positive fixed point (x∗,y∗), where (x∗,y∗) satisfy(8)x∗+θy∗1-x∗-y∗-δx∗-sx∗y∗=0,x∗=δ+γs.Now we study the stability of these fixed points. Note that the modules of eigenvalues of the characteristic equation at the fixed point determine the local stability of a fixed point(x,y).The Jacobian matrixJ of system (4) at any point (x,y) is given by(9)Jx,y=1+h1-2x-θ+1+sy-δhθ1-x-y-x+θy-sxhsy1+hsx-δ-γ.Let(10)a11=1+h1-x-y-x-θy-δ-sy;a12=hθ1-x-y-x+θy-sx;a21=hsy;a22=1+hsx-δ-γ,and then(11)Jx,y=a11a12a21a22.The characteristic equation of J(x,y) is(12)λ2+px,yλ+qx,y=0,where(13)px,y=-a11+a22,qx,y=a11a22-a12a21.Then we will discuss the local dynamics of three fixed points(0,0), (1-δ,0) and the positive fixed point (x∗,y∗). In order to discuss the stability of the fixed points of system (4), we get the following existence proposition of the fixed points by simple analysis [17–20].Case 1.
For(0,0), we have(14)J0,0=1+h1-δhθ01-hδ+γ.Since two eigenvalues of J(0,0) are λ1=1+h(1-δ)>1 and λ2=1-h(δ+γ)<1, then (0,0) is a saddle.Case 2.
A simple calculation shows the local dynamics of the fixed point(1-δ,0).Proposition 1.
The eigenvalues of the fixed point(1-δ,0) are λ1=1+h(δ-1), λ2=1+h[s(1-δ)-δ-γ]:(i)
(1-δ,0) is a source if 0<δ<1, δ+γ>s(1-δ) and h>max⁡2/1-δ,2/δ+γ-s(1-δ);(ii)
(1-δ,0) is a sink if 0<δ<1, δ+γ>s(1-δ) and 0<h<min⁡2/1-δ,2/δ+γ-s(1-δ);(iii)
(1-δ,0) is not hyperbolic if 0<δ<1, δ+γ>s(1-δ) and h=2/1-δ or h=2/δ+γ-s(1-δ);(iv)
(1-δ,0) is a saddle except for any values of parameters as (i) to (iii).Case 3.
The characteristic equation of the Jacobian matrixJ of system (4) at the positive fixed point (x∗,y∗) is written as(15)λ2+tr⁡Jx∗,y∗λ+det⁡Jx∗,y∗=0,where(16)tr⁡Jx∗,y∗=-2-hM,det⁡Jx∗,y∗=1+hM+h2N,M=1-2x∗-θ+1+sy∗-δ,N=sy∗x∗+θy∗+sx∗-θ1-x∗-y∗.Let(17)Fλ=λ2-2+hMλ+1+hM+h2N,and then(18)F1=h2N>0,F-1=4+2hM+h2N.Proposition 2.
Let(x∗,y∗) be the positive fixed point of (4):(i)
it is a sink if one of the following conditions holds:
(i.1)
-2N⩽M<0 and 0<h<-M/N;(i.2)
M<-2N and 0<h<-M-M2-4N/N;(ii)
it is a source if one of the following conditions holds:
(ii.1)
-2N⩽M<0 and h>-M/N;(ii.2)
M<-2N and h>-M+M2-4N/N;(ii.3)
M⩾0;(iii)
it is a saddle if the following condition holds:(19)M<-2N,-M-M2-4NN<h<-M+M2-4NN;(iv)
it is nonhyperbolic if one of the following conditions holds:
(iv.1)
M<-2N and h=-M±M2-4N/N and δ≠-2/M,-4/M;(iv.2)
-2N⩽M<0 and h=-M/N.
From the above discussion, we can easily obtain that if (iv.1) of Proposition2 holds, then one of the eigenvalues of the positive fixed point (x∗,y∗) is -1 and the other is neither 1 nor -1.
Let(20)F1=θ,δ,s,γ,h:M<-2N,h=-M-M2-4NN,θ,δ,s,γ,h>0and(21)F2=θ,δ,s,γ,h:M<-2N,h=-M+M2-4NN,θ,δ,s,γ,h>0.When parameters vary in a small neighborhood of F1 or F2, there may be flip bifurcation of the fixed point (x∗,y∗).
Then we can get that the eigenvalues of the positive fixed point(x∗,y∗) are a pair of conjugate complex numbers with modulus 1 if (iv.2) of Proposition 2 holds.
Let(22)F3=θ,δ,s,γ,h:-2N<M<0,h=-MN,θ,δ,s,γ,h>0.The fixed point (x∗,y∗) may undergo Neimark-Sacker bifurcation when parameters vary in a small neighborhood of F3.


## 3. Bifurcation Analysis

Based on the above analysis, we will mainly focus on the flip bifurcation of the positive fixed point(x∗,y∗) if parameters vary in a small neighborhood of F1 or F2 and the Neimark-Sacker bifurcation of (x∗,y∗) if parameters vary in a small neighborhood of F3, respectively, to investigate the local and global stability of system (4). We choose parameter h as a bifurcation parameter for studying the flip bifurcation and Neimark-Sacker bifurcation of (x∗,y∗) by using the center manifold theorem and bifurcation theory in [17–20, 33–36].


### 3.1. Flip Bifurcation

Taking parameters(θ,δ,s,γ,h1) arbitrarily from F1, we consider system (4) with (θ,δ,s,γ,h1), which is described by(23)xy⟶x+h1x+θy1-x-y-δx-sxyy+h1sxy-δ+γy.The map (23) has a unique positive fixed point (x∗,y∗), whose eigenvalues are λ1=-1, λ2=3+Mh1 with λ2≠1 by Proposition 2.Since(θ,δ,s,γ,h1)∈F1, h1=-M-M2-4N/N. Choose h∗ as a bifurcation parameter, we consider a perturbation of (23) as follows:(24)xy⟶x+h1+h∗x+θy1-x-y-δx-sxyy+h1+h∗sxy-δ+γy,where h∗≪1, which is a small perturbation parameter.LetX=x-x∗, Y=y-y∗. Then we transform the fixed point (x∗,y∗) of map (24) into the origin. By calculating we have(25)XY⟶a11X+a12Y+a13X2+a14XY+a15Y2+b1Xh∗+b2Yh∗+b3X2h∗+b4XYh∗+b5Y2h∗+OX+Y+h∗3a21X+a22Y+a23XY+c1Xh∗+c2Yh∗+c3XYh∗+OX+Y+h∗3,where(26)a11=1+h1-2x∗-θ+1y∗-δ-sy∗,a12=hθ1-x∗-y∗-x∗-θy∗-sx∗,a13=-h,a14=-h1+θ+s,a15=-θh,b1=1-2x∗-θ+1y∗-δ-sy∗,b2=θ1-x∗-y∗-x∗-θy∗-sx∗,b3=-1,b4=-1+θ+s,b5=-θ,a21=hsy∗,a22=1+hsx∗,a23=hs,c1=sy∗,c2=sx∗,c3=s,and h=h1.We construct an invertible matrixT as follows:(27)T=a12a12-1-a11λ2-a11,and for (25), use the translation(28)XY=Tx~y~,and then the map (25) becomes(29)x~y~⟶-100λ2x~y~+fX,Y,h∗gX,Y,h∗,where(30)fX,Y,h∗=a13λ2-a11a12λ2+1X2+a14λ2-a11+a12a23a12λ2+1XY+a15λ2-a11a12λ2+1Y2+b1λ2-a11-a12c1a12λ2+1h∗X+b2λ2-a11-a12c2a12λ2+1h∗Y-λ2-a11a12λ2+1h∗X2+b4λ2-a11-a12c3a12λ2+1h∗XY+b5λ2-a11a12λ2+1h∗Y2+OX+Y+h∗3,gX,Y,h∗=a13a11+1a12λ2+1X2+a14a11+1+a12a23a12λ2+1XY+a15a11+1a12λ2+1Y2+b1a11+1+a12c1a12λ2+1h∗X+b2a11+1+a12c2a12λ2+1h∗Y-a11+1a12λ2+1h∗X2+b4a11+1+a12c3a12λ2+1h∗XY+b5a11+1λ2-a11a12λ2+1h∗Y2+OX+Y+h∗3,and(31)X=a12x~+a12y~,Y=-1+a11x~+λ2-a11y~.From the center manifold theorem [33], we can determine the center manifold Wc(0,0,0) of (29) at the fixed point (0,0) in a small neighborhood of h∗=0, and then we obtain that there exists a center manifold Wc(0,0,0), which can be approximately represented as follows:(32)Wc0,0,0=x~,y~,h∗∈R3:y~=a1x~2+a2x~h∗+a3h∗2+Ox~+h∗3,where O((|x~|+|h∗|)3) is a function with order at least 3 in the variables, and(33)a1=-a12a13a11+1-a12a11+12a14a11+1-a12a23-a15a11+13a121-λ22,a2=a11+12b2+a11+1a12c2-a11+1a12b1-a122c1λ2+12a12,a3=0.Therefore, the map (29) which is restricted to the center manifold Wc(0,0,0) is defined:(34)F:x~⟶-x~+k1x~2+k2x~h∗+k3x~2h∗+k4x~h∗2+k5x~3+Ox~+h∗4,where(35)k1=λ2-a11a13λ2+1+a11+1a12a23-a14λ2-a11λ2+1+a11+12λ2-a11a15λ2+1a12,k2=λ2-a11b1-a12c1λ2+1-a11+1b2λ2-a11-a12c2λ2+1a12,k3=λ2-a11b1-a12c1λ2+1a1+b2λ2-a11-a12c2λ2+1a12λ2-a11a1+1+a112λ2-a11b5a12λ2+1+2a13a12a1λ2+1-4a15a21+a11λ2-a112a12λ2+1+λ2-a11a14-a12a23a12λ2+1a12λ2-a11-a111+a112a2,k4=λ2-a11b1-a12c1λ2+1a1+b2λ2-a11-a12c2λ2+1a12λ2-a11a2,k5=2a13a12a1λ2+1+λ2-a11a14-a12a23a12λ2+1a12λ2-a11-a111+a11a1-2a15a11+a11λ2-a112a12λ2+1.For map (34) in order to undergo a flip bifurcation, we require that two discriminatory quantities α1 and α2 are not zero, where (36)α1=∂2F∂x~∂h∗+12∂F∂h∗∂2F∂x~20,0=k2and(37)α2=16∂3F∂x~3+12∂2F∂x~220,0=k5+k12.On the basis of the above analysis and the theorem of [34–36], we obtain the following result.Theorem 3.
Ifα2≠0, then map (4) undergoes a flip bifurcation at the fixed point (x∗,y∗) when the parameter h varies in a small neighborhood of h1. Moreover, if α2>0 (resp., α2<0), the period-2 orbits that bifurcate from (x∗,y∗) are stable (resp., unstable).


### 3.2. Neimark-Sacker Bifurcation

In this section, we can give a similar argument with the flip bifurcation. Taking parameters(θ,δ,s,γ,h2) arbitrarily from F3, we consider system (4) with (θ,δ,s,γ,h2), which is described by (38)xy⟶x+h2x+θy1-x-y-δx-sxyy+h2sxy-δ+γy.The map (38) has a unique positive fixed point (x∗,y∗).Since(θ,δ,s,γ,h2)∈F3,h2=-M/N. We give a perturbation h¯∗ at h2. Then (38) becomes into the following form:(39)xy⟶x+h2+h¯∗x+θy1-x-y-δx-sxyy+h2+h¯∗sxy-δ+γy,where |h¯∗|≪1, which is a small perturbation parameter.LetX=x-x∗, Y=y-y∗. Then we transform the fixed point (x∗,y∗) of map (39) into the origin, and we have(40)XY⟶a11X+a12Y+a13X2+a14XY+a15Y2+OX+Y+h∗3a21X+a22Y+a23XY+OX+Y+h∗3,where a11,a12,a13,a14,a15,a21,a22,a23 are given in (26) by substituting h for h2+h¯∗.Note that the characteristic equation associated with the linearization of the map (40) at (X,Y)=(0,0) is following(41)λ2+ph¯∗λ+qh¯∗=0,where(42)ph¯∗=-2-h+h¯∗M,qh¯∗=1+h+h¯∗M+h+h¯∗2N.Since parameters(θ,δ,s,γ,h2)∈F3, the eigenvalues of (0,0) are a pair of complex conjugate numbers λ and λ¯ with modulus 1 by Proposition 2, where(43)λ,λ¯=-ph¯∗2±i24qh¯∗-p2h¯∗=1+Mh2+h¯∗2±ih2+h¯∗24N-M2,and we have(44)λ=qh¯∗,l=dλdh¯∗h¯∗=0=-M2>0.Moreover, it is required that whenh¯∗=0, λm, λ¯m≠1 (m≠1,2,3,4) which is equivalent to p(0)≠-2,0,1,2. Note that (θ,δ,s,γ,h2)∈F3. Thus, p(0)≠-2,2. We only need to require that p(0)≠0,1, which leads to(45)M2≠2N,3N.Hence, the eigenvaluesλ,λ¯ of fixed point (0,0) of (40) do not lie in the intersection of the unit circle with the coordinate axes when h¯∗=0 and (45) holds.Next, we discuss the normal form of (40) at h¯∗=0.Leth¯∗=0, μ=1+Mh2/2, ω=h2/24N-M2,(46)T=a120μ-a11-ω,and then T is invertible and use the translation(47)XY=Tx~y~for (40); then model (40) becomes the following form:(48)x~y~⟶μ-ωωμx~y~+f~x~,y~g~x~,y~,where(49)f~x~,y~=a13a12X2+a14a12XY+a15a12Y2+OX+Y+h∗3,g~x~,y~=a13μ-a11a12ωX2+a14μ-a11a12ω-a23ωXY+a15μ-a11a12ωY2+OX+Y+h∗3,X=a12x~,Y=μ-a11x~-ωy~.In order to undergo Neimark-Sacker bifurcation for (48), we require that the following discriminatory quantity is not zero [34–36]:(50)b=-Re⁡1-2λλ¯21-λξ20ξ11-12ξ112-ξ022+Re⁡λ¯ξ21h¯∗=0,where(51)ξ20=18f~x~x~-f~y~y~+2g~x~y~+ig~x~x~-g~y~y~-2f~x~y~,ξ11=14f~x~x~+f~y~y~+ig~x~x~+g~y~y~,ξ02=18f~x~x~-f~y~y~-2g~x~y~+ig~x~x~-g~y~y~+2f~x~y~,ξ21=116f~x~x~x~+f~x~y~y~+g~x~x~y~+g~y~y~y~+ig~x~x~x~+g~x~y~y~-f~x~x~y~-f~y~y~y~.Thus, some complicated calculation gives(52)b=-1321-μ2+ω2mf~x~x~2-g~x~x~2+g~y~y~2+2f~x~x~g~x~y~+2f~x~y~g~x~x~+f~x~y~g~y~y~+n2f~x~x~g~x~x~-2f~x~x~f~x~y~+2g~x~x~g~x~y~+2g~y~y~g~x~y~-1643f~x~x~2+3g~x~x~2+3g~y~y~2+4f~x~y~2+4g~x~y~2-4f~x~x~g~x~y~-4f~x~y~g~y~y~+4f~x~y~g~x~x~+2g~x~x~g~y~y~-4μf~x~x~x~+g~x~x~y~-4ωg~x~x~x~+g~x~y~y~-f~x~x~y~,where(53)m=μ2-3μ3+2μ4-ω2+μω2-2ω4,n=-ωω2+5μ2-2μ-4μ3-4μω2,f~x~x~=2a12a13+a14μ-a11+a15a12μ-a112,f~x~y~=a14ω-2a15ωμ-a11a12,f~y~y~=2a15ω2a12,f~x~x~x~=f~x~x~y~=f~x~y~y~=f~y~y~y~=0,g~x~x~=2a12a13μ-a11ω+a12μ-a11a14μ-a11a12ω-a23ω+a15μ-a113a12ω,g~x~y~=-a14μ-a11a12ω+a23a12-2a15μ-a112a12,g~y~y~=2a15μ-a11ωa12,g~x~x~x~=g~x~x~y~=g~x~y~y~=g~y~y~y~=0.Therefore, from the above analysis and the theorem in [34–36], we have the following result.Theorem 4.
If condition (45) holds and b≠0, then map (4) undergoes Neimark-Sacker bifurcation at the fixed point (x∗,y∗) when the parameter h varies in a small neighborhood of h2. Moreover, if b<0, (resp., b>0), then an attracting (resp., repelling) invariant closed curve bifurcates from the fixed point for h>h2 (resp., h<h2).


## 3.1. Flip Bifurcation

Taking parameters(θ,δ,s,γ,h1) arbitrarily from F1, we consider system (4) with (θ,δ,s,γ,h1), which is described by(23)xy⟶x+h1x+θy1-x-y-δx-sxyy+h1sxy-δ+γy.The map (23) has a unique positive fixed point (x∗,y∗), whose eigenvalues are λ1=-1, λ2=3+Mh1 with λ2≠1 by Proposition 2.Since(θ,δ,s,γ,h1)∈F1, h1=-M-M2-4N/N. Choose h∗ as a bifurcation parameter, we consider a perturbation of (23) as follows:(24)xy⟶x+h1+h∗x+θy1-x-y-δx-sxyy+h1+h∗sxy-δ+γy,where h∗≪1, which is a small perturbation parameter.LetX=x-x∗, Y=y-y∗. Then we transform the fixed point (x∗,y∗) of map (24) into the origin. By calculating we have(25)XY⟶a11X+a12Y+a13X2+a14XY+a15Y2+b1Xh∗+b2Yh∗+b3X2h∗+b4XYh∗+b5Y2h∗+OX+Y+h∗3a21X+a22Y+a23XY+c1Xh∗+c2Yh∗+c3XYh∗+OX+Y+h∗3,where(26)a11=1+h1-2x∗-θ+1y∗-δ-sy∗,a12=hθ1-x∗-y∗-x∗-θy∗-sx∗,a13=-h,a14=-h1+θ+s,a15=-θh,b1=1-2x∗-θ+1y∗-δ-sy∗,b2=θ1-x∗-y∗-x∗-θy∗-sx∗,b3=-1,b4=-1+θ+s,b5=-θ,a21=hsy∗,a22=1+hsx∗,a23=hs,c1=sy∗,c2=sx∗,c3=s,and h=h1.We construct an invertible matrixT as follows:(27)T=a12a12-1-a11λ2-a11,and for (25), use the translation(28)XY=Tx~y~,and then the map (25) becomes(29)x~y~⟶-100λ2x~y~+fX,Y,h∗gX,Y,h∗,where(30)fX,Y,h∗=a13λ2-a11a12λ2+1X2+a14λ2-a11+a12a23a12λ2+1XY+a15λ2-a11a12λ2+1Y2+b1λ2-a11-a12c1a12λ2+1h∗X+b2λ2-a11-a12c2a12λ2+1h∗Y-λ2-a11a12λ2+1h∗X2+b4λ2-a11-a12c3a12λ2+1h∗XY+b5λ2-a11a12λ2+1h∗Y2+OX+Y+h∗3,gX,Y,h∗=a13a11+1a12λ2+1X2+a14a11+1+a12a23a12λ2+1XY+a15a11+1a12λ2+1Y2+b1a11+1+a12c1a12λ2+1h∗X+b2a11+1+a12c2a12λ2+1h∗Y-a11+1a12λ2+1h∗X2+b4a11+1+a12c3a12λ2+1h∗XY+b5a11+1λ2-a11a12λ2+1h∗Y2+OX+Y+h∗3,and(31)X=a12x~+a12y~,Y=-1+a11x~+λ2-a11y~.From the center manifold theorem [33], we can determine the center manifold Wc(0,0,0) of (29) at the fixed point (0,0) in a small neighborhood of h∗=0, and then we obtain that there exists a center manifold Wc(0,0,0), which can be approximately represented as follows:(32)Wc0,0,0=x~,y~,h∗∈R3:y~=a1x~2+a2x~h∗+a3h∗2+Ox~+h∗3,where O((|x~|+|h∗|)3) is a function with order at least 3 in the variables, and(33)a1=-a12a13a11+1-a12a11+12a14a11+1-a12a23-a15a11+13a121-λ22,a2=a11+12b2+a11+1a12c2-a11+1a12b1-a122c1λ2+12a12,a3=0.Therefore, the map (29) which is restricted to the center manifold Wc(0,0,0) is defined:(34)F:x~⟶-x~+k1x~2+k2x~h∗+k3x~2h∗+k4x~h∗2+k5x~3+Ox~+h∗4,where(35)k1=λ2-a11a13λ2+1+a11+1a12a23-a14λ2-a11λ2+1+a11+12λ2-a11a15λ2+1a12,k2=λ2-a11b1-a12c1λ2+1-a11+1b2λ2-a11-a12c2λ2+1a12,k3=λ2-a11b1-a12c1λ2+1a1+b2λ2-a11-a12c2λ2+1a12λ2-a11a1+1+a112λ2-a11b5a12λ2+1+2a13a12a1λ2+1-4a15a21+a11λ2-a112a12λ2+1+λ2-a11a14-a12a23a12λ2+1a12λ2-a11-a111+a112a2,k4=λ2-a11b1-a12c1λ2+1a1+b2λ2-a11-a12c2λ2+1a12λ2-a11a2,k5=2a13a12a1λ2+1+λ2-a11a14-a12a23a12λ2+1a12λ2-a11-a111+a11a1-2a15a11+a11λ2-a112a12λ2+1.For map (34) in order to undergo a flip bifurcation, we require that two discriminatory quantities α1 and α2 are not zero, where (36)α1=∂2F∂x~∂h∗+12∂F∂h∗∂2F∂x~20,0=k2and(37)α2=16∂3F∂x~3+12∂2F∂x~220,0=k5+k12.On the basis of the above analysis and the theorem of [34–36], we obtain the following result.Theorem 3.
Ifα2≠0, then map (4) undergoes a flip bifurcation at the fixed point (x∗,y∗) when the parameter h varies in a small neighborhood of h1. Moreover, if α2>0 (resp., α2<0), the period-2 orbits that bifurcate from (x∗,y∗) are stable (resp., unstable).


## 3.2. Neimark-Sacker Bifurcation

In this section, we can give a similar argument with the flip bifurcation. Taking parameters(θ,δ,s,γ,h2) arbitrarily from F3, we consider system (4) with (θ,δ,s,γ,h2), which is described by (38)xy⟶x+h2x+θy1-x-y-δx-sxyy+h2sxy-δ+γy.The map (38) has a unique positive fixed point (x∗,y∗).Since(θ,δ,s,γ,h2)∈F3,h2=-M/N. We give a perturbation h¯∗ at h2. Then (38) becomes into the following form:(39)xy⟶x+h2+h¯∗x+θy1-x-y-δx-sxyy+h2+h¯∗sxy-δ+γy,where |h¯∗|≪1, which is a small perturbation parameter.LetX=x-x∗, Y=y-y∗. Then we transform the fixed point (x∗,y∗) of map (39) into the origin, and we have(40)XY⟶a11X+a12Y+a13X2+a14XY+a15Y2+OX+Y+h∗3a21X+a22Y+a23XY+OX+Y+h∗3,where a11,a12,a13,a14,a15,a21,a22,a23 are given in (26) by substituting h for h2+h¯∗.Note that the characteristic equation associated with the linearization of the map (40) at (X,Y)=(0,0) is following(41)λ2+ph¯∗λ+qh¯∗=0,where(42)ph¯∗=-2-h+h¯∗M,qh¯∗=1+h+h¯∗M+h+h¯∗2N.Since parameters(θ,δ,s,γ,h2)∈F3, the eigenvalues of (0,0) are a pair of complex conjugate numbers λ and λ¯ with modulus 1 by Proposition 2, where(43)λ,λ¯=-ph¯∗2±i24qh¯∗-p2h¯∗=1+Mh2+h¯∗2±ih2+h¯∗24N-M2,and we have(44)λ=qh¯∗,l=dλdh¯∗h¯∗=0=-M2>0.Moreover, it is required that whenh¯∗=0, λm, λ¯m≠1 (m≠1,2,3,4) which is equivalent to p(0)≠-2,0,1,2. Note that (θ,δ,s,γ,h2)∈F3. Thus, p(0)≠-2,2. We only need to require that p(0)≠0,1, which leads to(45)M2≠2N,3N.Hence, the eigenvaluesλ,λ¯ of fixed point (0,0) of (40) do not lie in the intersection of the unit circle with the coordinate axes when h¯∗=0 and (45) holds.Next, we discuss the normal form of (40) at h¯∗=0.Leth¯∗=0, μ=1+Mh2/2, ω=h2/24N-M2,(46)T=a120μ-a11-ω,and then T is invertible and use the translation(47)XY=Tx~y~for (40); then model (40) becomes the following form:(48)x~y~⟶μ-ωωμx~y~+f~x~,y~g~x~,y~,where(49)f~x~,y~=a13a12X2+a14a12XY+a15a12Y2+OX+Y+h∗3,g~x~,y~=a13μ-a11a12ωX2+a14μ-a11a12ω-a23ωXY+a15μ-a11a12ωY2+OX+Y+h∗3,X=a12x~,Y=μ-a11x~-ωy~.In order to undergo Neimark-Sacker bifurcation for (48), we require that the following discriminatory quantity is not zero [34–36]:(50)b=-Re⁡1-2λλ¯21-λξ20ξ11-12ξ112-ξ022+Re⁡λ¯ξ21h¯∗=0,where(51)ξ20=18f~x~x~-f~y~y~+2g~x~y~+ig~x~x~-g~y~y~-2f~x~y~,ξ11=14f~x~x~+f~y~y~+ig~x~x~+g~y~y~,ξ02=18f~x~x~-f~y~y~-2g~x~y~+ig~x~x~-g~y~y~+2f~x~y~,ξ21=116f~x~x~x~+f~x~y~y~+g~x~x~y~+g~y~y~y~+ig~x~x~x~+g~x~y~y~-f~x~x~y~-f~y~y~y~.Thus, some complicated calculation gives(52)b=-1321-μ2+ω2mf~x~x~2-g~x~x~2+g~y~y~2+2f~x~x~g~x~y~+2f~x~y~g~x~x~+f~x~y~g~y~y~+n2f~x~x~g~x~x~-2f~x~x~f~x~y~+2g~x~x~g~x~y~+2g~y~y~g~x~y~-1643f~x~x~2+3g~x~x~2+3g~y~y~2+4f~x~y~2+4g~x~y~2-4f~x~x~g~x~y~-4f~x~y~g~y~y~+4f~x~y~g~x~x~+2g~x~x~g~y~y~-4μf~x~x~x~+g~x~x~y~-4ωg~x~x~x~+g~x~y~y~-f~x~x~y~,where(53)m=μ2-3μ3+2μ4-ω2+μω2-2ω4,n=-ωω2+5μ2-2μ-4μ3-4μω2,f~x~x~=2a12a13+a14μ-a11+a15a12μ-a112,f~x~y~=a14ω-2a15ωμ-a11a12,f~y~y~=2a15ω2a12,f~x~x~x~=f~x~x~y~=f~x~y~y~=f~y~y~y~=0,g~x~x~=2a12a13μ-a11ω+a12μ-a11a14μ-a11a12ω-a23ω+a15μ-a113a12ω,g~x~y~=-a14μ-a11a12ω+a23a12-2a15μ-a112a12,g~y~y~=2a15μ-a11ωa12,g~x~x~x~=g~x~x~y~=g~x~y~y~=g~y~y~y~=0.Therefore, from the above analysis and the theorem in [34–36], we have the following result.Theorem 4.
If condition (45) holds and b≠0, then map (4) undergoes Neimark-Sacker bifurcation at the fixed point (x∗,y∗) when the parameter h varies in a small neighborhood of h2. Moreover, if b<0, (resp., b>0), then an attracting (resp., repelling) invariant closed curve bifurcates from the fixed point for h>h2 (resp., h<h2).


## 4. Existence of Marotto’s Chaos

In this section, we rigorously prove map (4) possesses chaotic behavior in the sense of Marotto’s definition [32].We first present Marotto’s chaos definitions and theorem which are quoted from [32].For any mapF:Rn→Rn, and any positive integer K, let FK represent the composition of F with itself K times. For a differentiable function F, let DF(Z) denote the Jacobian matrix of F evaluated at the point Z∈Rn, and DF(Z) its determinant. Let Br(Z) denote the closed ball in Rn of radius r centered at the point Z and Br0(Z) its interior. Also let Z be the usual Euclidean norm of Z in Rn.Definition 5.
LetF be differentiable in Br(Z0). The point Z0∈Rn is an expanding fixed point of F in Br(Z0), if F(Z0)=Z0 and all eigenvalues of DF(Z) exceed 1 in norm for all Z∈Br(Z0).Definition 6.
Assume thatZ0 is an expanding fixed point of F in Br(Z0) for some r>0; then Z0 is said to be a snap-back repeller of F if there exists a point Z¯∈Br(Z0) with Z¯≠Z0, FM(Z¯)=Z0 and DFM(Z¯)≠0 for some positive integer M.Theorem 7.
[32] If F possesses a snap-back repeller, then the map F is chaotic. That is, there exist(i)
a positive integerN such that for each integer P⩾N, F has a point of period P;(ii)
a “scrambled set” ofF, that is, an uncountable set S containing no periodic points of F such that
(a)
F[S]⊂S,(b)
for everyX,Y∈S with X≠Y(54)limk→∞sup⁡FkX-FkY>0,(c)
for everyX∈S and any periodic point Y of F(55)limk→∞sup⁡FkX-FkY>0,(iii)
an uncountable subsetS0 of S such that for every X,Y∈S0(56)limk→∞inf⁡FkX-FkY=0.
Now we theoretically give the condition of existence of chaotic phenomena for map (4) in the sense of Marotto’s definition of chaos.
SupposeZ0(x0,y0) be the fixed point of map (4). We firstly give the conditions such that the Z0 is a snap-back repeller. The eigenvalues associated with the fixed point Z0 are given by (57)λ1,2=-px0,y0±p2x0,y0-4qx0,y02,where (58)px,y=-2-h1-2x-θ+1+sy-δ+sx-δ-γ,qx,y=1+h1-2x-θ+1+sy-δ+sx-δ-γ+h21-2x-θ+1+sy-δsx-δ-γ-syθ1-x-y-x+θy-sx.
According to Definition5, we begin to find a neighborhood Br(Z0) of Z0 in which the norms of conjugate complex eigenvalues exceed 1 for all Z∈Br(Z0), which are equivalent to (59)p2x,y-4qx,y<0,qx,y-1>0.
Lets1x,y=p2x,y-4qx,y.
It is easy to see thats1(x,y)=h2[A1y2+B1y+C1], where (60)A1=θ+1+s2-8sθ,B1=2θ+1+s2x-γ-1-sx+4sθ,C1=1-2x-δ-sx-δ-γ2.If Δ1=B12-4A1C1=(4sθ)2+16sθ(θ+1+s)(2x-γ-1-sx)+32sθ(1-2x-sx+γ)2⩾0,
then the equations1(x,y)=0 has one real root with multiplicity 2 or two real roots denoted as y¯1=-B1-B12-4A1C1/2A1 and y¯2=-B1+B12-4A1C1/2A1.
AndΔ1⩾0; we can get (61)θ+s+12-8sθ>0sθ+2x-γ-1-sxθ+1+s+22x-γ-1-sx⩾0.
Ands1(x,y)<0 for x∈K1={x∣sθ+(2x-γ-1-sx)[θ+1+s+2(2x-γ-1-sx)]⩾0} and y∈D1=(y¯1,y¯2).
Lets2(x,y)=q(x,y)-1=h(A2y2+B2y+C2) where (62)A2=2hsθ,B2=hs2x+sx-sθ1-x-θ+1+ssx-δ-γ-θ+1+s,C2=1-2x-δ+sx-δ-γ+h1-2x-δsx-δ-γ.Under the conditions (63)x∈K2=x∣2+sδ+2γ-s>0,x>δ+γs,x<1-δ2or2+sδ+2γ-s<0,x>1-δ2,x<δ+γsand Δ2=B22-4A2C2⩾0, the equation s2(x,y)=0 has one real root with multiplicity 2 or two real roots denoted as y¯3=-B2-B22-4A2C2/2A2 and y¯4=-B2+B22-4A2C2/2A2.
Ands2(x,y)>0 for all x∈K2 and y∈D2=(-∞,y¯3)∪(y¯4,+∞).Lemma 8.
If the conditions (61) and (63) hold, and the y¯1<y¯3 or y¯4<y¯2, then D1∩D2≠⌀.
Moreover, if one of above conditions holds and they-coordinate of fixed point Z0(x0,y0), y0, satisfies y0∈Dy0∗⊂D1∩D2, then Z0(x0,y0) is expanding fixed point of (4) in UZ0={(x,y)|x∈K1∩K2,y∈Dy0∗}.
Due to Definition  2 of snap-back repeller, we need to find one pointZ∈Br(Z0) such that Z≠Z0, FM(Z)=Z0, DF(Z)≠0 for some positive integer M.
In fact, we have(64)x1+hx1+θy11-x1-y1-δx1-sx1y1=x2,y1+hsx1y1-δ+γy1=y2,and(65)x2+hx2+θy21-x2-y2-δx2-sx2y2=x0,y2+hsx2y2-δ+γy2=y0.
Now aF2 map has been constructed to map the point Z1(x1,y1) to the fixed point Z0(x0,y0) after two iterations if there are solutions different from Z0 for (64) and (65). By the calculation, the solutions different from Z0 for (65) satisfy the following equation:(66)E1x24+E2x23+E3x22+E4x2+E5=0y2=y01+hsx1-δ-γ,where(67)E1=h3s2,E2=2h2s1-hδ-hγ-h3s2,E3=h1-hδ-hγ2+h2sy0-h2sθy0+h2sδ+h2s2y0-2h2s1-hδ-hγ,E4=hhδ+hγ-11-hδ-hγ-y0-θy0-δ-sy0-h2sθy0-1,E5=hθy01-hδ-hγ-hθy02-x0.
Substitutingx2 and y2 into (64) and solving x1,y1, one gets (68)E1′x14+E2′x13+E3′x12+E4′x1+E5′=0y1=y21+hsx1-δ-γwhere(69)E1′=h3s2,E2′=2h2s1-hδ-hγ-h3s2,E3′=h1-hδ-hγ2+h2sy2-h2sθy2+h2sδ+h2s2y2-2h2s1-hδ-hγ,E4′=hhδ+hγ-11-hδ-hγ-y2-θy2-δ-sy2-h2sθy2-1,E5′=hθy21-hδ-hγ-hθy22-x2.
Obviously, if the condition in Lemma8 is satisfied, the solutions of (64) and (65) satisfies (66) and (68). And the condition of Z1(x1,y1),Z2(x2,y2)≠Z0(x0,y0),Z1(x1,y1)∈UZ0, and DF2(Z1)≠0 holds; then Z0 is a snap-back repeller in UZ0. Thus, the following theorem is established.Theorem 9.
Assume thatZ0(x0,y0) is an expanding fixed point in UZ0, if the solutions (x1,y1),(x2,y2) of (64) and (65) satisfy (x1,y1),(x2,y2)≠(x0,y0),(x1,y1)∈UZ0 and DF2(Z1)≠0. Then Z0(x0,y0) is a snap-back repeller of map (4), and hence map (4) is chaotic in the sense of Marotto.Next, we give specific values of the parameters for illustrating the existence of conditions in Theorem9.Example 10.
Forθ=0.2, s=2.25, δ=0.005, γ=0.465, h=2.5, system (4) has one fixed point Z0(x0,y0)=(0.376,0.2981). Based on Lemma 8 and Theorem 9, we find that region of Z0 is UZ0={(x,y)|0.4975<x<2, y¯1<y<y¯3}⫋{(x,y)|x∈K1∩K2, y∈D1∩D2}, and a point Z1(x1,y1)=(1.19535,0.0309312) such that F2(Z1)=Z0 and DF2(Z1)≠0, where (70)y¯1=1.725x+8.3085-1.725x+8.30852-33.211.465-4.25x216.605y¯3=5.2125-21.725x-6.76354.5.By simple calculation, we can get Z0(x0,y0),Z1(x1,y1)∈UZ0. Thus, Z0(x0,y0) is a snap-back repeller.


## 5. Numerical Simulations

In the section, we use Maximum Lyapunov exponents, the bifurcation diagrams, and phase portraits for system (4) to confirm the above theoretical analysis and show the new interesting complex dynamical behaviors. The bifurcation parameters are considered in the following three cases:(i)
varyingh in the range 2.7⩽h⩽3.5, and fixing θ=0.6, s=0.2, δ=0.0011, γ=0.103,(ii)
varyingh in the range 1.5⩽h⩽2.6, and fixing θ=0.2, s=1.25, δ=0.005, γ=0.465,(iii)
varyingh in the range 1.5⩽h⩽1.57, and fixing θ=0.2, s=4.17, δ=0.033, γ=0.31.For Case (i). Let θ=0.6, s=0.2, δ=0.0011, and γ=0.103, and then system (4) has only one positive fixed point (x∗,y∗)=(0.5205,0.4220). After a simple calculation, we obtain α1=-0.2195≠0 and α2=-13.4227<0. We can discover that system (4) with the above coefficients (θ,s,δ,γ,h)=(0.6,0.2,0.0011,0.103,2.857)∈F1 satisfies all the conditions of the Theorem 3. Figure 1 shows the bifurcation diagram in the h-x plane with the parameters given by case (i), from which we observe that the fixed point (0.5205,0.4220) is stable for h<2.857 and loses its stability at h=2.857. We also observe that there is cascade of period-doubling. The maximum Lyapunov exponents corresponding to bifurcation diagram are computed. Figure 2 shows the phase portrait of period-1, 2, 4, and 8 orbits and chaotic sets for different values of h. The maximum Lyapunov exponents corresponding to h=3.35 are larger than 0 that confirm the existence of the chaotic sets.Figure 1
Lyapunov exponents and bifurcation diagram in theh-x plane with the parameters given by case (i).
(a) (b)Figure 2
Phase portraits for various values ofh in Figure 1.
(a)
h = 2.7 (b)
h = 2.9 (c)
h = 3.2 (d)
h = 3.33 (e)
h = 3.34 (f)
h = 3.35For Case (ii). System (4) has only one positive fixed point (0.376,0.2981) when θ=0.2, s=1.25, δ=0.005, and γ=0.465. After calculation of the positive fixed point of system (4), the Neimark-Sacker bifurcation emerges at h=1.5569 and its eigenvalues are λ,λ¯=0.6199±i0.7847. For h=1.5569, we have |λ,λ¯|=1,l=0.24415>0, b=0.5998. This shows the correctness of Theorem 4. From Figure 3, we observe that the fixed point (0.376,0.2981) of system (4) is stable for h<1.5569, that it loses its stability at h=1.5569, and that an invariant circle appears when the parameter h exceeds 1.5569. The maximum Lyapunov exponents corresponding to Figure 3 are calculated and plotted in these. For h>1.5669, some Lyapunov exponents are bigger than 0 and some are smaller than 0, so there exist stable fixed point or stable period windows in the chaotic region. In general the positive Lyapunov exponent is considered to be one of the characteristics implying the existence of chaos [37–40]. The phase portraits which are associated with Figure 3 are displayed in Figure 4, which clearly depicts how a smooth invariant circle bifurcates from the stable fixed point (0.376,0.2981) and periodic orbits and attracting chaotic sets.Figure 3
Lyapunov exponents and bifurcation diagram in theh-x space for h=1.5–2.6 with the parameters given by case (ii).
(a) (b)Figure 4
Phase portraits for various values ofh in Figure 3.
(a)
h = 1.54 (b)
h = 1.5555 (c)
h = 1.56 (d)
h = 1.9 (e)
h = 2 (f)
h = 2.1 (g)
h = 2.2 (h)
h = 2.25 (i)
h = 2.26 (j)
h = 2.31 (k)
h = 2.32 (l)
h = 2.4For Case (iii). Vary h in the range 1.5⩽h⩽1.57, and fix θ=0.2, s=4.17, δ=0.033, and γ=0.31. We know that the Neimark-Sacker bifurcation emerges only from (x∗,y∗)=(0.0823,0.2494) at h=1.5107, and its eigenvalues are λ1,2=-0.5945±i0.8041. For h=1.5107, we have |λ1,2|=1, l=0.2684>0, and b=-3.3752. This shows the correctness of Theorem 4. From Figure 5, we observe that the fixed point (x∗,y∗)=(0.0823,0.2494) of system (4) is stable for h<1.5107, that it loses its stability at h=1.5107, and that an invariant circle appears when the parameter h exceeds 1.5107. The maximum Lyapunov exponents corresponding to Figure 5 are calculated and plotted in these. Figure 6 shows phase portraits in the s-x plane with the parameters given by case (iii).Figure 5
Lyapunov exponents and bifurcation diagram in theh-x plane for h=1.5–1.57 with the parameters given by case (iii).
(a) (b)Figure 6
Phase portraits for various values ofh in Figure 5.
(a)
h = 1.5 (b)
h = 1.512 (c)
h = 1.513 (d)
h = 1.52


## 6. Chaos Control

In this section, we apply the state feedback control method [24, 41, 42] to stabilize chaotic orbits at an unstable fixed point of (4).Consider the following controlled form of system (4): (71)F:xy⟶x+hx+θy1-x-y-δx-sxy+uy+hsxy-δ+γywith the following feedback control law as the control force: (72)u=-k1x-x∗-k2y-y∗,where k1 and k2 are the feedback gain and (x∗,y∗) is the positive point of model (4). The Jacobian matrix J of the controlled system (71) and (72) evaluated at the point (x∗,y∗) is given by(73)Jx,y=a11-k1a12-k2a21a21,where a11,a12,a21,a22 are given in (26). The characteristic equation of the Jacobian matrix is (74)λ2-a11+a22-k1λ+a22a11-k1-a21a12-k2=0.Assume that the eigenvalues are given by λ1 and λ2; then (75)λ1+λ2=a11+a22-k1and (76)λ1λ2=a22a11-k1-a21a12-k2.The lines of marginal stability are determined by solving the equationλ1=±1 and λ1λ2=1. These conditions guarantee that the eigenvalues λ1 and λ2 have modulus less than 1.Assume thatλ1λ2=1, and from (76) we have (77)l1:k1a22-k2a21=a11a22-a12a21-1.Assume thatλ1=1, and from (75) and (76) we have (78)l2:k11-a22-k2a21=a11+a22-1-a11a22+a12a21.Assume thatλ1=-1, and from (76) we have (79)l3:k11+a22-k2a21=a11+a22+1+a11a22-a12a21.The stable eigenvalues lie within a triangular region by linesl1, l2, and l3 (see Figure 7).Figure 7
The bounded region for the eigenvalues of the controlled system (48), (50) in the (k1,k2) plane for θ=0.2, s=1.25, δ=0.005, γ=0.465, and h=2.3.We have performed some numerical simulations to see how the state feedback method controls the unstable fixed point. Parameter values are fixed asθ=0.2, s=1.25, δ=0.005, γ=0.465, and h=2.3. The initial value (0.3,0.1), and the feedback gain k1=0.91,k2=1.5. It is shown in Figure 8 that a chaotic trajectory is stabilized at the fixed point (x∗,y∗)=(0.376,0.2981).Figure 8
(a) The time responses for the statex of the controlled system (48), (50) in the (n,x) plane for θ=0.2, s=1.25, δ=0.005, γ=0.465, h=2.3, k1=0.91, and k2=1.5. The initial value is (0.3, 0.1). (b) The time responses for the state y of the controlled system (48), (50) in the (n,y) plane.
(a) (b)


## 7. Conclusion

In this paper, we mainly investigated the existence and stability of the nonnegative fixed points of system (4) and flip bifurcation and Neimark-Sacker bifurcation under some conditions at the unique positive fixed point (x∗,y∗) by using manifold theorem and bifurcation theory and also proved the chaos in the sense of Marotto. Our main results are given in Propositions 1, 2 and Theorems 3, 4, and 9 and numerical simulations in Section 5. In the details of the result, when the integral step size h is chosen as a bifurcation parameter, the discrete-time parasite-host system exists with much richer nonlinear dynamical behaviors. Seen from Figures 1–6, there exist period-1, 2, 4, and 8 orbits, attracting invariant cycles, and even stranger chaotic attractors. These results demonstrate that the integral step size h makes differences corresponding to the local and global stability of the discrete-time parasite-host system. Moreover, Figure 1 shows that the system exists in flip bifurcation, which proves the correctness of Theorem 3. In addition, Figure 1 demonstrates the system exists with other complex dynamics including stable coexistence, period-doubling bifurcations. Figures 3 and 5 show that the system exists with Neimark-Sacker bifurcation, which proves the correctness of Theorem 4. The phase portraits corresponding to Figures 1, 3, and 5 are shown in Figures 2, 4, and 6. These figures show that, along with step h increase, the stability properties of the equilibrium will be lost and the flip bifurcation, Neimark-Sacker bifurcation, and chaos phenomenon will appear. Specifically, we have stabilized the chaotic orbits at an unstable fixed point using the feedback control method.However, in this paper we only present the numerical results. In our future work, in order to investigate the more biological significance of system (4), we expect to obtain some real data to illustrate the validity of our theoretical results and to study the effect of bifurcations and chaos on a discrete parasite-host model.



---

*Source: 9275474-2017-01-23.xml*