# Dynamic Analysis of an Impulsively Controlled Predator-Prey Model with Holling Type IV Functional Response

Authors: Yanzhen Wang; Min Zhao
Journal: Discrete Dynamics in Nature and Society
 (2012)

---

## Abstract

The dynamic behavior of a predator-prey model with Holling type IV functional response is investigated with respect to impulsive control strategies. The model is analyzed to obtain the conditions under which the system is locally asymptotically stable and permanent. Existence of a positive periodic solution of the system and the boundedness of the system is also confirmed. Furthermore, numerical analysis is used to discover the influence of impulsive perturbations. The system is found to exhibit rich dynamics such as symmetry-breaking pitchfork bifurcation, chaos, and nonunique dynamics.


---

## Body

## 1. Introduction

In recent years, impulsive control strategies in predator-prey models have become a major field of inquiry. Many authors have studied the dynamics of predator-prey models with impulsive control strategies [1–10]. Much research has also been done on three-species food chain systems with impulsive perturbations [11–17]. Holling-type functional responses are well known, and therefore many authors have studied the dynamics of predator-prey models with different Holling-type functional responses with respect to an impulsive control strategy. For example, the dynamics of a predator-prey model with Holling type I functional response with respect to an impulsive control strategy have been reported in [18]. The dynamics of a predator-prey model with Holling type II functional response with an impulsive control strategy were presented in [19, 20]. In [21, 22], the results of studies of the dynamics of a predator-prey model with Holling type VI functional response with respect to an impulsive control strategy were presented. The following predator-prey model with Holling type VI functional response with respect to an impulsive control strategy was proposed in [23]:
(1.1)dx(t)dt=x(t)(a-bx(t))-c1x(t)y(t)1+e1(x(t))2,dy(t)dt=c2x(t)y(t)1+e1(x(t))2-d1y(t),t≠(n+τ-1)T,t≠nT,x(t+)=(1-p1)x(t),y(t+)=(1-p2)y(t),t=(n+τ-1)T,x(t+)=x(t),y(t+)=y(t)+q,(x(0+),y(0+))=(x0,y0),t=nT,
where x(t) and y(t) are the respective densities of prey and predator at time t, the constant a is the intrinsic growth rate of the prey population, b is the coefficient of intraspecies competition, c1 is the per-capita rate of predation of the predator, e1 is the half-saturation constant, which depicts a critical concentration of the nutrient composition to maintain the normal growth of the predator, d1 is the death rate of the predator, and c2 is the rate of conversion of a consumed prey to a predator.More recently, the author of [24] has developed the following three-species Holling type IV system by introducing spraying of pesticides and periodic constant release of mid-level predators at different fixed times:(1.2)dx(t)dt=x(t)(a-bx(t))-c1x(t)y(t)1+e1(x(t))2,dy(t)dt=c2x(t)y(t)1+e1(x(t))2-c3y(t)z(t)1+e2(y(t))2-d1y(t),dz(t)dt=c4y(t)z(t)1+e2(y(t))2-d2z(t),t≠(n+τ-1)T,t≠nT,Δx(t)=-p1x(t),Δy(t)=-p2y(t),Δz(t)=-p3z(t),t=(n+τ-1)T,x(t+)=x(t),y(t+)=y(t)+q,z(t+)=z(t)t=nT,(x(0+),y(0+))=(x0,y0),
where c3 is the per-capita rate of predation of the top predator, e2 is the half-saturation constant, d2 is the death rate of the top predator, and c4 is the rate of conversion of consumed mid-level predators to the top predator.Assume that the top predator also eats the prey, or in other words, that the relationship between the top predator and the mid-level predator is not only that of predator and prey, but also that of competitors. To represent this, a predator-prey model with Holling type IV functional response with respect to an impulsive control strategy can be constructed as follows:(1.3)dx(t)dt=x(t)(a-bx(t))-c1x(t)y(t)1+e1(x(t))2-c5x(t)z(t)1+e3(x(t))2,dy(t)dt=c2x(t)y(t)1+e1(x(t))2-c3y(t)z(t)1+e2(y(t))2-d1y(t),dz(t)dt=c4y(t)z(t)1+e2(y(t))2+c6x(t)z(t)1+e3(x(t))2-d2z(t),t≠nT,Δx(t)=0,Δy(t)=-qy(t)+p,Δz(t)=0,t=nT,
where x(t), y(t), and z(t) are the respective densities of the prey, top predator, and mid-level predator at time t, c5 is the per-capita rate of predation of the top predator, e3 is the half-saturation constant, and c6 is the rate of conversion of consumed prey to the top predator. The meanings of the other parameters are the same as in (1.1) and (1.2). Δx(t)=x(t+)-x(t), Δy(t)=y(t+)-y(t), and Δz(t)=z(t+)-z(t), 0≤q<1, represent the fractions of prey and predator which die because of harvesting or for other reasons, p is the magnitude of immigration or stocking of the predator, and T is the period of impulsive immigration or stocking of the predator. Here it is assumed that the rate of conversion of consumed prey to predator is smaller than the per-capita rate of predation of the predator.The rest of this paper is organized as follows. Section2 presents a mathematical analysis of the model. Section 3 describes some numerical simulations. In the last section, a brief discussion is provided, and conclusions are drawn.


## 2. Mathematical Analysis

First, some useful notations and statements will be provided for use in subsequent proofs. The following definitions will be useful.LetR+=[0,+∞),R+3={X=(x(t),y(t),z(t))∈R3∣X≥0}. Denote by f=(f1,f2,f3) the map defined by the right-hand sides of the first, second, third, and fourth equations of system (1.3). LetV0={V:R+×R+3→R+}; then V is continuous on (nT,(n+1)T]×R+3,n∈N, lim(t,μ)→(nT+,X)V(t,μ)=V(nT+,X) exists, and V is locally Lipschitzian in X.Definition 2.1.
LetV∈V0; then for (nT,(n+1)T]×R+3, the upper right derivative of V(t,X) with respect to the impulsive differential system (1.3) can be defined as
(2.1)D+V(t,X)=lim⁡h→0+sup⁡1h[V(t+h,X+hf(t,X))-V(t,X)].
The solution of system (1.3) is a piecewise continuous function X:R+×R+3, X(t) is continuous on (nT,(n+1)T],n∈N, and X(nT+)=limt→nT+X(T) exists. Obviously, the smoothness properties of f guarantee the global existence and uniqueness of a solution of system (1.3); for details, see [25, 26].Definition 2.2.
System (1.3) is permanent if there exists M≥m>0 such that, for any solution (x(t),y(t),z(t)) of system (1.3) with a positive initial value X0, m≤limt→∞infx(t)≤limt→∞supx(t)≤M, m≤limt→∞infy(t)≤limt→∞supy(t)≤M, and m≤limt→∞infz(t)≤limt→∞supz(t)≤M.Lemma 2.3.
Assume thatX(t) is a solution of system (1.3) with X(0+)≥0; then X(t)≥0 for all t≥0. Furthermore, X(t)>0, t>0 if X(0+)>0.Lemma 2.4 (see [23]).
LetV∈V0 and assume that
(2.2)D+V(t,X)≤g(t,V(t,X)),t≠nT,V(t,X(t+))≤Φn(V(t,X(t))),t=nT,
where g:R+×R+→R is continuous on (nT,(n+1)T] for u∈R+2,n∈N, lim(t,y)→(nT+)g(t,v)=g(nT+,u) exists, and ϕni(i=1,2):R+→R+ is nondecreasing. Let r(t) be a maximal solution of the scalar impulsive differential equation:
(2.3)du(t)dt=g(t,u(t)),t≠nT,u(t+)=Φn(u(t)),t=nT,u(0+)=u0
existing on (0,+∞]. Then V(0+,X0)≤u0 implies that V(t,X(t))≤r(t),t≥0, where X(t) is any solution of system (1.3). If certain smoothness conditions on g can be established to guarantee the existence and uniqueness of solutions for (2.3), then r(t) is exactly a unique solution of (2.3).Now consider a special case of Lemma2.4. Let PC(R+,R)(PC1(R+,R)) denote the class of real piecewise continuous functions defined on R+.Lemma 2.5 (see [23]).
Letu(t)∈PC1(R+,R) and let it satisfy the inequalities:
(2.4)u′(t)≤f(t)u(t)+h(t),t≠τk,t>0,u(τk+)≤αku(τk)+βk,k≥0,u(0+)≤u0,
where f,h∈PC1(R+,R) and αk≥0, βk and u0 are constants, and τk(k≥0) is a strictly increasing sequence of positive real numbers. Then, for t>0,
(2.5)u(t)≤u0(∏0<τk<tαk)exp⁡(∫0tf(s)ds)+∫0t(∏s≤τk<tαk)exp⁡(∫stf(γ))h(s)ds+∑0<τk<t(∏τk<τj<tαj)exp⁡(∫τktf(γ)dγ)βk.For convenience, some basic properties can be defined for the following subsystems of system (1.3):(2.6)dy(t)dt=-d1y(t),t≠nT,y(t+)=(1-q)y(t)+p,t=nT,y(0+)=y0.
Then the following lemma results.Lemma 2.6.
For a positive periodic solutiony*(t) of system (2.6) and a solution y(t) of system (2.6) with initial value y0=y(0+)≥0,
(2.7)|y(t)-y*(t)|⟶0,t⟶∞,
where
(2.8)y*(t)=pexp⁡(-d1(t-nT))1-(1-q)exp⁡(-d1T),t∈(nT,(n+1)T],n∈N,y*(0+)=p1-(1-q)exp⁡(-d1T),y(t)=(1-q)n+1(y(0+)-(p1-(1-q)exp⁡(-d1T)))exp⁡⁡(-d1T)+y*(t).Next, some main theorems will be proposed.Theorem 2.7 (boundedness).
There exists a constantM such that x(t)≤M, y(t)≤M, and z(t)≤M for each solution X=(x(t),y(t),z(t)) of system (1.3) for all sufficiently large t.Proof.
Let(x(t),y(t),z(t)) be a solution of system (1.3); let u(t)=x(t)+y(t)+z(t),(t≥0). Then
(2.9)u′(t)=x(t)(a-bx(t))+(c2-c1)x(t)y(t)1+e1(x(t))2+(c4-c3)y(t)z(t)1+e2(y(t))2+(c6-c5)x(t)z(t)1+e3(x(t))2-d1y(t)-d2z(t).
From the first part, it is known that c2-c1<0,c4-c3<0,c6-c5<0, so u′(t)≤x(t)(a-bx(t))-d1y(t)-d2z(t). Choosing 0<C<min{d1,d2}, then u′(t)+Cu(t)≤-b(x(t))2+(a+C)x(t)<β0,(β0=(a+C)2/4b). The following can then be obtained:
(2.10)u′(t)+Cu(t)≤β0,t≠nT,u(nT+)≤(1-q)u(nT)+p.
From Lemma 2.5,
(2.11)u(t)≤u0(1-q)[t/kT]exp⁡(∫0t-Cds)+∫0t(1-q)[(t-s)/kT]exp⁡(∫st-Cdγ)β0ds+∑j=1[t/kT](1-q)[(t-s)/kT]exp⁡(∫kTt-Cdγ)p≤u(0+)exp⁡(-Ct)+β0C(1-exp⁡(-Ct))+pexp⁡(CT)exp⁡(CT)-1.
When t→∞, (2.11) →β0/C+pexp(CT)/(exp(CT)-1)<∞. It is clear that u(t) is bounded for sufficiently large t. Therefore, x(t), y(t), and z(t) are bounded. This completes the proof.Next the stability of a prey and top-predator eradication periodic solution will be examined.Theorem 2.8.
The solution(0,y*(t),0) is locally asymptotically stable if T<c1Γ/a and T>(c4/d2)((θ1-θ2)/d1e2), where Γ=p(1-exp(-d1T))/d1(1-(1-q)exp(-d1T)), θ1=e2p/Λ, θ2=(e2p/Λ)exp(-d1T), and Λ=e2(1-(1-q)exp(-d1T)).Proof.
The local stability of the periodic solution(0,y*(t),0,0) may be determined by considering the behavior of small-amplitude perturbations of the solution. Define
(2.12)x(t)=u(t),y(t)=v(t)+y*(t),z(t)=w(t).
Substituting (2.12) into system (1.3), it is possible to obtain a linearization of the system as follows:
(2.13)du(t)dt=(a-c1y*(t))u(t),dv(t)dt=c2y*(t)u(t)-d1v(t)-c3y*(t)1+e2(y*(t))2w(t)dw(t)dt=(-d2+c4y*(t)1+e2(y*(t))2)w(t),t≠nT,Δu(t)=0,Δv(t)=-qy(t)+pΔw(t)=0t=nT,
which can be written as
(2.14)(u(t)v(t)w(t))=ϕ(t)(u(0)v(0)w(0)),0≤t≤T,
where the ϕ(t) satisfies
(2.15)dϕ(t)dt=(a-c1y*(t)00c2y*(t)-d1-c3y*(t)1+e2(y*(t))200-d2+c4y*(t)1+e2(y*(t))2)
with ϕ(0)=I, where I is the identity matrix and
(2.16)(u(nT+)v(nT+)w(nT+))=(100010001)(u(nT)v(nT)w(nT)).
Hence, the stability of the periodic solution (0,y*(t),0) is determined by the eigenvalues of
(2.17)θ=(100010001)ϕ(t).
If the absolute values of all eigenvalues are less than one, the periodic solution (0,y*(t),0) is locally stable. Then all eigenvalues of ϕ can be denoted by λ1,λ2, and λ3:
(2.18)λ1=exp⁡(∫0T(a-c1y*(t)dt)),λ2=exp⁡(-d1T)<1,λ3=exp⁡(∫0T(-d2+c4y*(t)1+e2(y*(t))2)dt).
Therefore,
(2.19)∫0Ty*(t)dt=p(1-exp⁡(-d1T))d1(1-(1-q)exp⁡(-d1T))=ΔΓ,∫0Ty*(t)1+e2(y*(t))2dt=1d1e2(θ1-θ2),
where θ1=tan-1(e2p/Λ), θ2=tan-1((e2p/Λ)exp(-d1T)), and Λ=e2(1-(1-q)exp(-d1T)). Clearly, |λ2|<1 and |λ1|<1 only if T<c1Γ/a, and |λ3|<1 only if T>(c4/d2)((θ1-θ2)/d1e2). According to the Floquet theory of impulsive differential equations, the periodic solution (0,y*(t),0) is locally stable. This completes the proof.Theorem 2.9.
The solution(0,y*(t),0) is said to be globally stable if T>c4Γ/d2.Proof.
By Theorem2.7, there exists a constant M>0 such that x(t)≤M, y(t)≤M, z(t)≤M for each solution X=(x(t),y(t),z(t)) of system (1.3) for all sufficiently large t. From the first equation in system (1.3), it is possible to obtain
(2.20)x′(t)≤ax(t)-c1x(t)y(t)1+e1M2=(a-c1y(t)1+e1M2)x(t).
Then,
(2.21)dydt=c2x(t)y(t)1+e1(x(t))2-c3y(t)z(t)1+e2(y(t))2-d1y(t)≤-(d1-c2M)y(t),t≠nTΔy(t)=-qy(t)+p,t=nT
and using (2.21),
(2.22)y(t)≥y1*(t)-ε=pexp⁡(-(d1-c2M)(t-nT))1-(1-q)exp⁡(-(d1-c2M)T)-ε.
Denoting λ as: λ=(pexp(-(d1-c2M)(t-nT)))/(1-(1-q)exp(-(d1-c2M)T))-ε, by Lemmas 2.4 and 2.6, there exists a t1>0, and it is possible to select ε>0. If ε is small enough, then y(t)≥y1*(t)-ε for all t≥t1 and satisfying a<c1λ/(1+e1M2). Therefore, (2.20) ≤(a-c1λ/(1+e1M2))x(t); if a<c1λ/(1+e1M2), then (2.20) <0. As t→∞, x(t)→0; this implies that there exist values of ε1 and T1 such that x(t)<ε1 for t≥T1, and when T>c4Γ/d2, these satisfy δ=exp(-d2T+c4Γε1+c4ε1T+c6ε1T)<1. From the second equation in system (1.3), y′(t)≤y(t)(-d1+c2x(t))≤y(t)(-d1+c2ε1). Let y1(t) be the solution of the following equation:
(2.23)y1′(t)=-(d1-c2ε1)y1(t),t≠nT,y1(t+)=(1-q)y1(t)+p,t=nT,y1(0+)=y0.
From Lemmas 2.4 and 2.6, and the third equation in system (1.3),
(2.24)z′(t)≤z(t)(-d2+c4y1(t)+c6x(t))≤z(t)(-d2+c4y1*(t)+c6ε1),
where y1*(t) is the periodic solution of (2.23):
(2.25)y*(t)=pexp⁡(-(d1-c2ε1)(t-nT))1-(1-q)exp⁡(-(d1-c2ε1)T).
Integrating (2.24) over (nT,(n+1)T],
(2.26)z((n+1)T)≤z(nT+)exp⁡∫nT(n+1)T(-d2+c4y1*(t)+c4ε1+c6ε1)dt≤z(nT)δ,
where δ=exp(-d2T+c4Γε1+c4ε1T+c6ε1T) and
(2.27)Γε1=p(exp⁡(Td1)-exp⁡(Tc2ε1))(1-q)(1-d1)exp⁡(Tc2ε1)-(c2ε1-d1)exp⁡(Td1).
When δ<1, z((n+1)T)≤z(0)δn→0 as n→∞. For t∈(nT,(n+1)T],
(2.28)z(t)≤z(nT+)≤z(nT)∫nTt(-d2+c4y1*(t)+c4ε1+c6ε1)dt≤z(0)δn,
which implies that z(t)→0 as t→∞. Therefore, it can be assumed that z(t)<ε2 for t>0. From the second equation in system (1.3), y′(t)≥y(t)(-d1-c3z(t))≥y(t)(-d1-c3ε2). Let y2(t)  and y2*(t) be solutions of the following equation:
(2.29)y2′(t)=-(d1+c3ε2)y2(t),t≠nT,y2(t+)=(1-q)y2(t)+p,t=nT,y2(0+)=y0.
From Lemmas 2.4 and 2.6, y2(t)≤y(t)≤y1(t) and y1(t)→y1*(t), y2(t)→y2*(t) as t→∞. Also note that y1*(t)→y*(t) and y2*(t)→y*(t) as ε1→0 and ε2→0. Therefore, y(t)→y*(t) as t→∞. This completes the proof.Theorem 2.10.
System (1.3) is permanent if
(2.30)T>c1Γa-bM0-c5M0>c1Γa,T<c4d2θ1-θ2d1e2.Proof.
Let(x(t),y(t),z(t)) be a solution of system (1.3). From Theorem 2.7, there exists a constantM>0 such that x(t)≤M, y(t)≤M, z(t)≤M for each solution X=(x(t),y(t),z(t)) of system (1.3) for all sufficiently large t. Let M0=max{M,M/c3}; then x(t)≤M0, y(t)≤M0, z(t)≤M0, and y′(t)≥-(d1+M0)y(t).
By Lemmas2.4 and 2.6, y(t)≥u*(t)-ε,(ε>0), where
(2.31)u*(t)=pexp⁡(-(d1+M0)(t-nT))1-(1-q)exp⁡(-(d1+M0)T).
Then, y(t)≥(pexp(-(d1+M0)))/(1-(1-q)exp(-(d1+M0)T))-ε≡m0, for sufficiently large t. Therefore, it is necessary only to find an m2>0 such that x(t)≥m2 and z(t)≥m2 for sufficiently large t. This can be done in the following two steps.
Step 1.
Choosem1>0,ε1>0 small enough that when (2.30) and (2.11) hold,
(2.32)σ≡exp⁡(aT-bM0T-c1Γε1-c3ε1T-c5M0T)<1,ρ≡exp⁡(-d2T-c4ε1T+c4(d1+c3m1)e2(θ1ε1-θ2ε1))<1,
where
(2.33)Γε1=p(exp⁡(Td1)-exp⁡(Tc2ε1))(1-q)(1-d1)exp⁡(Tc2ε1)-(d1+c3m1)exp⁡(Td1),θ1ε1=tan-1(e2(-pexp⁡(-(d1+c3m1)(n-1)T)+ε1exp⁡((d1+c3m1)T)-(1-q))Λm1),θ2ε1=tan-1(e2(-p-ε1(1-q)exp⁡(-(d1+c3m1)T))Λm1),Λm1=e2(1-(1-q)exp⁡(-(d1+c3m1)T)).
This step will show that x(t1)≥m1 and z(t1)≥m1 for some t1>0. Assuming the contrary, the following system results:
(2.34)u′(t)=(a-bM0-c1v(t)-c5M0)u(t),v′(t)=(-d1-c1m1)v(t),w′(t)=(-d2+c4v(t)1+e2(v(t))2)w(t)t≠nT,Δu(t)=0,Δv(t)=-qv(t)+p,Δw(t)=0.t=nT,
By Lemma 2.4, x(t)≥u(t), y(t)≥v(t), and z(t)≥w(t). By Lemma 2.6, v*(t)+ε1≥v(t)≥v*(t)-ε1, v*(t)=(pexp(-(d1+c1m1)(t-nT)))/(1-(1-q)exp(-(d1+c1m1)T)),(t∈(nT,(n+1)T]). Thus,
(2.35)u′(t)≥(a-bM0-c5M0-c1ε1-c1v*(t))u(t),w′(t)≥(-d2+c4(v*(t)-ε1)1+e2(v*(t)-ε1)2)w(t)≥(-d2-c4ε1+c4v*(t)1+e2(v*(t)-ε1)2).
Integrating (2.35) over (nT,(n+1)T] yields
(2.36)u((n+1)T)≥u(nT+)exp⁡∫nT(n+1)T(a-c1v*(t)-c1ε1-bM0-c5M0)dt≥u(nT)σ,w((n+1)T)≥w(nT+)exp⁡∫nT(n+1)T(-d2-c4ε1+c4v*(t)1+e2(v*(t)-ε1)2)dt≥w(nT)ρ.
Therefore,x((n+k)T)≥u((n+k)T)≥u(nT)σk and w((n+k)T)≥w((n+k)T)≥w(nT)σk.
Therefore,σk→∞ and ρk→∞ as k→∞. This implies that x(t)→∞ and z(t)→∞ as t→∞, which contradicts the boundedness of x(t) and z(t).
Step 2.
Ifx(t)≥m1 and z(t)≥m1 for all t≥t1, then the proof is complete. If not, let t*=inft≥t1{x(t)<m1,z(t)<m1}; then x(t)≥m1 and z(t)≥m1 for t∈[t1,t*) and x(t*)=m1,z(t*)=m1. By Step 1, there exists a t′>t* such that z(t′)≥m1. Set t2=inft>t*{x(t)≥m1,z(t)≥m1}; then x(t)<m1 and z(t)<m1 for t∈(t*,t2), and z(t2)=m1. This process can be continued by repeating Step 1. If this process stops after a finite number of repetitions, the proof is complete. If not, there exists an interval sequence [tn,tn+1],(n∈N), such that x(t)≥m1 and z(t)≥m1, t∈[tn,tn+1],(n∈N). Let T′=sup{tn+1-tn∣n∈N}. If T′=∞, there must exist a subsequence {tni} such that tni+1-tni→∞ as ni→∞. From Step 1, this can lead to a contradiction with the boundedness of x(t) and z(t). Therefore, T′<∞. Then,
(2.37)x(t)≥x(tn)exp⁡∫tnt(a-c1v*(s)-c1ε1-bM0-c5M0)ds≥m1exp⁡(-(bM0+c5M0)T′)≡m2,z(t)≥z(tn)exp⁡∫tnt(-d2-c4ε1+c4v*(s)1+e2(v*(s)-ε1)2)ds≥m1exp⁡(-d2T′)≡m3.
Let m4=min{m2,m3}; then liminft→∞x(t)≥m4, and liminft→∞z(t)≥m4. This completes the proof.


## 3. Numerical Analysis

### 3.1. Bifurcation

To study the dynamics of system (1.3), a period T was chosen, and the impulsive control parameter p was defined as the bifurcation parameter. The bifurcation diagram provides a summary of the essential dynamic behavior of the system [27–29].First, the influence of the periodT will be investigated. The bifurcation diagrams are shown in Figure 1. The results are in accordance with Theorem 2.8. Next, the influence of the impulsive control parameter p will be examined. A time series of the system response is shown in Figure 2. It is apparent that the solution (0,y*(t),0), which is said to be globally stable, is in agreement with Theorem 2.9. The bifurcation diagrams are shown in Figure 3.Bifurcation diagram of system (1.3) with initial conditions x(0)=0.1, y(0)=0.1, z(0)=0.1, and 0≤T≤50; a=1.85, b=0.65, c1=0.65, c2=0.55, c3=0.45, c4=0.15, c5=0.55, c6=0.8, e1=0.6, e2=0.6, d1=0.1, d2=0.3, P=0.1, and q=0.8.
(a)(b)(c)Dynamics of system (1.3) with a=1.72, b=0.52, c1=0.83, c2=0.64, c3=0.23, c4=0.32, c5=0.76, c6=0.48, e1=0.56, e2=0.62, d1=0.12, d2=0.25, T=13, p=4.2, and q=0.9. (a) time series of prey, (b) time series of intermediate predator, and (c) time series of top predator.
(a)(b)(c)Bifurcation diagram of system (1.3) with initial conditions x(0)=0.1, y(0)=0.1, z(0)=0.1, and 0≤p≤5; a=1.72, b=0.52, c1=0.83, c2=0.64, c3=0.23, c4=0.32, c5=0.76, c6=0.48, e1=0.56, e2=0.62, d1=0.12, d2=0.25, T=13, and q=0.9.
(a)(b)(c)To see the dynamics of system (1.3) clearly, a phase diagram with a different value of the parameter p corresponding to the bifurcation diagrams in Figure 3 is shown in Figure 4.Periodic behavior and chaos corresponding to Figure3: (a) phase diagram for p=0.05, (b) phase diagram for p=0.5, (c) phase diagram for p=2, (d) phase diagram for p=2.75, (e) phase diagram for p=2.85, (f) phase diagram for p=3.05, (g) phase diagram for p=3.15, (h) phase diagram for p=3.4, and (i) phase diagram for p=4.3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figures1 and 3 illustrate the complex dynamics of system (1.3), such as period-doubling cascades, symmetry-breaking pitchfork bifurcations, chaos, and nonunique dynamics. Because all the bifurcation diagrams are similar, only one of them will be explained. Part (a) of Figure 3 will be taken as an example. When p∈[0,0.414], the dynamics of the system are not obvious, but as p increases, the dynamics become more evident. The system enters into periodic windows with a chaotic band, as shown in Figures 4(a) (not obviously chaotic) and 4(b) (T-periodic solution). When p is between 0.414 and 2.689, chaos is intense, as shown in Figure 4(c). When p moves beyond 2.689, the chaotic behavior disappears, and periodic windows appear, as shown in Figures 4(d) (4T-periodic solution) and Figure 4(e) (2T-periodic solution). When p∈[3.019,3.076], the chaotic attractor increases in strength, and chaos reappears, as shown in Figure 4(f). For p greater than 3.076, periodic windows appear, as shown in Figure 4(g). For p in the interval 3.258 to 3.486, chaos reemerges, as shown in Figure 4(h). As the value of p increases further, the system enters into a stable state, as shown in Figure 4(i).


### 3.2. The Largest Lyapunov Exponent

The largest Lyapunov exponent is always calculated to detect whether a system is exhibiting chaotic behavior. The largest Lyapunov exponent takes into account the average exponential rates of divergence or convergence of nearby orbits in phase space [30]. A positive largest Lyapunov exponent indicates that the system is chaotic. If the largest Lyapunov exponent is negative, then periodic windows or a stable state must exist. Using the largest Lyapunov exponent, it is possible to see when the system is chaotic, at what time periodic windows disappear, and when the system is stable. The largest Lyapunov exponents corresponding to Figures 1 and 3 were calculated and are plotted in Figures 5 and 6.Figure 5
The largest Lyapunov exponents (LLE) corresponding to Figure1.Figure 6
The largest Lyapunov exponents (LLE) corresponding to Figure3.


### 3.3. Strange Attractors and Power Spectra

To understand the qualitative nature of strange attractors, power spectra can be used [31]. From the discussion in Section 3.2, it is known that the largest Lyapunov exponent for strange attractor (a) is 0.475 and that for strange attractor (b) is 0.0144. This means that these are both chaotic attractors, and the fact that the exponent of (a) is larger than that of (b) means that the chaotic dynamics of (a) are more intense than those of (b). Considering the power spectra of these attractors, the spectrum of strange attractor (b) is composed of strong broadband components and sharp peaks, as shown in Figure 7(d). By contrast, in the spectrum of the strong chaotic attractor (a), it is not easy to distinguish any sharp peaks, as can be seen in Figure 7(c). By means of power spectra, it is possible to determine that (b) originates in a strong limit cycle, but that (a) comes from a set of weak limit cycles.Strange attractors and power spectra: (a) strange attractor whenp=2.2, (b) strange attractor when p=3.039, (c) power spectrum of attractor (a), (d) power spectrum of attractor (b).
(a)(b)(c)(d)


## 3.1. Bifurcation

To study the dynamics of system (1.3), a period T was chosen, and the impulsive control parameter p was defined as the bifurcation parameter. The bifurcation diagram provides a summary of the essential dynamic behavior of the system [27–29].First, the influence of the periodT will be investigated. The bifurcation diagrams are shown in Figure 1. The results are in accordance with Theorem 2.8. Next, the influence of the impulsive control parameter p will be examined. A time series of the system response is shown in Figure 2. It is apparent that the solution (0,y*(t),0), which is said to be globally stable, is in agreement with Theorem 2.9. The bifurcation diagrams are shown in Figure 3.Bifurcation diagram of system (1.3) with initial conditions x(0)=0.1, y(0)=0.1, z(0)=0.1, and 0≤T≤50; a=1.85, b=0.65, c1=0.65, c2=0.55, c3=0.45, c4=0.15, c5=0.55, c6=0.8, e1=0.6, e2=0.6, d1=0.1, d2=0.3, P=0.1, and q=0.8.
(a)(b)(c)Dynamics of system (1.3) with a=1.72, b=0.52, c1=0.83, c2=0.64, c3=0.23, c4=0.32, c5=0.76, c6=0.48, e1=0.56, e2=0.62, d1=0.12, d2=0.25, T=13, p=4.2, and q=0.9. (a) time series of prey, (b) time series of intermediate predator, and (c) time series of top predator.
(a)(b)(c)Bifurcation diagram of system (1.3) with initial conditions x(0)=0.1, y(0)=0.1, z(0)=0.1, and 0≤p≤5; a=1.72, b=0.52, c1=0.83, c2=0.64, c3=0.23, c4=0.32, c5=0.76, c6=0.48, e1=0.56, e2=0.62, d1=0.12, d2=0.25, T=13, and q=0.9.
(a)(b)(c)To see the dynamics of system (1.3) clearly, a phase diagram with a different value of the parameter p corresponding to the bifurcation diagrams in Figure 3 is shown in Figure 4.Periodic behavior and chaos corresponding to Figure3: (a) phase diagram for p=0.05, (b) phase diagram for p=0.5, (c) phase diagram for p=2, (d) phase diagram for p=2.75, (e) phase diagram for p=2.85, (f) phase diagram for p=3.05, (g) phase diagram for p=3.15, (h) phase diagram for p=3.4, and (i) phase diagram for p=4.3.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figures1 and 3 illustrate the complex dynamics of system (1.3), such as period-doubling cascades, symmetry-breaking pitchfork bifurcations, chaos, and nonunique dynamics. Because all the bifurcation diagrams are similar, only one of them will be explained. Part (a) of Figure 3 will be taken as an example. When p∈[0,0.414], the dynamics of the system are not obvious, but as p increases, the dynamics become more evident. The system enters into periodic windows with a chaotic band, as shown in Figures 4(a) (not obviously chaotic) and 4(b) (T-periodic solution). When p is between 0.414 and 2.689, chaos is intense, as shown in Figure 4(c). When p moves beyond 2.689, the chaotic behavior disappears, and periodic windows appear, as shown in Figures 4(d) (4T-periodic solution) and Figure 4(e) (2T-periodic solution). When p∈[3.019,3.076], the chaotic attractor increases in strength, and chaos reappears, as shown in Figure 4(f). For p greater than 3.076, periodic windows appear, as shown in Figure 4(g). For p in the interval 3.258 to 3.486, chaos reemerges, as shown in Figure 4(h). As the value of p increases further, the system enters into a stable state, as shown in Figure 4(i).


## 3.2. The Largest Lyapunov Exponent

The largest Lyapunov exponent is always calculated to detect whether a system is exhibiting chaotic behavior. The largest Lyapunov exponent takes into account the average exponential rates of divergence or convergence of nearby orbits in phase space [30]. A positive largest Lyapunov exponent indicates that the system is chaotic. If the largest Lyapunov exponent is negative, then periodic windows or a stable state must exist. Using the largest Lyapunov exponent, it is possible to see when the system is chaotic, at what time periodic windows disappear, and when the system is stable. The largest Lyapunov exponents corresponding to Figures 1 and 3 were calculated and are plotted in Figures 5 and 6.Figure 5
The largest Lyapunov exponents (LLE) corresponding to Figure1.Figure 6
The largest Lyapunov exponents (LLE) corresponding to Figure3.


## 3.3. Strange Attractors and Power Spectra

To understand the qualitative nature of strange attractors, power spectra can be used [31]. From the discussion in Section 3.2, it is known that the largest Lyapunov exponent for strange attractor (a) is 0.475 and that for strange attractor (b) is 0.0144. This means that these are both chaotic attractors, and the fact that the exponent of (a) is larger than that of (b) means that the chaotic dynamics of (a) are more intense than those of (b). Considering the power spectra of these attractors, the spectrum of strange attractor (b) is composed of strong broadband components and sharp peaks, as shown in Figure 7(d). By contrast, in the spectrum of the strong chaotic attractor (a), it is not easy to distinguish any sharp peaks, as can be seen in Figure 7(c). By means of power spectra, it is possible to determine that (b) originates in a strong limit cycle, but that (a) comes from a set of weak limit cycles.Strange attractors and power spectra: (a) strange attractor whenp=2.2, (b) strange attractor when p=3.039, (c) power spectrum of attractor (a), (d) power spectrum of attractor (b).
(a)(b)(c)(d)


## 4. Conclusions and Remarks

In this paper, the dynamic behavior of a predator-prey model with Holling type IV functional response with respect to an impulsive control strategy has been investigated. The conditions for locally asymptotically stable and globally stable periodic solutions and for system permanence have been determined. It has been determined that an impulsive control strategy changes the dynamic behavior of the model. Complex dynamic patterns also have been observed in continuous-time predator-prey or three-species food-chain models [32–34]. Numerical simulations were performed to obtain bifurcation diagrams with respect to the period T and the impulsive control parameter p. Using computer-based simulation, phase diagrams were generated to reveal details of the bifurcation behavior, such as period-doubling cascades, symmetry-breaking pitchfork bifurcations, chaos, and nonunique dynamics. The largest Lyapunov exponents were also simulated to verify the chaotic dynamics of the system. In addition, power spectra were used to understand the qualitative nature of strange attractors. On the basic of numerical simulations, the main difference of this paper from other papers is that the top predator z cannot be persistent although the mid-level predator y is persistent. This is because the mid-level predator y is easy to escape the top predator z, which is a variable to the competitive. It can be concluded that an impulsive control strategy is an effective method to control the system dynamics of a predator-prey model.



---

*Source: 141272-2012-02-07.xml*