# Complex Dynamics of a Discrete-Time Predator-Prey System with Ivlev Functional Response

Authors: Hunki Baek
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

The dynamics of a discrete-time predator-prey system with Ivlev functional response is investigated in this paper. The conditions of existence for flip bifurcation and Hopf bifurcation in the interior ofR+2 are derived by using the center manifold theorem and bifurcation theory. Numerical simulations are presented not only to substantiate our theoretical results but also to illustrate the complex dynamical behaviors of the system such as attracting invariant circles, periodic-doubling bifurcation leading to chaos, and periodic-halving phenomena. In addition, the maximum Lyapunov exponents are numerically calculated to confirm the dynamical complexity of the system. Finally, we compare the system to discrete systems with Holling-type functional response with respect to dynamical behaviors.


---

## Body

## 1. Introduction

In population dynamics, continuous-time dynamical systems, such as a Lotka-Volterra system, have been used to understand the interaction between ecological species [1–8]. In recent years, discrete-time population systems also come to the fore due to the following reasons [9–16]: Firstly, discrete-time systems are more suitable than continuous-time systems to describe populations with nonoverlapping generations. Secondly, they can produce more complex and rich dynamical behaviors than continuous-time systems. Finally, we can get more accurate numerical simulation results from discrete-time systems compared with continuous-time systems.In order to investigate the dynamical relationship between two species, predator and prey, one of important components we should pay attention to is the predator’s rate of feeding on the prey, which is called functional response. In this context, Holling [17] suggested three different kinds of functional response for different kinds of species to model the phenomena of predation, which are called Holling’s type I, II, and III. These functional responses are not only monotonically increasing, but also uniformly bounded functions in the first quadratic. The authors [18, 19] introduced another Holling-type functional response called Holling-type IV functional response or Monod-Haldane functional response. Ivlev [20] also suggested a well-known functional response, called Ivlev functional response, in the following form;(1)px=h1-exp⁡-cx,where h, c are positive constants and represent the maximum rate of predation and the decrease in motivation to hunt, respectively.There are many works to investigate predator-prey system with Ivlev-type functional response. The results indicated that predator-prey systems with Ivlev-type have wide applicabilities in ecology, such as dynamics in predator-prey system [21–25], host-parasitoid system [26, 27], and animal coat pattern [28].Taking into account this functional response, we will consider the following predator-prey system with Ivlev functional response:(2)dxdt=ax1-xb-h1-exp⁡-cxy,dydt=-Dy+k1-exp⁡-cxy,where x(t) and y(t) are functions of time representing the population densities of the prey and the predator, respectively. All constants are positive. The parameter a means the growth rate of the prey and b is the carrying capacity. The constant D denotes the death rate of the predator, and k is the rate of conversion of a consumed prey to a predator.Now, we will apply the forward Euler scheme to system (2) to obtain the following discrete-time predator-prey system with Ivlev functional response:(3)x⟶x+δax1-xb-h1-exp⁡-cxy,y⟶y+δ-Dy+k1-exp⁡-cxy,where δ>0 is the time step size.Until now, many researchers have concerned themselves with discrete-time predator-prey systems with Holling-type functional response ([9, 12–16, 29, 30]). Most of these works have dealt with the existence and stability of fixed points and bifurcation phenomena such as flip bifurcation, saddle-node bifurcation, or Hopf bifurcation. On the contrary, in this paper, we will take into account a discrete predator-prey system with Ivlev functional response which has not been studied well until now.The purpose of this paper is to analyze dynamical behaviors of system (3) and to observe a flip bifurcation and Hopf bifurcation by using the center manifold theorem and bifurcation theory. Moreover, numerical simulations are presented not only to substantiate our theoretical results but also to illustrate the complex dynamics of system (3) which are difficult to explain theoretically.This paper is organized as follows. In Section2, we investigate the existence and local stability of the fixed points of system (3). In Section 3, we show that there exist some values of parameters such that system (3) undergoes a flip bifurcation and Hopf bifurcation. In Section 4, we display the numerical simulation to provide numerical evidence for theoretical results and to exhibit numerically various dynamical behaviors including periodic doubling and chaotic phenomena. In addition, we calculate the maximum Lyapunov exponents to confirm the dynamical complexity of system (3). In Section 5, we summarize the results of this paper and compare system (3) to the discrete systems with Holling-type functional response with respect to dynamical behaviors.


## 2. The Existence and Stability of the Fixed Points

In this section, from the point of view of biology, we will focus on the dynamical behaviors of system (3) in the closed first quadrant R+2. In order to determine the fixed points of system (3) we have to solve the nonlinear system given by(4)x=x+δax1-xb-h1-exp⁡-cxy,y=y+δ-Dy+k1-exp⁡-cxy.Simple calculation yields that there are at most three fixed points:(i)
E0(0,0) is the origin;(ii)
E1(b,0) is the axial fixed point in the absence of the predator;(iii)
E2(x∗,y∗) is the interior fixed point, where(5)x∗=1clog⁡kk-D,y∗=akhcD1-1bclog⁡kk-Dlog⁡kk-Dexists if and only if the following condition is satisfied:(6)D<k1-exp⁡-bc.We will study the stability of these fixed points by computing the variation matrix corresponding to each fixed point. The Jacobian matrix of system (3) at the state variable (x,y) is given by(7)Jx,y=1+aδ-2aδxb-chδyexp⁡-cxhδexp⁡-cx-1ckδyexp⁡-cx1-Dδ+kδ1-exp⁡-cx.The characteristic equation of the Jacobian matrix J(x,y) can be obtained as(8)λ2-px,yλ+qx,y=0,where p(x,y) is the trace and q(x,y) is the determinant of the Jacobian matrix J(x,y) which are defined as(9)px,y=2+δa-D+k-2axb-δchy+kexp⁡-cx,qx,y=chkδ2y1-exp⁡-cxexp⁡-cx+1-Dδ+kδ1-exp⁡-cx1+aδb-2xb-chδyexp⁡-cx.In order to determine the stability of the fixed points of system (3), we give the following lemma which can be easily proved by the relations between roots and coefficients of the quadratic equation (8) [9, 30].Lemma 1.
LetF(λ)=λ2-Bλ+C. Suppose that F(1)>0 and λ1 and λ2 are the two roots of F(λ)=0. Then(1)
λ1<1 and λ2<1 if and only if F(-1)>0 and C<1;(2)
λ1<1 and λ2>1 (or λ1>1 and λ2<1) if and only if F(-1)<0;(3)
λ1>1 and λ2>1 if and only if F(-1)>0 and C>1;(4)
λ1=-1 and λ2≠1 if and only if F(-1)=0 and B≠-2,0;(5)
λ1 and λ2 are complex and λ1=λ2=1 if and only if B2-4C<0 and C=1.Letλ1 and λ2 be the two roots of (8), which are called eigenvalues of the fixed point (x,y). We recall some definitions of topological types for a fixed point (x,y) [12, 14, 32]. A fixed point (x,y) is called asink if λ1<1 and λ2<1, so the sink is locally asymptotically stable. (x,y) is called asource if λ1>1 and λ2>1, so the source is locally unstable. (x,y) is called asaddle if λ1<1 and λ2>1 (or λ1>1 and λ2<1). And (x,y) is callednonhyperbolic if either λ1=1 or λ2=1.Theorem 2.
For the fixed pointE0(0,0), one has the following topological types:(1)
E0 is a source if δ>2/D;(2)
E0 is a saddle if 0<δ<2/D;(3)
E0 is nonhyperbolic if δ=2/D.Proof.
It is easy to get that the Jacobian matrixJ at E0 is given by(10)JE0=1+aδ001-Dδ.Hence the eigenvalues of the matrix are 1+aδ and 1-Dδ. Therefore, it is easy to show that (1)–(3) hold.Theorem 3.
For the fixed pointE1(b,0), one has the following topological types:(1)
D>k(1-exp(-bc)):
(i)
E1 is a sink if δ<min⁡{2/a,2/(D+k(exp⁡(-bc)-1))},(ii)
E1 is a source if δ>max⁡{2/a,2/(D+k(exp⁡(-bc)-1))},(iii)
E1 is nonhyperbolic if either δ=2/a or δ=2/(D+k(exp⁡(-bc)-1)).(2)
D<k(1-exp⁡(-bc)):
(i)
E1 is a saddle if δ<2/a,(ii)
E1 is a source if δ>2/a,(iii)
E1 is nonhyperbolic if δ=2/a.(3)
D=k(1-exp⁡(-bc)): E1 is always nonhyperbolic.Proof.
The Jacobian matrixJ at E1 can be obtained as(11)JE1=1-aδhδexp⁡-bc-101-Dδ+kδ1-exp⁡-bc.Then we have the eigenvalues of the characteristic equation (8) as follows:(12)λ1=1-aδ,λ2=1-Dδ+kδ1-exp⁡-bc.Therefore elementary calculation yields that the theorem holds.It follows from Theorem3 that for the fixed point E1(b,0), if (a,b,c,D,k,δ)∈F1, where(13)F1=a,b,c,D,k,δ:δ=2a,D≠k1-exp-bc,δ≠2D+kexp⁡-bc-1,b>0,c>0,D>0,k>0,δ>0,then one of the eigenvalues of the fixed point E1(b,0) is -1 and the other is neither -1 nor 1. Thus system (3) can undergo a flip bifurcation due to the stability of the fixed point E1(b,0) when parameters vary in the small neighborhood of F1. In fact, if the parameters of system (3) are in F1, then a center manifold of system (3) at E1(b,0) is y=0, and if system (3) restricted to this center manifold is the same as the logistic system, dx/dt=ax(1-x/b). Thus in this case the predator becomes extinction; on the other hand, the prey undergoes the period-doubling bifurcation to chaos in the sense of Li-Yorke for suitable bifurcation parameter δ.Theorem 4.
IfD<k(1-exp⁡(-bc)) and bc-log⁡k/k-DD-k-Dlog⁡k/k-D<Dlog⁡(k/(k-D)) are satisfied, then system (3) has a unique positive fixed point E2(x∗,y∗). Moreover, one has the following topological types for the positive fixed point E2:(1)
E2 is a sink if one of the following conditions is satisfied:
(i)
-2Δ1≤Δ2 and δ<-Δ2/Δ1,(ii)
Δ2<-2Δ1 and δ<(-Δ2-Δ22-4Δ1)/Δ1.(2)
E2 is a source if one of the following conditions is satisfied:
(i)
-2Δ1≤Δ2 and δ>-Δ2/Δ1,(ii)
Δ2<-2Δ1 and δ>(-Δ2+Δ22-4Δ1)/Δ1.(3)
E2 is a saddle if Δ2<-2Δ1 and (-Δ2-Δ22-4Δ1)/Δ1<δ<(-Δ2+Δ22-4Δ1)/Δ1.(4)
E2 is nonhyperbolic if one of the following conditions is satisfied:
(i)
-2Δ1<Δ2 and δ=-Δ2/Δ1,(ii)
Δ2<-2Δ1 and δ=(-Δ2±Δ22-4Δ1)/Δ1.
HereΔ1=ak-D/bcbc-log⁡k/k-Dlog⁡(k/(k-D)) and Δ2=(a/bcD)((k-D)log⁡k/k-D2+bc-2D-bcklog⁡(k/(k-D))+bcD).Proof.
LetΔ1=ak-D/bcbc-log⁡k/k-Dlog⁡(k/(k-D)) and Δ2=(a/bcD)((k-D)log⁡k/k-D2+bc-2D-bcklog⁡(k/(k-D))+bcD). Since D<k(1-exp⁡(-bc)), we have log(k/(k-D))<bc, which implies that Δ1>0. Also, simple calculations yield that Δ2=(a/bcD)(bc-log⁡k/k-DD-k-Dlog⁡k/k-D-Dlog⁡(k/(k-D))). Thus Δ2<0 by the hypotheses. From elementary calculation, we obtain that p(x∗,y∗)=Δ2δ+2 and q(x∗,y∗)=Δ1δ2+Δ2δ+1 in (8). Now consider F(λ)=λ2-Bλ+C as in Lemma 1, where B=Δ2δ+2 and C=Δ1δ2+Δ2δ+1. Then we can easily see that(14)F1=ak-Dδ2bcbc-log⁡kk-Dlog⁡kk-D=Δ1δ2>0,F-1=Δ1δ2+2Δ2δ+4.It follows form Lemma 1 that one of roots of F(λ)=0 is -1 and the other is neither -1 nor 1 if and only if Δ2<-2Δ1 and δ=(-Δ2±Δ22-4Δ1)/Δ1, and the roots of F(λ)=0 are a pair of conjugate complex numbers with modulus 1 if and only if -2Δ1<Δ2 and δ=-Δ2/Δ1. Thus, it is easy to see that the theorem holds.Now consider the following sets;(15)F2=a,b,c,D,k,δ∣δ=-Δ2-Δ22-4Δ1Δ1,Δ2<-2Δ1,a>0,b>0,c>0,D>0,k>0,D<k1-exp⁡-bc,bc-log⁡kk-DD-k-Dlog⁡kk-D<Dlog⁡kk-D,F3=a,b,c,D,k,δ∣δ=-Δ2+Δ22-4Δ1Δ1,Δ2<-2Δ1,a>0,b>0,c>0,D>0,k>0,D<k1-exp⁡-bc,bc-log⁡kk-DD-k-Dlog⁡kk-D<Dlog⁡kk-D,FH=a,b,c,D,k,δ∣δ=-Δ2Δ1,Δ2>-2Δ1,a>0,b>0,c>0,D>0,k>0,D<k1-exp⁡-bc,bc-log⁡kk-DD-k-Dlog⁡kk-D<Dlog⁡kk-D.From Theorem 4  (4), we see that if (a,b,c,D,k,δ)∈F2(orF3) then one of the eigenvalues of the Jacobian matrix J(E2) at the fixed point E2(x∗,y∗) is -1 and the other is neither -1 nor 1, and if (a,b,c,D,k,δ)∈FH, then the eigenvalues of the matrix J(E2) are a pair of conjugate complex numbers with modulus 1.In the following section, we will investigate a flip bifurcation and Hopf bifurcation of the fixed pointE2(x∗,y∗) when parameters vary in the small neighborhood of F2 (or F3) and FH, respectively.


## 3. Bifurcation Analysis

In this section, we will discuss two kinds of bifurcation phenomena of system (3), flip bifurcation and Hopf bifurcation, according to the parameter δ. System (3) has no saddle-node bifurcation since the system has at most one positive fixed point E2(x∗,y∗) as mentioned in Section 2.First we investigate a flip (or period-doubling) bifurcation of system (3) at the fixed point E2(x∗,y∗) when parameters vary in a small neighborhood of F2. Similar arguments can be applied to the case F3.Taking parameters(a,b,c,D,k,δ1)∈F2, where δ1=(-Δ2-Δ22-4Δ1)/Δ1, Δ1=ak-D/bcbc-log⁡k/(k-D)log⁡(k/(k-D)), and Δ2=(a/bcD)((k-D)log⁡k/k-D2+bc-2D-bcklog⁡(k/(k-D))+bcD), we consider system (3) with (a,b,c,D,k,δ1) which is described by(16)x⟶x+δ1ax1-xb-h1-exp⁡-cxy,y⟶y+δ1-Dy+k1-exp⁡-cxy.Then system (16) has a unique fixed point E2(x∗,y∗) whose eigenvalues are λ1=-1 and λ2 with λ2≠1.Considering a perturbationδ∗ at δ1, a perturbed system of system (16) is obtained as follows:(17)x⟶x+δ1+δ∗ax1-xb-h1-exp⁡-cxy≡fx,y,δ∗,y⟶y+δ1+δ∗-Dy+k1-exp⁡-cxy≡gx,y,δ∗,where δ∗≪1, which is a small perturbation parameter.In order to transform the fixed pointE2(x∗,y∗) of system (17) into the origin, let u=x-x∗ and v=y-y∗, and then, expanding f and g as a Taylor series at (u,v,δ∗)=(0,0,0) to the second order, system (17) becomes(18)uv⟶a11u+a12v+a13uv+a14u2+e1δ∗+e2uδ∗+e3vδ∗+Ou+v+δ∗3a21u+a22v+a23uv+a24u2+e4δ∗+e5uδ∗+e6vδ∗+Ou+v+δ∗3,where(19)a11=1+aδ1-2aδ1x∗b-cδ1hy∗exp⁡-cx∗,a12=δ1hexp⁡-cx∗,a13=-cδ1hexp⁡-cx∗,a14=-aδ1b+12c2δ1hy∗exp⁡-cx∗,a21=cδ1ky∗exp⁡-cx∗,a22=1-Dδ1-δ1kexp⁡-cx∗,a23=ckδ1exp⁡-cx∗,a24=-12c2δ1ky∗exp⁡-cx∗,e1=-h+ax∗-ax∗2b+hy∗exp⁡-cx∗,e2=a-2ax∗b-chy∗exp⁡-cx∗,e3=hexp⁡-cx∗,e4=k-Dy∗-ky∗exp⁡-cx∗,e5=cky∗exp⁡-cx∗,e6=-D-kexp⁡-cx∗.We construct an invertible matrix(20)T=a12a12-1-a11λ2-a11and using the transformation(21)uv=Tx¯y¯for system (18), then system (18) can be written as(22)x¯y¯=-100λ2x¯y¯+f1x¯,y¯,δ∗g1x¯,y¯,δ∗,where(23)f1x¯,y¯,δ∗=λ2-a11a14-a12a24a121+λ2u2+λ2-a11a13-a12a23a121+λ2uv+λ2-a11e2-a12e5a121+λ2uδ∗+λ2-a11e3-a12e6a121+λ2vδ∗+λ2-a11e1-a12e4a121+λ2δ∗+Ou+v+δ∗3,g1x¯,y¯,δ∗=a11+1a14+a12a24a121+λ2u2+a11+1a13+a12a23a121+λ2uv+a11+1e2+a12e5a121+λ2uδ∗+a11+1e3+a12e6a121+λ2vδ∗+a11+1e1+a12e4a121+λ2δ∗+Ou+v+δ∗3.And u=a12(x¯+y¯),v=-(1+a11)x¯+(λ2-a11)y¯.It follows from the center manifold theorem [32] that there exists a center manifold Wc(0,0) of system (18) at the fixed point (0,0) in a small neighborhood of δ∗=0, which can be approximately represented as follows:(24)Wc0,0=x¯,y¯∣y¯=r1x¯2+r2x¯δ∗+r3δ∗2+Ox¯+δ∗3,where(25)r1=1+a112a13-a121+a11a14-1+a11a23+a12a24λ22-1,r2=1+a112e3-a122e5-1+a11a12e2-e6a12λ2+12+21+a112a13-a121+a11a14-1+a11a23+a12a24a11e1+a12e4-e1λ2a12λ2-1λ2+13,r3=-1+a112e3+a122e5+1+a11a12e2-e6a11e1+a12e4-e1λ2a122λ2-1λ2+13+1+a112a13-a121+a11a14-1+a11a23+a12a24a11e1+a12e4-e1λ22a122λ2-1λ2+14.Thus we can obtain the map, which is restricted to the center manifold Wc(0,0):(26)G:x¯⟶p1δ∗+p2x¯2+p3x¯δ∗+p4xδ∗2+p5x2δ∗+p6x¯3+p7δ∗3+Ox¯+δ∗4,where(27)p1=h1,p2=a122h1-a11+1a12h2,p3=a12h3-a11+1h4,p4=2a122r3h1-a11+1a12r3h2+a12r2h3+a12r3h2λ2-a11+r2h4λ2-a11,p5=2a122r2h1-a11+1a12r2h2+r1a12h3+a12r2h2λ2-a11+r1h4λ2-a11,p6=2r1a122h1-r1a11+1a12h2+r1a12h2λ2-a11,p7=r3a12h3,h1=λ2-a11a14-a12a24a121+λ2,h2=λ2-a11a13-a12a23a121+λ2,h3=λ2-a11e2-a12e5a121+λ2,h4=λ2-a11e3-a12e6a121+λ2,h5=λ2-a11e1-a12e4a121+λ2.In order for the map (26) to undergo flip bifurcation, we require that two discriminatory quantities α1 and α2 are not zero, where(28)α1=∂2G∂x¯∂δ∗+12∂G∂2G∂δ∗∂x¯20,0=p1p2,(29)α2=16∂3G∂x¯3+12∂2G∂x¯220,0=p22+p6.From the above analysis and the theorems in [32] we have the following theorem.Theorem 5.
Ifα1≠0 and α2≠0, then system (3) undergoes a flip bifurcation at the fixed point E2(x∗,y∗) when δ varies in a small neighborhood of δ∗. Moreover, if α2>0 (resp. α2<0), then the period-2 points that bifurcate from E2(x∗,y∗) are stable (resp. unstable).Next we discuss Hopf bifurcation of system (3) at E2(x∗,y∗) when parameters of system (3) vary in a small neighborhood of FH. For this, taking parameters (a,b,c,D,k,δ2)∈FH, where δ2=-Δ2/Δ1, Δ1=ak-D/bcbc-log⁡k/k-Dlog⁡(k/(k-D)), and Δ2=(a/bcD)((k-D)log⁡k/k-D2+bc-2D-bcklog⁡(k/(k-D))+bcD), we consider system (3) with (a,b,c,D,k,δ2) which is described by(30)x⟶x+δ2ax1-xb-h1-exp⁡-cxy,y⟶y+δ2-Dy+k1-exp⁡-cxy.Then it follows from Theorem 4 that the roots of the characteristic equation (8) associated with system (30) at E2(x∗,y∗) are a pair of conjugate complex numbers with modulus 1.Now consider a perturbation system of system (30) by choosing the bifurcation parameter δ∗∗ as follows:(31)x⟶x+δ2+δ∗∗ax1-xb-h1-exp⁡-cxy,y⟶y+δ2+δ∗∗-Dy+k1-exp⁡-cxy,where δ∗∗≪1, which is a small perturbation parameter.Letu1=x-x∗ and v1=y-y∗ in system (31) to transform the equilibrium E2(x∗,y∗) into the origin. Thus system (31) can be rewritten as(32)u1v1⟶f2u1,v1,δ∗∗g2u1,v1,δ∗∗,where(33)f2u1,v1,δ∗∗=u1+δ2+δ∗∗bab-u1-x∗u1+x∗-bhv1+y∗+bhv1+y∗exp⁡-cu1+x∗,(34)g2u1,v1,δ∗∗=v1+δ2+δ∗∗v1+y∗-D+k1-exp-cu1+x∗.Note that the characteristic equation associated with the linearization of the map (32) at (u1,v1)=(0,0) is given by(35)λ2-pδ∗∗λ+qδ∗∗=0,where(36)pδ∗∗=2+Δ2δ2+δ∗∗,qδ∗∗=1+Δ1δ2+δ∗∗2+Δ2δ2+δ∗∗,Δ1=ak-Dbcbc-log⁡kk-Dlog⁡kk-D,Δ2=abcDk-Dlogkk-D2+bc-2D-bcklog⁡kk-D+bcD.The roots of the characteristic equation (35) are(37)λ,λ¯=12pδ∗∗±pδ∗∗2-4qδ∗∗and then we have(38)λ=qδ∗∗1/2.Also, it follows from Theorem 4 that λ=λ¯=1 when δ∗∗=0. Since (a,b,c,D,k,δ2)∈FH, we have Δ1δ2=-Δ2. Therefore,(39)l=dλdδ∗∗δ∗∗=0=-12Δ2>0.In addition, it is required that, whenδ∗∗=0, λm,λ¯m≠1 for m=1,2,3,4, which is equivalent to p(0)≠-2,0,1,2. Since p2(0)-4q(0)<0 and q(0)=1, p(0)≠-2,2. Thus we only show that p(0)≠0,1, which leads to p(0)=2+Δ2δ2≠0,1. Then -Δ2δ2≠1,2 and therefore(40)-aδ2bcDk-Dlog⁡kk-D2-δ2bc-2D-bcklog⁡kk-D+bcD≠1,2.Hence the eigenvalues λ, λ¯ do not lie in the intersection of the unit circle with the coordinate axes when δ∗∗=0 and (40) holds.Next we investigate the normal form of system (32) when δ∗∗=0. Let(41)α=1+12δ2abcDk-Dlogkk-D2+bc-2D-bcklogkk-D+bcD,β=12δ24Δ1-Δ22,T=01βα.By expanding f2 and g2 in system (32) as Taylor series at (u1,v1)=(0,0) to the third order and use the transformation(42)u1v1=Tx~y~for system (32), system (32) becomes(43)x~y~=α-ββαx~y~+Px~,y~Qx~,y~,where(44)Px~,y~=αβa~23-αa~13x~y~+-αa~13+a~23+1βa~24-αa~14x~2+αβa~25-αa~15x~2y~+-αa~15+a~25+1βa~26-αa~16x~3+Ox~+y~4,Qx~,y~=αa~13x~y~+a~14+βa~13x~2+αa~15x~2y~+a~16+βa~15x~3+Ox~+y~4,a~13=-cδ2hexp⁡-cx∗,a~14=-aδ2b+12c2δ1hy∗exp⁡-cx∗,a~15=12c2δ2hy∗exp⁡-cx∗,a~16=-16c3δ2hy∗exp⁡-cx∗,a~23=ckδ2∗exp⁡-cx∗,a~24=-12c2δ2ky∗exp⁡-cx∗,a~25=-12c2δ2kexp⁡-cx∗,a~26=16c3δ2ky∗exp⁡-cx∗.In order to undergo Hopf bifurcation, it is required that the following discriminatory quantity A is not zero [33]:(45)A=-Re1-2λλ¯21-λξ11ξ20-12ξ112-ξ022+Reλ¯ξ21,where(46)ξ20=18Px~x~0,0-Py~y~0,0+2Qx~y~0,0+iQx~x~0,0-Qy~y~0,0-2Px~y~0,0=14αa~13-αββa~13+a~14+a~24β+a~23+iβa~13+a~14+αβαa~13-a~23,ξ11=14Px~x~0,0+Py~y~0,0+iQx~x~0,0+Qx~y~0,0=12-αββa~13+a~14+a~24β+a~23+i14α+2βa~13+2a~14,ξ02=18Px~x~0,0-Py~y~0,0-2Qx~y~0,0+iQx~x~0,0-Qy~y~0,0+2Px~y~0,0=14-αa~13-αββa~13+a~14+a~24β+a~23+iβa~13+a~14-αβαa~13-a~23,ξ21=116Px~x~x~0,0+Px~y~y~0,0+Qx~x~y~0,0+Qy~y~y~0,0+iQx~x~x~0,0+Qx~y~y~0,0-Px~x~y~0,0-Py~y~y~0,0=18αa~15-3αββa~15+a~16+3a~25+a~26β+i3βa~15+a~16+αβαa~15-a~25.From the above analysis and theorems in [33], we have the following result.Theorem 6.
If condition (40) holds and A≠0, then system (3) undergoes Hopf bifurcation at the fixed point E2(x∗,y∗) when the parameter δ∗∗ varies in the small neighborhood of the origin. Moreover, if A<0 (resp., A>0), then an attracting (resp., repelling) invariant closed curve bifurcates from E2(x∗,y∗) for δ∗∗>0 (resp., δ∗∗<0).


## 4. Numerical Simulations

In this section, we present some bifurcation diagrams and phase portraits of system (3) to substantiate the above theoretical results and illustrate various dynamical behaviors via numerical simulations. In particular, we calculate the maximum Lyapunov exponents numerically to show the dynamical complexity of system (3). The bifurcation parameters are considered in the following two cases:(i)
Varyingδ in range 3.9≤δ≤4.35, and fixing a=0.8, b=0.3, c=0.5, D=0.1, h=0.3, k=1.(ii)
Varyingδ in range 2.35≤δ≤4.15, and fixing a=0.1, b=1, c=1, D=0.4, h=0.6, k=1.For case (i), based on Theorem4, we can know that system (4) has only one positive fixed point E2(x∗,y∗)=(0.2107,1.6723) since D<k(1-exp⁡(-bc))=0.1393 and bc-log⁡k/k-DD-k-Dlog⁡k/k-D-Dlog⁡(k/(k-D))=-0.0103<0. It is obtained from tedious calculation that α1=p1p2=121.9182 and α2=p22+p6=21.2815 in (28) and (29), respectively. Thus it follows from Theorem 5 that a flip bifurcation emerges from the fixed point E2 at δ1=3.9613 when (a,b,c,D,k,δ1)=(0.8,0.3,0.5,0.1,1,3.9613) in F2. The bifurcation diagrams of system (3) are given in Figures 1(a) and 1(b) for showing this phenomenon. From Figure 1, we can observe that the fixed point E2 of system (3) is stable for δ<3.9613 and loses its stability at δ=3.9613. We can also see that there is a cascade of period-doubling due to the flip bifurcation. In fact, there are orbits of period 1,2,4,8 for δ∈(3.9,4.27), which is illustrated in Figure 3. Moreover, chaotic phase portraits can be observed as shown in Figure 4. In order to confirm the existence of chaotic trajectories of system (3), the maximum Lyapunov exponents corresponding to Figure 1 are computed in Figure 2.Figure 1
(a) Bifurcation diagram of system (3) in the (δ,x) plane when a=0.8, b=0.3, c=0.5, D=0.1, h=0.3, k=1, and 3.9≤δ≤4.35. The initial value is (0.2117, 1.6733). (b) Bifurcation diagram of system (3) in the (δ,y) plane.
(a)
(b)Figure 2
Maximum Lyapunov exponents of system (3) corresponding to Figure 1.Figure 3
(a) Phase portrait of a period-1 orbit of system (3) when a=0.8, b=0.3, c=0.5, D=0.1, h=0.3, k=1, and δ=3.95. (b) Phase portrait of a period-2 orbit when δ=3.97. (c) Phase portrait of a period-4 orbit when δ=4.15. (d) Phase portrait of a period-8 orbit when δ=4.27.
(a)
δ
=
3.95
(b)
δ
=
3.97
(c)
δ
=
4.15
(d)
δ
=
4.27Figure 4
Chaotic phase portraits of system (3) when a=0.8, b=0.3, c=0.5, D=0.1, h=0.3, k=1. (a) δ=4.33. (b) δ=4.36.
(a)
δ
=
4.33
(b)
δ
=
4.36For case (ii), it follows from Theorem4 that system (4) has only one positive fixed point E2(x∗,y∗)=(0.5108,0.1041), which is sink when δ<δ2=-Δ2/Δ1=2.6444, since D<k(1-exp⁡(-bc))=0.6321, bc-log⁡k/k-DD-k-Dlog⁡k/k-D-Dlog⁡(k/(k-D))=-0.1586<0, and Δ2<-2Δ1. After elementary calculation, we can obtain that -(aδ2/bcD)((k-D)log⁡k/k-D2-δ2bc-2D-bcklog⁡(k/(k-D))+bcD)=-1.6182 in (40) and A=-3.4722 in (45), which means that the hypotheses of Theorem 6 hold. Thus from Theorem 6, Hopf bifurcation emerges from the fixed point (0.5108,0.1041) at δ2=2.6444. In order to investigate various dynamical phenomena for this case, the bifurcation diagrams of system (3) are displayed in Figures 5(a) and 5(c). It is observed from Figure 5(a) that the fixed point (0.5108,0.1041) of system (3) is sink for δ<δ2, and an attracting invariant close curve appears, as shown in Figure 7, when the parameter δ exceeds δ2. After the invariant circle occurs, its diameter is bigger and bigger according the value of δ. In addition, period-doubling and period-halving phenomena are also observed in Figures 5(b) and 5(d) which are the local amplification of Figures 5(a) and 5(c), respectively, for δ∈[4.04,4.09]. In fact, if the value of δ exceeds 4.048, as shown in Figure 8, orbits of period 25 and 50 appear. Also, periodic-halving phenomenon is detected in Figure 9. The maximum Lyapunov exponents corresponding to Figure 5(a) are calculated and plotted in Figure 6 to display the dynamical complexity of system (3).Figure 5
(a) Bifurcation diagram of system (3) in the (δ,x) plane when a=0.1, b=1, c=1, D=0.4, h=0.6, k=1, and 2.35≤δ≤4.15. The initial values is (0.5118,0.1051). (b) Local amplification corresponding to (a) for δ∈[4.04,4.09]. (c) Bifurcation diagram of system (3) in the (δ,y) plane. (d) Local amplification corresponding to (c) for δ∈[4.04,4.09].
(a)
(b)
(c)
(d)Figure 6
(a) Maximum Lyapunov exponents of system (3) corresponding to Figure 5(a). (b) Local amplification corresponding to (a) for δ∈[4.04,4.09].
(a)
2.35
≤
δ
≤
4.15
(b)
4.04
≤
δ
≤
4.09Figure 7
(a) Phase portrait of a stable fixed point of system (3) when a=0.1, b=1, c=1, D=0.4, h=0.6, k=1, and δ=2.6. (b) Phase portrait of system (3) after Hopf bifurcation when δ=2.72. (c) Phase portrait of system (3) after Hopf bifurcation when δ=2.9. Here, the initial point is (0.6108, 0.2041).
(a)
δ
=
2.6
(b)
δ
=
2.72
(c)
δ
=
2.9Figure 8
(a) Phase portrait of a period-25 orbit of system (3) when a=0.1, b=1, c=1, D=0.4, h=0.6, k=1, and δ=4.05. (b) Phase portrait of a period-50 orbit when δ=4.056.
(a)
δ
=
4.05
(b)
δ
=
4.056Figure 9
(a) Phase portrait of a period-50 orbit of system (3) when a=0.1, b=1, c=1, D=0.4, h=0.6, k=1, and δ=4.072. (b) Phase portrait of a period-25 orbit when δ=4.075.
(a)
δ
=
4.072
(b)
δ
=
4.075


## 5. Conclusion and Discussion

In summary, we considered the discrete-time predator-prey system with the step sizeδ obtained from a continuous predator-prey system with Ivlev functional response. Firstly, we found all fixed points of the system and investigated their stability via the eigenvalues of Jacobian matrix. Based on this analysis, we provided the parameter sets for which the system undergoes a flip bifurcation and Hopf bifurcation according to the value δ by using the center manifold theorem and bifurcation theory. In addition, we demonstrated numerical simulations including phase portraits and bifurcation diagrams of the system in order to validate our mathematical results and show that the system has complex dynamical behaviors such as a cascade of periodic doubling, an attracting invariant circle, or chaotic attractors. Also, we computed the Lyapunov exponents to confirm the dynamical complexity of the system.In fact, it follows from Theorem5 that system (3) can have a flip bifurcation, which results in periodic doubling and chaotic phenomena as shown in Figures 1 and 4. Here the chaotic phenomenon indicates biologically that the biological system (3) is not stable and, particularly, if the prey is in chaotic, then the predator will eventually go extinct or go to a stable equilibrium. Also Theorem 6 guarantees the existence of a closed curve as shown in Figure 7, which means biologically that the predator coexists with the prey with at least two-frequency oscillation.There is an interesting thing that the parameterh, the maximum rate of predation, has no effect on qualitative dynamical behaviors including the stability of the fixed points and the existence of bifurcation phenomena. For instant, if we take h=0.1 and h=0.2 instead of h=0.3 and h=0.6 in cases (i) and (ii), respectively, we can obtain qualitatively similar bifurcation diagrams in Figure 10 to Figures 1 and 5.Figure 10
(a) Bifurcation diagram of system (3) in the (δ,x) plane when a=0.1, b=1, c=1, D=0.4, h=0.1, k=1, and 3.9≤δ≤4.35. (b) Bifurcation diagram of system (3) in the (δ,y). (c) Bifurcation diagram of system (3) in the (δ,x) plane when a=0.1, b=1, c=1, D=0.4, h=0.2, k=1, and 1≤δ≤3. (d) Bifurcation diagram of system (3) in the (δ,y).
(a)
(b)
(c)
(d)As mentioned in the introduction, there are many works about discrete-time predator-prey systems with Holling-type functional responses. Now we will compare the dynamics of these systems with that of system (3). If a predator-prey system has stable positive fixed points, then all species of the system could coexist, which is a biologically significant fact from the point of view of biodiversity. In this context, we will focus on positive fixed points. We summarize the dynamical behavior of systems around positive fixed points in Table 1 by letting LAS, SNB, FB, and HB stand for locally asymptotically stable, saddle-node bifurcation, flip bifurcation, and Hopf bifurcation, respectively.Table 1
This table lists the stability and bifurcations of systems according to functional response.
Systems
Existence
LAS
SNB
FB
HB
References
System (2)
Yes
Yes
No
No
No
[23]
System (3)
Yes
Yes
No
Yes
Yes
–
Lotka-Volterra
Yes
Yes
No
Yes
Yes
[15]
Holling II
Yes
Yes
No
Yes
Yes
[29]
Holling III
Yes
Yes
No
Yes
Yes
[13]
Holling IV
Yes
Yes
Yes
Yes
Yes
[31]Comparing the first row to the second row of Table1, we figure out that the discrete-time predator-prey system exhibits much richer dynamics than the continuous-time system. From the second row to the fifth row, we assert that the discrete-time predator-prey systems have similar dynamical behaviors around their positive fixed points even if their functional responses are different from each other. Moreover, we conclude that the discrete-time systems have various dynamical behaviors regardless of the functional responses.It can be inferred from the results obtained in this paper that the Ivlev functional response may be a cause of chaos-like oscillation of populations and can make the behavior of population systems more complex. Nevertheless, identifying complicated, possibly chaotic dynamics in population data has remained a major challenge in ecological studies [34, 35]. However, this work could be very helpful for the biologists who work with discrete-time prey-predator systems.



---

*Source: 8635937-2018-02-08.xml*