# A Holling Type II Pest and Natural Enemy Model with Density Dependent IPM Strategy

Authors: Xia Wang; Yuan Tian; Sanyi Tang
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

Resource limitations and density dependent releasing of natural enemies during the pest control and integrated pest management will undoubtedly result in nonlinear impulsive control. In order to investigate the effects of those nonlinear control strategies on the successful pest control, we have proposed a pest-natural enemy system concerning integrated pest management with density dependent instant killing rate and releasing rate. In particular, the releasing rate depicts how the number of natural enemy populations released was guided by their current density at the fixed moment. The threshold condition which ensures the existence and global stability of pest-free periodic solution has been discussed first, and the effects of key parameters on the threshold condition reveal that reducing the pulse period does not always benefit pest control; that is, frequent releasing of natural enemies may not be beneficial to the eradication of pests when the density dependent releasing method has been implemented. Moreover, the forward and backward bifurcations could occur once the pest-free periodic solution becomes unstable, and the system could exist with very complex dynamics. All those results confirm that the control actions should be carefully designed once the nonlinear impulsive control measures have been taken for pest management.


---

## Body

## 1. Introduction

Integrated pest management (IPM), also known as integrated pest control (IPC), is a wide approach that integrates multiple control measures for economic control of pests. IPM aims to maintain the pest population below the economic injury level (EIL) by taking biological, chemical, and physical control tactics. During the pest control, regular monitor is critically important which provides some cues on how and when to apply IPM strategies [1–5].Note that the outbreaks of pest often cause serious ecological problems and bring large economic losses. The traditional methods to reduce the harm to agriculture due to insect pests include chemical pesticides and biological control measures [6–8]. Chemical pesticide is a measure directly acting on the pest. Biological control measure releasing natural enemies is a measure to increase the number of predators to control the population of the pest. There are several different releasing methods including inoculative release and inundative release which are different tactics of biological control that affect the target pest in different ways. Augmentative control includes the periodic introduction of natural enemies [3, 9]. With inundative release, predators are collected, mass-reared, and periodically released in large numbers into the pest area, which can result in an immediate reduction in the pest populations.However, pesticide pollution has been recognized as a major health hazard to human, it can seriously change the ecological system, and it may cause damage to nontarget organisms. So, to combine these control measures together and to maximize control effects at minimum levels of control resources, a number of impulsive differential equations have been developed to investigate the dynamics of pest and enemies and to analyze the optimal control schemes [3, 4, 10–13].In 2003, Liu and Chen studied the dynamics of the Holling type 2 predator-prey model by introducing releasing predators at fixed moments [14], and the model is as follows:(1)dxtdt=rxt1-xtK-axtyt1+ωxt,dytdt=cxtyt1+ωxt-dyt,t≠nT,xt+=xt,yt+=yt+σ,t=nT,where x and y are densities of the prey and predator populations at time t, respectively; r denotes the intrinsic growth rate of prey; K represents the environmental capacity of prey; d is the death rate of the predator; ax(t)/1+ωx(t) is the Holling type 2 functional response; c/a is the rate of conversing prey into predator; T is the pulse period; σ is the release amount of predator at time points t=nT. From then on, the predator-prey model considering both pulse spraying pesticides and releasing natural enemies have been extensively studied by many researchers [3, 4, 13, 15–21]. Note that for all those studies the pulse functions introduced to the models could be summarized as the following generalized linear functions:(2)xt+=1-p1xt,yt+=1-p2yt+σ,t=nT,where p1 and p2 are the fraction of pests and predators which have been killed instantaneously after pesticide applications.The key assumptions from the linear pulse functions shown in (2) are as follows: (a) the instant killing rate p1 is a constant, which shows that no matter how large the number of pest populations is, the total pest killed linearly depends on its current density [22, 23]; (b) the adverse effect of pesticides on the natural enemies is also linear dependence on its current density; (c) the releasing number of natural enemies is a constant no matter how large the natural enemies remain in the field; that is, the density of natural enemies is not carefully monitored before IPM strategy is applied. However, the killing rate and releasing constant should strictly depend on a number of factors including current densities of both populations, resources limitation, and management purposes. For example, in reality, due to the resource limitation, the more proper and usual way of releasing of enemies may be the way that is guided by the density of the current natural enemies. Therefore, taking those factors into account, we propose the following model with density dependent pulse control actions:(3)dxtdt=rxt1-xtK-axtyt1+ωxt,dytdt=cxtyt1+ωxt-dyt,t≠nT,xt+=1-δxtxt+hxt,yt+=yt+τ1+θyt,t=nT,where δ>0 and h>0 represent the maximal fatality rate and the half-saturation constant for the prey, τ>0 is the maximal release amount of the predator, and θ is a shape parameter. We assume that the control measures are conducted every T time, and thus when t=nT the amounts of the prey and the predator are updated to 1-δx(t)/x(t)+hx(t) and y(t)+τ/1+θy(t) instantaneously. Note that the function τ/1+θy(t) representing the amount of predators releasing at t=nT is a decreasing function with respect to y(t). Note that the amount of natural enemies releasing at t=nT depends on its density y(nT), which indicates that the larger the density of natural enemy is, the smaller the natural enemy has been released at time t=nT, and obviously this strategy is more effective than the constant releasing method [3, 4, 16, 24].The main purpose of the present paper is to focus on how the density dependent killing rate and density dependent releasing strategy affect the dynamics of the system and the effectiveness of the control measures. Based on this purpose, we first investigate the threshold condition which guarantees the existence and global stability of pest-free periodical solution. The sensitivity analyses for the key parameters on the threshold condition depict that reducing the pulse period may not be of benefit to the eradication of pests when the density dependent releasing method has been implemented, and this interesting result revealed by the proposed model shows that incorrect control methods may result in pest outbreak and resurgence. Furthermore, bifurcation analyses indicate that the forward and backward bifurcations could occur once the pest-free periodic solution becomes unstable, and the system could exist with two stable positive periodic solutions and even more complex dynamics. The important results found in this work confirm that the control actions should be carefully designed once the nonlinear impulsive control measures have been taken for pest management.


## 2. Mathematical Analysis and Threshold Conditions

In this section, we first demonstrate the existence of the pest-free periodic solution of model (3) and focus on its stability.


### 2.1. Existence and Stability of the Pest-Free Periodic Solution

In order to investigate the existence of the pest-free periodic solution, we first focus on the pest-free set{(x,y)∈R+2,x=0}, which is invariantly related to model (3). Within this set, model (3) becomes(4)dytdt=-dyt,t≠nT,yt+=yt+τ1+θyt,t=nT,y0+=y0.Solving the system in interval (nT,(n+1)T], one yields (5)yn+1T+=ynT+e-dT+τ1+θynT+e-dT.Denote y(nT+)=yn; then the equation can be rewritten as follows:(6)yn+1=yne-dT+τ1+θyne-dT≐fyn.This equation is the so-called stroboscopic map of model (4). It describes the relations of the number of natural enemies between any two successive pulse points. So, the existence of the positive steady state of this equation implies the existence of a positive periodic solution of model (4). Thus, the steady state of (6) is discussed first. The derivative of f(yn) with respect to yn can be given as follows: (7)f′yn=e-dT1-τθ1+θyne-dT2.Denotey∗ to be the fixed point of the stroboscopic map (6); then y∗ satisfies the equation f(y∗)=y∗e-dT+τ/1+θy∗e-dT=y∗. So, y∗ satisfies the following equation:(8)A1y∗2+A2y∗+A3=0,where A1=θe-dT(1-e-dT), A2=1-e-dT, and A3=-τ. Denote Δ=A22-4A1A3; then we have (9)Δ=1-e-dT2+4τθe-dT1-e-dT>0.So, (6) has a unique positive root:(10)y∗=-1-e-dT+Δ2θe-dT,which is stable provided (11)f′ynyn=y∗=e-dT1-τθ1+θy∗e-dT2<1.Obviously,e-dT(1-τθ/(1+θy∗e-dT)2)<1. Moreover, (1+θy∗e-dT)2=τθe-dT+1/2(e-dT-1)2+1/2(e-dT-1)Δ+e-dT+Δ, which indicates that (12)τθ1+θy∗e-dT2=τθτθe-dT+1/2e-dT-12+1/2e-dT-1Δ+e-dT+Δ<edT.Therefore, (13)e-dT1-τθ1+θy∗e-dT2>e-dT-1>-1and consequently y∗ is stable. Then, we can conclude that model (4) has a unique nontrivial positive periodic solution, denoted by yp(t), which can be expressed as follows:(14)ypt=y∗e-dt-nT,fornT<t≤n+1T.Further, yp(t) is locally stable because y∗ is stable [25]. Then, we can obtain the following theorem by proving the global attractivity of y∗.Theorem 1.
Model (4) has a positive periodic solution yp(t) with period T and for every solution y(t) of (4) one has |y(t)-yp(t)|→0 as t→+∞, where yp(t)=y∗e-d(t-nT) (nT<t≤(n+1)T) and y∗=-(1-e-dT)+Δ/2θe-dT.Proof.
Letf(y)=ye-dT+τ/1+θye-dT. To study the properties of the function f(y), we first solve the equation df(y)/dy=e-dT1-τθ/(1+θye-dT)2=0, with respect to y, and get two stationary points (15)y1=-τθ-1θedT<0,y2=τθ-1θedT.The second derivative of the function f(y) is d2f/dy2=2τ(θe-dT)2/(1+θye-dT)3. It is obvious that d2f/dy2(y1)<0 and d2f/dy2(y2)>0. So, y1 is the local maximum and y2 is the local minimum of the function f(y) with y2>y1 and y1<0.
Thus, ify2>0, the function f must be a convex function and if y2<0f must be monotonous in the first quadrant. So, based on the sign of y2 and the positional relations between point y2 and y∗, we consider the following three possible cases.
Case A
(
y
2
<
0
<
y
∗
). For this case the function f(y) is a monotonically increasing function for y>y2. For any y∈[0,y∗), according to y<f(y)<y∗ we know that fk(y) is monotonically increasing as k increases, and fk(y)→y∗(k→∞).
For anyy>y∗, we have fk(y)>y∗ for all k, and then according to f(y)<y we conclude that fk(y) is monotonically decreasing as k increases, which means that fk(y)→y∗(k→∞).
Case B
(
0
<
y
2
<
y
∗
). For this case, if y∈[y2,y∗) or y>y∗, then the conclusions corresponding to the Case A are true; that is, we have fk(y)→y∗(k→∞).
Moreover, for anyy∈[0,y2), there must be f(y)>y2, and according to the conclusion of Case A we obtain f1+k(y)→y∗(k→∞).
Case C
(
y
2
>
y
∗
>
0
). For this case, it is easy to know that the function f(y) is a monotonically decreasing function for y∈[0,y2], and it is an increasing function for y∈(y2,∞).
It is easy to see that there must exist a positive integerl such that fl(y)∈[y∗,y2] for y∈[0,y∗) or y∈(y2,∞). Thus if fk(y)→y∗(k→∞) for any y∈(y∗,y2], then fn(y)→y∗ (n→∞) for any y∈[0,∞) and the results follow. Hence in the following we will focus on fk(y)→y∗(k→∞) for any y∈(y∗,y2].
According to the definition of the functionf and by simple calculations, we have (16)fy=ye-dT+τ1+θye-dT,ffy=fye-dT+τ1+θfye-dT,f′y=e-dT1-τθ1+θye-dT2,f′fy=e-dT1-τθ1+θfye-dT2.
Now we first show thatf′(f(y))>-f(f(y))/f(y) for any y∈(y∗,y2]. In fact, for y∈(y∗,y2], f′(f(y))>-f(f(y))/f(y)⇔f′(f(y))f(y)>-f(f(y)), since (17)f′fyfy=e-dT1-τθ1+θfye-dT2fy,-ffy=-fye-dT-τ1+θfye-dT,2fye-dT>0>-τ1+θfye-dT2,which indicates that (18)f′fyfy>-ffy.
Next, we show thatf′(y)>-f(y)/y for any y∈(y∗,y2]. In fact, for y∈(y∗,y2] we have f′(y)>-f(y)/y⇔f′(y)y>-f(y). Since (19)f′yy=e-dT1-τθ1+θye-dT2y,-fy=-ye-dT-τ1+θye-dT,2q2e-dTy>0>-τ1+θye-dT2,thus we have (20)f′fyy>-fy.
It follows fromf′(f(y))>-f(f(y))/f(y) that we know (21)f′fyf′yy<-f′yyfyffyand f′(y)y>-f(y) for f′(y)>-f(y)/y; that is, -f′(y)y/f(y)<1 holds. Furthermore, we have 0<-f′(y)y/f(y)<1 due to -f′(y)y/f(y)>0. All those confirm that f′fyf′(y)y<-f′(y)y/f(y)f(f(y))<f(f(y)); namely,(22)f′fyf′yy<ffy.
Denoteg(y)=f(f(y))/y; then we have g′(y)=f′fyf′(y)y-f(f(y))/y2. Then, g′(y)<0 and g(y) is a monotonically decreasing function at interval y∈(y∗,y2] due to (22).
Therefore,g(y∗)=f(f(y∗))/y∗=f(y∗)/y∗=y∗/y∗=1, and g(y)<1 for y∈(y∗,y2]; that is, f(f(y))<y.
Based on the above discussions we conclude thatf2k(y) is a monotonically decreasing function as k increases with f2k(y)→y∗(k→∞). Similarly, we can prove that f2k+1(y) is a monotonically decreasing function as k increases with f2k(y)→y∗(k→∞). Thus, fn(y)→y∗(n→∞) for any y∈[0,∞) holds.
Combined with all cases together we know that the unique equilibriumy∗ is globally stable if it exists, and consequently the periodic solution yp(t) of model (4) is globally stable. This completes the proof.The stability of the positive periodic solution is determined by the stability of the positive equilibrium of the difference equation made by sequence of impulsive points; thus the positive periodic solutionyp(t) of model (4) is globally stable according to Theorem 1. Therefore, a general expression of unique pest-free periodic solution of model (3) over the interval nT<t≤(n+1)T for all n∈N can be denoted by(23)0,ypt=0,y∗e-dt-nT.In the following, we focus on the stability of the pest-free periodic solution for model (3). Firstly, define(24)x¯t,y¯t=xt,yt-0,ypt,and then (x¯(t),y¯(t)) represents the deviations form of the pest-free solution. According to model (3), the derivations of (x¯(t),y¯(t)) are (25)x¯˙=rx¯-rKx¯2-ax¯1+ωx¯y¯+ypt,y¯˙=cx¯1+ωx¯y¯+ypt-dy¯.Assuming that (x¯(t),y¯(t)) are small enough, then the linear approximation of the deviation system around the periodic solution (0,yp(t)) can be obtained as follows:(26)x¯˙=r-ayptx¯,y¯˙=cyptx¯-dy¯.Then, we have the following main result.Theorem 2.
The pest-free periodic solution(0,yp(t)) of model (3) is locally stable in the first quadrant provided that (27)R1≐rdTay∗1-e-dT<1.Proof.
LetΦ(t) be the fundamental matrix of (26); then (28)ΦT=e∫0Tr-ayptdt0∗e-dT,and Φ(0)=I is the identity matrix. The term ∗ is not required in the following analysis.
The resetting impulsive conditions of (3) become (29)x¯nT+y¯nT+=1001-τθ1+θypnT2x¯nTy¯nT=BnTx¯nTy¯nT.Thus, the Floquet multipliers of the matrix (30)M=BTΦTare(31)λ1=e∫0Tr-ayptdt,λ2=1-τθ1+θypT2e-dT.It follows from the literature [26] that the pest-free solution is locally stable if and only if |λ1|<1 and |λ2|<1. The former inequality is equal to R1≐rdT/ay∗(1-e-dT)<1 and the latter inequality is true if and only if R2≐1-τθ/(1+θy∗e-dT)2e-dT<1, which is true naturally based on the analysis in the previous section. Therefore, the pest-free periodic solution is locally stable when R1<1.


### 2.2. Boundness and Permanence of System (3)

According to the proof of Theorem2, the pest-free periodic solution is unstable when R1>1. So, next we aim to prove the permanence of the system if R1>1. First we prove the boundness of system (3).Theorem 3.
There exists a constantM0>0 such that x(t)≤M0 and y(t)≤M0 for each solution of system (3) with all t large enough.Proof.
Supposex(t) and y(t) are any solution of system (3). Let V(t)=(c/a)x(t)+y(t). Then, when t≠nT, (32)dVtdt+λVt=car+λx-craKx2-λ-dy.When 0<λ≤d, the right hand of the equation is bounded; namely, (33)dVtdt+λVt≤M1,with M1>0. When t=nT, (34)Vt+=Vt+τ1+θyt-caδxth+xt≤Vt+τ.
Considering the system(35)dVtdt≤-λVt+M1,t≠nT,Vt+≤Vt+τ,t=nT,we have (36)Vt≤V0+-M1λe-λt+τ1-e-nλT1-e-λTe-λt-nT+M1λ.So, V(t) is uniformly ultimately bounded. Hence, there exists a constant M0>0 such that x(t)≤M0 and y(t)≤M0 for each solution of system (3) with all t large enough. This completes the proof.Theorem 4.
System (3) is permanent provided R1>1.Proof.
It follows from Theorem1 that there exists an m1>0 such that y≥m1 for t large enough. Thus, in order to show the permanence of model (3) we shall next find a constant m2>0, such that x≥m2 for t large enough.
Considering the second equation of system (3), we have(37)dytdt≤-dyt,t≠nT,yt+=yt+τ1+θyt,t=nT.Therefore, according to the comparison theorem and the analyses in the previous subsection, we have (38)yt≤ypt+ε1for t>t1 and ε1>0. Then, by the first equation of system (3) we have(39)dxtdt≥rxt1-xtK-aypt+ε1xt,t≠nT,xt+=1-δxtxt+hxt,t=nT.Considering the auxiliary system(40)dutdt=rut1-utK-aypt+ε1ut,t≠nT,ut+=1-δutut+hut,t=nTand using the same method as system (4), we have (41)un+1+=un+11-δun+1un+1+hwith(42)un+1=un+e∫0Tr-aypt+ε1dt1+un+∫0Tr/Ke∫0tr-aypξ+ε1dξdt.Let M1=e∫0Tr-aypt+ε1dt and M2=∫0Tr/Ke∫0t(r-a(yp(ξ)+ε1))dξdt; then un+1=un+M1/1+M2un+. Therefore, we have the following different equation with respect to un+:(43)un+1+=un+M11+M2un+1-δun+M1/1+M2un+h+un+M1/1+M2un+≐Fun+.Let u∗ be the steady state of the above difference equation, and by simple calculation we have (44)a1u∗2+a2u∗+a3=0,where a1=M22+M2M1, a2=M2+M1+M2h-M2M1-M12+δM12, and a3=h(1-M1). Obviously, a1>0. Moreover, a3<0 when R1>1 for ε1 sufficiently small, which shows that the above difference equation exists with a unique positive equilibrium u∗. Therefore, system (40) has only one positive periodic solution when R1>1, denoted by u¯(t).
To study the stability of this periodic solution, we first investigate the properties of the functionF(u)=M1u/1+M2u(1-δM1u/h(1+M2u)+M1u). Taking the derivation of F(u) with respect to u, one yields (45)F′u=M1hM2+M1u+h2-2δM12hM2+δM23u2-2δM12hu1+M2u2hM2+M1u+h2.Let F′(u)=0; then we get the following equation: (46)b1u2+b2u+b3=0,where b1=M1M22h2+2(1-δ)hM12M2+(1-δ)M13,b2=2h2M1M2+2(1-δ)hM12, and b3=h2M1. It follows from 0<δ<1 that b1>0,b2>0, and b3>0. Moreover, we can prove that b22-4b1b3=4M14h2δ(1-δ)>0. Therefore, the equation F′(u)=0 does not exist with any positive root. Thus, F′(u)>0 for all u>0; that is, the function F(u) is a monotonically increasing function for all u>0. Further, it is easy to show that F′(0)=M1 according to the expression of F′(u) with M1>1 for R1>1 and ε1 sufficiently small. All those confirm that the intersection angle of the function F(u) and the horizontal coordinate is larger than π/4. Therefore, the function image is as shown in Figure 1. Then, according to u<F(u)<u∗ we have the fact that Fk(u) is monotonically increasing as k increases for any u∈[0,u∗), and Fk(u)→u∗(k→∞); according to F(u)<u we have the fact that Fk(u) is monotonically decreasing as k increases for any u>u∗, and Fk(u)→u∗(k→∞). Therefore, the equilibrium u∗ is globally stable.
Thus, the positive periodic solutionu¯(t) is globally stable. Then, according to the comparison theory, x(t)≥u(t)≥u¯(t)-ε2 for t≥t2. Let m2=min(u¯(t))-ε2, ε2>0; then we have x≥m2 for t large enough. Thus, system (3) is permanent provided R1>1. This completes the proof.Figure 1
Illustration of the global stability ofu∗.


### 2.3. Bifurcation of Positive Periodic Solution

In this section, we investigate the existence and stability of positive periodic solution of system (3) by employing the bifurcation analysis related to the pest-free periodic solution X0=(0,yp(t)) [26].Theorem 5.
As the parameterT passes through the critical value T0 satisfying e∫0T0[r-ayp(ξ)]dξ=1, a nontrivial positive periodic solution appears near X0. The bifurcation is supercritical provided BC<0, while it is subcritical if BC>0, where (47)B=-r-aypT0+adT0y∗T0C2e-dT01-C2e-dT0e∫0T0r-aypξdξ,(48)C=2δh-1e2∫0T0r-aypξdξ-2aT01-C1C21-C2e-dT0e∫0T0r-aypξdξ-C3.Proof.
To use the theory of [26], we employ the same notations as those in [26]; then (49)F1x,y=cxy1+ωx-dy,F2x,y=rx1-xK-axy1+ωx,Θ1x,y=y+τ1+θy,Θ2x,y=1-δxx+hx.By simple calculations, we get (50)∂Φ1T0,X0∂y=e-dT0>0,∂Φ2T0,X0∂x=e∫0T0r-aypξdξ>0,∂Φ1T0,X0∂x=c∫0T0ypue-dT0-ue∫0ur-aypξdξdu≐C1>0,a0′=1-C2e-dT0,d0′=1-e∫0T0r-aypξdξ,b0′=1-C1C2,where C2=1-τθ/(1+θyp(T0))2 and T0 is the root of d0′=0, satisfying the condition cy∗(e-d(T0-nT)-1)/adT0=1. Obviously, we have the fact that C2<1, which shows that a0′>0. Further, we have (51)∂Φ1T0,X0∂T¯=-dy∗T0e-dT0,∂2Φ2T0,X0∂x∂y=-aT0e∫0T0r-aypξdξ<0,∂2Φ2T0,X0∂x∂T¯=r-aypT0e∫0T0r-aypξdξ,∂2Φ2T0,X0∂x2=T0-2rK+2aωypT0e∫0T0r-aypξdξ-∫0T0ae∫uT0r-aypξdξ∫0ue-du-vcypve∫0T0r-aypξdξdvdu≐C3.∂Θ1(T0,X0)/∂y=C2, ∂Θ1(T0,X0)/∂x=0, ∂Θ2/∂x=1, ∂Θ2(T0,X0)/∂y=0, ∂2Θ1(T0,X0)/∂x∂y=-2δh-1, and ∂2Θ1(T0,X0)/∂x∂y=0. Thus, (52)B=-r-aypT0+adT0y∗T0C2e-dT01-C2e-dT0e∫0T0r-aypξdξ,C=2δh-1e2∫0T0r-aypξdξ-2aT01-C1C21-C2e-dT0e∫0T0r-aypξdξ-C3.Therefore, according to Theorem 2 in literature [26], the supercritical branch occurs at the point T0 if BC<0 and the subcritical branch occurs if BC>0. This completes the proof.


## 2.1. Existence and Stability of the Pest-Free Periodic Solution

In order to investigate the existence of the pest-free periodic solution, we first focus on the pest-free set{(x,y)∈R+2,x=0}, which is invariantly related to model (3). Within this set, model (3) becomes(4)dytdt=-dyt,t≠nT,yt+=yt+τ1+θyt,t=nT,y0+=y0.Solving the system in interval (nT,(n+1)T], one yields (5)yn+1T+=ynT+e-dT+τ1+θynT+e-dT.Denote y(nT+)=yn; then the equation can be rewritten as follows:(6)yn+1=yne-dT+τ1+θyne-dT≐fyn.This equation is the so-called stroboscopic map of model (4). It describes the relations of the number of natural enemies between any two successive pulse points. So, the existence of the positive steady state of this equation implies the existence of a positive periodic solution of model (4). Thus, the steady state of (6) is discussed first. The derivative of f(yn) with respect to yn can be given as follows: (7)f′yn=e-dT1-τθ1+θyne-dT2.Denotey∗ to be the fixed point of the stroboscopic map (6); then y∗ satisfies the equation f(y∗)=y∗e-dT+τ/1+θy∗e-dT=y∗. So, y∗ satisfies the following equation:(8)A1y∗2+A2y∗+A3=0,where A1=θe-dT(1-e-dT), A2=1-e-dT, and A3=-τ. Denote Δ=A22-4A1A3; then we have (9)Δ=1-e-dT2+4τθe-dT1-e-dT>0.So, (6) has a unique positive root:(10)y∗=-1-e-dT+Δ2θe-dT,which is stable provided (11)f′ynyn=y∗=e-dT1-τθ1+θy∗e-dT2<1.Obviously,e-dT(1-τθ/(1+θy∗e-dT)2)<1. Moreover, (1+θy∗e-dT)2=τθe-dT+1/2(e-dT-1)2+1/2(e-dT-1)Δ+e-dT+Δ, which indicates that (12)τθ1+θy∗e-dT2=τθτθe-dT+1/2e-dT-12+1/2e-dT-1Δ+e-dT+Δ<edT.Therefore, (13)e-dT1-τθ1+θy∗e-dT2>e-dT-1>-1and consequently y∗ is stable. Then, we can conclude that model (4) has a unique nontrivial positive periodic solution, denoted by yp(t), which can be expressed as follows:(14)ypt=y∗e-dt-nT,fornT<t≤n+1T.Further, yp(t) is locally stable because y∗ is stable [25]. Then, we can obtain the following theorem by proving the global attractivity of y∗.Theorem 1.
Model (4) has a positive periodic solution yp(t) with period T and for every solution y(t) of (4) one has |y(t)-yp(t)|→0 as t→+∞, where yp(t)=y∗e-d(t-nT) (nT<t≤(n+1)T) and y∗=-(1-e-dT)+Δ/2θe-dT.Proof.
Letf(y)=ye-dT+τ/1+θye-dT. To study the properties of the function f(y), we first solve the equation df(y)/dy=e-dT1-τθ/(1+θye-dT)2=0, with respect to y, and get two stationary points (15)y1=-τθ-1θedT<0,y2=τθ-1θedT.The second derivative of the function f(y) is d2f/dy2=2τ(θe-dT)2/(1+θye-dT)3. It is obvious that d2f/dy2(y1)<0 and d2f/dy2(y2)>0. So, y1 is the local maximum and y2 is the local minimum of the function f(y) with y2>y1 and y1<0.
Thus, ify2>0, the function f must be a convex function and if y2<0f must be monotonous in the first quadrant. So, based on the sign of y2 and the positional relations between point y2 and y∗, we consider the following three possible cases.
Case A
(
y
2
<
0
<
y
∗
). For this case the function f(y) is a monotonically increasing function for y>y2. For any y∈[0,y∗), according to y<f(y)<y∗ we know that fk(y) is monotonically increasing as k increases, and fk(y)→y∗(k→∞).
For anyy>y∗, we have fk(y)>y∗ for all k, and then according to f(y)<y we conclude that fk(y) is monotonically decreasing as k increases, which means that fk(y)→y∗(k→∞).
Case B
(
0
<
y
2
<
y
∗
). For this case, if y∈[y2,y∗) or y>y∗, then the conclusions corresponding to the Case A are true; that is, we have fk(y)→y∗(k→∞).
Moreover, for anyy∈[0,y2), there must be f(y)>y2, and according to the conclusion of Case A we obtain f1+k(y)→y∗(k→∞).
Case C
(
y
2
>
y
∗
>
0
). For this case, it is easy to know that the function f(y) is a monotonically decreasing function for y∈[0,y2], and it is an increasing function for y∈(y2,∞).
It is easy to see that there must exist a positive integerl such that fl(y)∈[y∗,y2] for y∈[0,y∗) or y∈(y2,∞). Thus if fk(y)→y∗(k→∞) for any y∈(y∗,y2], then fn(y)→y∗ (n→∞) for any y∈[0,∞) and the results follow. Hence in the following we will focus on fk(y)→y∗(k→∞) for any y∈(y∗,y2].
According to the definition of the functionf and by simple calculations, we have (16)fy=ye-dT+τ1+θye-dT,ffy=fye-dT+τ1+θfye-dT,f′y=e-dT1-τθ1+θye-dT2,f′fy=e-dT1-τθ1+θfye-dT2.
Now we first show thatf′(f(y))>-f(f(y))/f(y) for any y∈(y∗,y2]. In fact, for y∈(y∗,y2], f′(f(y))>-f(f(y))/f(y)⇔f′(f(y))f(y)>-f(f(y)), since (17)f′fyfy=e-dT1-τθ1+θfye-dT2fy,-ffy=-fye-dT-τ1+θfye-dT,2fye-dT>0>-τ1+θfye-dT2,which indicates that (18)f′fyfy>-ffy.
Next, we show thatf′(y)>-f(y)/y for any y∈(y∗,y2]. In fact, for y∈(y∗,y2] we have f′(y)>-f(y)/y⇔f′(y)y>-f(y). Since (19)f′yy=e-dT1-τθ1+θye-dT2y,-fy=-ye-dT-τ1+θye-dT,2q2e-dTy>0>-τ1+θye-dT2,thus we have (20)f′fyy>-fy.
It follows fromf′(f(y))>-f(f(y))/f(y) that we know (21)f′fyf′yy<-f′yyfyffyand f′(y)y>-f(y) for f′(y)>-f(y)/y; that is, -f′(y)y/f(y)<1 holds. Furthermore, we have 0<-f′(y)y/f(y)<1 due to -f′(y)y/f(y)>0. All those confirm that f′fyf′(y)y<-f′(y)y/f(y)f(f(y))<f(f(y)); namely,(22)f′fyf′yy<ffy.
Denoteg(y)=f(f(y))/y; then we have g′(y)=f′fyf′(y)y-f(f(y))/y2. Then, g′(y)<0 and g(y) is a monotonically decreasing function at interval y∈(y∗,y2] due to (22).
Therefore,g(y∗)=f(f(y∗))/y∗=f(y∗)/y∗=y∗/y∗=1, and g(y)<1 for y∈(y∗,y2]; that is, f(f(y))<y.
Based on the above discussions we conclude thatf2k(y) is a monotonically decreasing function as k increases with f2k(y)→y∗(k→∞). Similarly, we can prove that f2k+1(y) is a monotonically decreasing function as k increases with f2k(y)→y∗(k→∞). Thus, fn(y)→y∗(n→∞) for any y∈[0,∞) holds.
Combined with all cases together we know that the unique equilibriumy∗ is globally stable if it exists, and consequently the periodic solution yp(t) of model (4) is globally stable. This completes the proof.The stability of the positive periodic solution is determined by the stability of the positive equilibrium of the difference equation made by sequence of impulsive points; thus the positive periodic solutionyp(t) of model (4) is globally stable according to Theorem 1. Therefore, a general expression of unique pest-free periodic solution of model (3) over the interval nT<t≤(n+1)T for all n∈N can be denoted by(23)0,ypt=0,y∗e-dt-nT.In the following, we focus on the stability of the pest-free periodic solution for model (3). Firstly, define(24)x¯t,y¯t=xt,yt-0,ypt,and then (x¯(t),y¯(t)) represents the deviations form of the pest-free solution. According to model (3), the derivations of (x¯(t),y¯(t)) are (25)x¯˙=rx¯-rKx¯2-ax¯1+ωx¯y¯+ypt,y¯˙=cx¯1+ωx¯y¯+ypt-dy¯.Assuming that (x¯(t),y¯(t)) are small enough, then the linear approximation of the deviation system around the periodic solution (0,yp(t)) can be obtained as follows:(26)x¯˙=r-ayptx¯,y¯˙=cyptx¯-dy¯.Then, we have the following main result.Theorem 2.
The pest-free periodic solution(0,yp(t)) of model (3) is locally stable in the first quadrant provided that (27)R1≐rdTay∗1-e-dT<1.Proof.
LetΦ(t) be the fundamental matrix of (26); then (28)ΦT=e∫0Tr-ayptdt0∗e-dT,and Φ(0)=I is the identity matrix. The term ∗ is not required in the following analysis.
The resetting impulsive conditions of (3) become (29)x¯nT+y¯nT+=1001-τθ1+θypnT2x¯nTy¯nT=BnTx¯nTy¯nT.Thus, the Floquet multipliers of the matrix (30)M=BTΦTare(31)λ1=e∫0Tr-ayptdt,λ2=1-τθ1+θypT2e-dT.It follows from the literature [26] that the pest-free solution is locally stable if and only if |λ1|<1 and |λ2|<1. The former inequality is equal to R1≐rdT/ay∗(1-e-dT)<1 and the latter inequality is true if and only if R2≐1-τθ/(1+θy∗e-dT)2e-dT<1, which is true naturally based on the analysis in the previous section. Therefore, the pest-free periodic solution is locally stable when R1<1.


## 2.2. Boundness and Permanence of System (3)

According to the proof of Theorem2, the pest-free periodic solution is unstable when R1>1. So, next we aim to prove the permanence of the system if R1>1. First we prove the boundness of system (3).Theorem 3.
There exists a constantM0>0 such that x(t)≤M0 and y(t)≤M0 for each solution of system (3) with all t large enough.Proof.
Supposex(t) and y(t) are any solution of system (3). Let V(t)=(c/a)x(t)+y(t). Then, when t≠nT, (32)dVtdt+λVt=car+λx-craKx2-λ-dy.When 0<λ≤d, the right hand of the equation is bounded; namely, (33)dVtdt+λVt≤M1,with M1>0. When t=nT, (34)Vt+=Vt+τ1+θyt-caδxth+xt≤Vt+τ.
Considering the system(35)dVtdt≤-λVt+M1,t≠nT,Vt+≤Vt+τ,t=nT,we have (36)Vt≤V0+-M1λe-λt+τ1-e-nλT1-e-λTe-λt-nT+M1λ.So, V(t) is uniformly ultimately bounded. Hence, there exists a constant M0>0 such that x(t)≤M0 and y(t)≤M0 for each solution of system (3) with all t large enough. This completes the proof.Theorem 4.
System (3) is permanent provided R1>1.Proof.
It follows from Theorem1 that there exists an m1>0 such that y≥m1 for t large enough. Thus, in order to show the permanence of model (3) we shall next find a constant m2>0, such that x≥m2 for t large enough.
Considering the second equation of system (3), we have(37)dytdt≤-dyt,t≠nT,yt+=yt+τ1+θyt,t=nT.Therefore, according to the comparison theorem and the analyses in the previous subsection, we have (38)yt≤ypt+ε1for t>t1 and ε1>0. Then, by the first equation of system (3) we have(39)dxtdt≥rxt1-xtK-aypt+ε1xt,t≠nT,xt+=1-δxtxt+hxt,t=nT.Considering the auxiliary system(40)dutdt=rut1-utK-aypt+ε1ut,t≠nT,ut+=1-δutut+hut,t=nTand using the same method as system (4), we have (41)un+1+=un+11-δun+1un+1+hwith(42)un+1=un+e∫0Tr-aypt+ε1dt1+un+∫0Tr/Ke∫0tr-aypξ+ε1dξdt.Let M1=e∫0Tr-aypt+ε1dt and M2=∫0Tr/Ke∫0t(r-a(yp(ξ)+ε1))dξdt; then un+1=un+M1/1+M2un+. Therefore, we have the following different equation with respect to un+:(43)un+1+=un+M11+M2un+1-δun+M1/1+M2un+h+un+M1/1+M2un+≐Fun+.Let u∗ be the steady state of the above difference equation, and by simple calculation we have (44)a1u∗2+a2u∗+a3=0,where a1=M22+M2M1, a2=M2+M1+M2h-M2M1-M12+δM12, and a3=h(1-M1). Obviously, a1>0. Moreover, a3<0 when R1>1 for ε1 sufficiently small, which shows that the above difference equation exists with a unique positive equilibrium u∗. Therefore, system (40) has only one positive periodic solution when R1>1, denoted by u¯(t).
To study the stability of this periodic solution, we first investigate the properties of the functionF(u)=M1u/1+M2u(1-δM1u/h(1+M2u)+M1u). Taking the derivation of F(u) with respect to u, one yields (45)F′u=M1hM2+M1u+h2-2δM12hM2+δM23u2-2δM12hu1+M2u2hM2+M1u+h2.Let F′(u)=0; then we get the following equation: (46)b1u2+b2u+b3=0,where b1=M1M22h2+2(1-δ)hM12M2+(1-δ)M13,b2=2h2M1M2+2(1-δ)hM12, and b3=h2M1. It follows from 0<δ<1 that b1>0,b2>0, and b3>0. Moreover, we can prove that b22-4b1b3=4M14h2δ(1-δ)>0. Therefore, the equation F′(u)=0 does not exist with any positive root. Thus, F′(u)>0 for all u>0; that is, the function F(u) is a monotonically increasing function for all u>0. Further, it is easy to show that F′(0)=M1 according to the expression of F′(u) with M1>1 for R1>1 and ε1 sufficiently small. All those confirm that the intersection angle of the function F(u) and the horizontal coordinate is larger than π/4. Therefore, the function image is as shown in Figure 1. Then, according to u<F(u)<u∗ we have the fact that Fk(u) is monotonically increasing as k increases for any u∈[0,u∗), and Fk(u)→u∗(k→∞); according to F(u)<u we have the fact that Fk(u) is monotonically decreasing as k increases for any u>u∗, and Fk(u)→u∗(k→∞). Therefore, the equilibrium u∗ is globally stable.
Thus, the positive periodic solutionu¯(t) is globally stable. Then, according to the comparison theory, x(t)≥u(t)≥u¯(t)-ε2 for t≥t2. Let m2=min(u¯(t))-ε2, ε2>0; then we have x≥m2 for t large enough. Thus, system (3) is permanent provided R1>1. This completes the proof.Figure 1
Illustration of the global stability ofu∗.


## 2.3. Bifurcation of Positive Periodic Solution

In this section, we investigate the existence and stability of positive periodic solution of system (3) by employing the bifurcation analysis related to the pest-free periodic solution X0=(0,yp(t)) [26].Theorem 5.
As the parameterT passes through the critical value T0 satisfying e∫0T0[r-ayp(ξ)]dξ=1, a nontrivial positive periodic solution appears near X0. The bifurcation is supercritical provided BC<0, while it is subcritical if BC>0, where (47)B=-r-aypT0+adT0y∗T0C2e-dT01-C2e-dT0e∫0T0r-aypξdξ,(48)C=2δh-1e2∫0T0r-aypξdξ-2aT01-C1C21-C2e-dT0e∫0T0r-aypξdξ-C3.Proof.
To use the theory of [26], we employ the same notations as those in [26]; then (49)F1x,y=cxy1+ωx-dy,F2x,y=rx1-xK-axy1+ωx,Θ1x,y=y+τ1+θy,Θ2x,y=1-δxx+hx.By simple calculations, we get (50)∂Φ1T0,X0∂y=e-dT0>0,∂Φ2T0,X0∂x=e∫0T0r-aypξdξ>0,∂Φ1T0,X0∂x=c∫0T0ypue-dT0-ue∫0ur-aypξdξdu≐C1>0,a0′=1-C2e-dT0,d0′=1-e∫0T0r-aypξdξ,b0′=1-C1C2,where C2=1-τθ/(1+θyp(T0))2 and T0 is the root of d0′=0, satisfying the condition cy∗(e-d(T0-nT)-1)/adT0=1. Obviously, we have the fact that C2<1, which shows that a0′>0. Further, we have (51)∂Φ1T0,X0∂T¯=-dy∗T0e-dT0,∂2Φ2T0,X0∂x∂y=-aT0e∫0T0r-aypξdξ<0,∂2Φ2T0,X0∂x∂T¯=r-aypT0e∫0T0r-aypξdξ,∂2Φ2T0,X0∂x2=T0-2rK+2aωypT0e∫0T0r-aypξdξ-∫0T0ae∫uT0r-aypξdξ∫0ue-du-vcypve∫0T0r-aypξdξdvdu≐C3.∂Θ1(T0,X0)/∂y=C2, ∂Θ1(T0,X0)/∂x=0, ∂Θ2/∂x=1, ∂Θ2(T0,X0)/∂y=0, ∂2Θ1(T0,X0)/∂x∂y=-2δh-1, and ∂2Θ1(T0,X0)/∂x∂y=0. Thus, (52)B=-r-aypT0+adT0y∗T0C2e-dT01-C2e-dT0e∫0T0r-aypξdξ,C=2δh-1e2∫0T0r-aypξdξ-2aT01-C1C21-C2e-dT0e∫0T0r-aypξdξ-C3.Therefore, according to Theorem 2 in literature [26], the supercritical branch occurs at the point T0 if BC<0 and the subcritical branch occurs if BC>0. This completes the proof.


## 3. Numerical Investigations

Bifurcation analyses shown above indicate that system (3) could exist with complex dynamics once the pest-free periodic solution becomes unstable, which is a challenge for theoretical studies. Therefore, in this section we will focus on the complex dynamics by using numerical techniques.


### 3.1. Backward and Forward Bifurcation

Our theoretical analysis shows that whenR1>1 system (3) is permanent and when R1<1 the pest-free periodic solution is locally stable in the first quadrant. Moreover, Theorem 5 indicates the existence of nontrivial positive periodic solution under certain conditions. Thus, to verify these theoretical results, we first take some numerical simulations under certain parameter spaces and initial values. Figure 2 shows the bifurcation diagram of system (3) with respect to the threshold condition R1 with d=0.6 and d=0.8.Figure 2
Backward (forward) bifurcation of system (3) with respect to threshold value R1. The parameter values are as follows: T=6,τ=30,θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,c=0.23, and d=0.6 in (a) and d=0.8 in (b).
(a)
(b)It follows from Figure2(a) that when R1<1 there exists a stable pest-free periodic solution, and when R1>1 the pest-free periodic solution becomes unstable and a stable positive periodic solution appears, which is called a forward bifurcation at R1=1. Figure 2(b) indicates that when R1<0.78(T<0.406) there exists only one pest-free periodic solution and it is stable. When R1>1, there exists at least one stable positive periodic solution. Specifically, when 0.78<R1<1, the pest-free periodic and the positive periodic solutions can coexist which are bistable; namely, when the value of R1 passes through 1 the supercritical bifurcation occurs or backward bifurcation occurs. Both of these two figures confirm that when R1<1 the pest-free periodic solution is locally stable and when R1>1 system (3) is permanent. Moreover, there may also exist a locally stable positive periodic solution when R1<1, and this interesting result could be produced by the density dependent pulse control strategies.As shown in Figure2(b), there may be two locally stable periodic solutions when R1<1. Therefore, the initial values of both the pest and predator populations are crucial. To see how different initial values affect the successful pest control, we plot the basins of attraction of two locally stable periodic solutions of system (3), as shown in Figure 2(b). It follows from Figure 3 that the pest can be eradicated when the initial values of the predator are very small or large. However, when the initial values of the predator have been chosen between 2 and 30, the pest can be eradicated only if the initial values of the pest are very small.Figure 3
Basins of attraction of system (3). The parameter values are as follows: T=6,r=1.5;a=0.3;K=50;ω=0.3;c=0.23;d=0.8;τ=30;h=0.8;δ=0.85;θ=0.5.The theoretical results only give the existence and global stability of the pest-free periodic solution or local stability of the positive periodic solution. To further study the dynamics of system (3), the bifurcation diagrams for the pest population with respect to bifurcation parameters T, τ, and δ are displayed in Figures 4, 5, and 6. The resulting bifurcation diagrams clearly show rich dynamics of system (3). Specifically, Figure 5 shows the existence of period-halving bifurcations, and Figure 4 depicts the existence of period-doubling bifurcations, while Figure 6 indicates that several types of attractor may coexist. All those complex dynamics may result in complexity of pest control.Figure 4
Bifurcation diagram of system (3) with respect to T. The parameter values are as follows: τ=1,θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.Figure 5
Bifurcation diagram of system (3) with respect to parameter τ. The other parameter values are as follows: T=6θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.Figure 6
Bifurcation diagram of system (3) with respect to parameter δ. The other parameter values are as follows: T=6,τ=1,θ=0.5,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.


### 3.2. Effects of Density Dependent Releasing Strategy

According to the main results shown in Theorem4, system (3) is permanent provided R1>1. Thus, to control the pest population, we need to know how the parameters affect the value of the threshold R1. Obviously, the threshold condition is determined by the key parameters including the pulse frequency, the maximum fatality rate, and the half-saturation constant. Figure 7 shows the relationship between the threshold R1 and the pulse period T. It follows from this figure that the value of R1 could increase as pulse period T increases, as shown in Figure 7(a). However, it also can decrease as T increases, as shown in Figure 7(b). All those indicate that the relations between the threshold R1 and pulse period T are quite complex, which depicts the threshold condition R1 is nonmonotonic with respect to T. Thus, reducing the pulse period T is not always beneficial for pest control. That is, frequent releasing of natural enemies may not be beneficial to the eradication of pests when the amount of predator released depends on the predator population, which is one of the novel results obtained in the present work compared with the previous results obtained on this topic [3–5, 13, 27].Figure 7
(a) Relationship between the thresholdR1 and the pulse period T with θ=0.5 (the dashed line) and θ=4 (the solid line). (b) Behavior of the prey when T=3,6, and 9 for the case θ=4. The values of R1 when T=3,6, and 9 are 1.0928, 0.9245, and 1.1056, respectively. Other parameter values are r=0.8,a=0.3,τ=20,d=0.9.
(a)
(b)The releasing term in model (3) is formulated by τ/1+θy, which is a predator density dependent function and depicts how the number of natural enemy populations released was guided by its current density at the fixed moment. The parameter θ is a shape parameter which indicates the dependence of releasing amount and the predator population. Thus, to investigate the effects of the predator population guided releasing method on the pest control, the variations of the threshold R1 with respect to parameters T, θ, and τ have been addressed, as shown in Figure 8.Figure 8
(a) Contour plots ofR1 with respect to θ and τ for T=3. (b) Contour plots of R1 with respect to τ and T when θ=4. (c) Contour plots of R1 with respect to θ and T when τ=20. The other parameters are fixed as r=0.8, d=0.9, and a=0.3.
(a)
(b)
(c)Figure8(a) shows the contour plots of R1 with respect to θ and τ. It follows from this figure that increasing τ can lead to an obvious decrease of R1, no matter how much the value of θ is, while decreasing θ can also lead to a slight decrease of R1 and the effect is more evident when the value of τ is large. Figure 8(b) shows the contour plots of R1 with respect to T and τ. This figure indicates that both T and τ together can affect the threshold condition R1 significantly. Note that for relatively small τ the threshold condition R1 is a decreasing function as T decreases. However, for a relatively large τ, decreasing period T may first lead to a complex effect on R1; that is, R1 is nonmonotonic with respect to T, which is similar to the important results shown in Figure 7.Similarly, this type of interesting results can also be found in Figure8(c) for the contour plots of R1 with respect to θ and T. Figure 8(c) indicates that if T is small, decreasing θ can result in an obvious decrease of R1, while if T is relatively large (larger than 8), decreasing θ almost has no effect on R1. Thus, Figure 8 reveals that the key control factors including the pulse frequency, pulse amount, and the dependence of pulse amount and the current amount of the predator could affect each other. For example, when the control measure is implemented frequently (i.e., T is small) or the releasing amount is large (i.e., τ is large), effect of the dependence between the pulse amount and the current amount of the predator is very small. Moreover, the complex relationships between T and R1 usually appear when τ and θ are relatively large. All those confirm the control measures should be carefully designed in practice. Otherwise, the paradox and hormesis could occur when pesticides are applied in the field.


## 3.1. Backward and Forward Bifurcation

Our theoretical analysis shows that whenR1>1 system (3) is permanent and when R1<1 the pest-free periodic solution is locally stable in the first quadrant. Moreover, Theorem 5 indicates the existence of nontrivial positive periodic solution under certain conditions. Thus, to verify these theoretical results, we first take some numerical simulations under certain parameter spaces and initial values. Figure 2 shows the bifurcation diagram of system (3) with respect to the threshold condition R1 with d=0.6 and d=0.8.Figure 2
Backward (forward) bifurcation of system (3) with respect to threshold value R1. The parameter values are as follows: T=6,τ=30,θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,c=0.23, and d=0.6 in (a) and d=0.8 in (b).
(a)
(b)It follows from Figure2(a) that when R1<1 there exists a stable pest-free periodic solution, and when R1>1 the pest-free periodic solution becomes unstable and a stable positive periodic solution appears, which is called a forward bifurcation at R1=1. Figure 2(b) indicates that when R1<0.78(T<0.406) there exists only one pest-free periodic solution and it is stable. When R1>1, there exists at least one stable positive periodic solution. Specifically, when 0.78<R1<1, the pest-free periodic and the positive periodic solutions can coexist which are bistable; namely, when the value of R1 passes through 1 the supercritical bifurcation occurs or backward bifurcation occurs. Both of these two figures confirm that when R1<1 the pest-free periodic solution is locally stable and when R1>1 system (3) is permanent. Moreover, there may also exist a locally stable positive periodic solution when R1<1, and this interesting result could be produced by the density dependent pulse control strategies.As shown in Figure2(b), there may be two locally stable periodic solutions when R1<1. Therefore, the initial values of both the pest and predator populations are crucial. To see how different initial values affect the successful pest control, we plot the basins of attraction of two locally stable periodic solutions of system (3), as shown in Figure 2(b). It follows from Figure 3 that the pest can be eradicated when the initial values of the predator are very small or large. However, when the initial values of the predator have been chosen between 2 and 30, the pest can be eradicated only if the initial values of the pest are very small.Figure 3
Basins of attraction of system (3). The parameter values are as follows: T=6,r=1.5;a=0.3;K=50;ω=0.3;c=0.23;d=0.8;τ=30;h=0.8;δ=0.85;θ=0.5.The theoretical results only give the existence and global stability of the pest-free periodic solution or local stability of the positive periodic solution. To further study the dynamics of system (3), the bifurcation diagrams for the pest population with respect to bifurcation parameters T, τ, and δ are displayed in Figures 4, 5, and 6. The resulting bifurcation diagrams clearly show rich dynamics of system (3). Specifically, Figure 5 shows the existence of period-halving bifurcations, and Figure 4 depicts the existence of period-doubling bifurcations, while Figure 6 indicates that several types of attractor may coexist. All those complex dynamics may result in complexity of pest control.Figure 4
Bifurcation diagram of system (3) with respect to T. The parameter values are as follows: τ=1,θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.Figure 5
Bifurcation diagram of system (3) with respect to parameter τ. The other parameter values are as follows: T=6θ=0.5,δ=0.85,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.Figure 6
Bifurcation diagram of system (3) with respect to parameter δ. The other parameter values are as follows: T=6,τ=1,θ=0.5,h=0.8,r=1.5,a=0.3,K=50,ω=0.3,d=0.8,c=0.5.


## 3.2. Effects of Density Dependent Releasing Strategy

According to the main results shown in Theorem4, system (3) is permanent provided R1>1. Thus, to control the pest population, we need to know how the parameters affect the value of the threshold R1. Obviously, the threshold condition is determined by the key parameters including the pulse frequency, the maximum fatality rate, and the half-saturation constant. Figure 7 shows the relationship between the threshold R1 and the pulse period T. It follows from this figure that the value of R1 could increase as pulse period T increases, as shown in Figure 7(a). However, it also can decrease as T increases, as shown in Figure 7(b). All those indicate that the relations between the threshold R1 and pulse period T are quite complex, which depicts the threshold condition R1 is nonmonotonic with respect to T. Thus, reducing the pulse period T is not always beneficial for pest control. That is, frequent releasing of natural enemies may not be beneficial to the eradication of pests when the amount of predator released depends on the predator population, which is one of the novel results obtained in the present work compared with the previous results obtained on this topic [3–5, 13, 27].Figure 7
(a) Relationship between the thresholdR1 and the pulse period T with θ=0.5 (the dashed line) and θ=4 (the solid line). (b) Behavior of the prey when T=3,6, and 9 for the case θ=4. The values of R1 when T=3,6, and 9 are 1.0928, 0.9245, and 1.1056, respectively. Other parameter values are r=0.8,a=0.3,τ=20,d=0.9.
(a)
(b)The releasing term in model (3) is formulated by τ/1+θy, which is a predator density dependent function and depicts how the number of natural enemy populations released was guided by its current density at the fixed moment. The parameter θ is a shape parameter which indicates the dependence of releasing amount and the predator population. Thus, to investigate the effects of the predator population guided releasing method on the pest control, the variations of the threshold R1 with respect to parameters T, θ, and τ have been addressed, as shown in Figure 8.Figure 8
(a) Contour plots ofR1 with respect to θ and τ for T=3. (b) Contour plots of R1 with respect to τ and T when θ=4. (c) Contour plots of R1 with respect to θ and T when τ=20. The other parameters are fixed as r=0.8, d=0.9, and a=0.3.
(a)
(b)
(c)Figure8(a) shows the contour plots of R1 with respect to θ and τ. It follows from this figure that increasing τ can lead to an obvious decrease of R1, no matter how much the value of θ is, while decreasing θ can also lead to a slight decrease of R1 and the effect is more evident when the value of τ is large. Figure 8(b) shows the contour plots of R1 with respect to T and τ. This figure indicates that both T and τ together can affect the threshold condition R1 significantly. Note that for relatively small τ the threshold condition R1 is a decreasing function as T decreases. However, for a relatively large τ, decreasing period T may first lead to a complex effect on R1; that is, R1 is nonmonotonic with respect to T, which is similar to the important results shown in Figure 7.Similarly, this type of interesting results can also be found in Figure8(c) for the contour plots of R1 with respect to θ and T. Figure 8(c) indicates that if T is small, decreasing θ can result in an obvious decrease of R1, while if T is relatively large (larger than 8), decreasing θ almost has no effect on R1. Thus, Figure 8 reveals that the key control factors including the pulse frequency, pulse amount, and the dependence of pulse amount and the current amount of the predator could affect each other. For example, when the control measure is implemented frequently (i.e., T is small) or the releasing amount is large (i.e., τ is large), effect of the dependence between the pulse amount and the current amount of the predator is very small. Moreover, the complex relationships between T and R1 usually appear when τ and θ are relatively large. All those confirm the control measures should be carefully designed in practice. Otherwise, the paradox and hormesis could occur when pesticides are applied in the field.


## 4. Discussion

Considering the effects of resources limitation and pest management purpose, we proposed a Holling type II pest-natural enemy model with density guided releasing strategy and killing rate, which can provide a more nature understanding and a better portrayal of the reality. However, the nonlinear pulse functions not only result in difficulties for dynamic analyses, but also can reveal some important results related to pest control compared with those obtained in the previous papers [3, 4, 9, 22, 23].Theoretical results show that the pest-free periodic solution is locally stable whenR1<1 and the system is permanent when R1>1. Besides, theoretical analysis also indicates that when T passes through the critical value T0, a nontrivial positive periodic solution appears near the pest-free periodic solution. Further, extensive numerical investigations imply that there may be a locally stable positive periodic solution even if R1<1. This results in a backward bifurcation and is quite different from the previous works [3, 17], in which the pulse functions are linear. This also indicates that the pest may not be eradicated even if R1 is controlled to be below 1. All those results demonstrate that the dynamic behavior of model (3) may be dramatically affected by small changes in the value of initial densities of pest and natural enemy under resource limitation. Moreover, as shown in Figure 3, small or large initial values of the natural enemies may lead to the eradication of the pest. Medium level of the initial values of the natural enemies usually can not suppress the pest population and often lead to the outbreak of the pest.The threshold conditionR1 depends on all key parameters related to the pest control, such as releasing period, maximum releasing constant, and shape parameters of the nonlinear functions. To analyze the effects of those parameters, we first investigate the relationship between R1 and the pulse period. Our results show that R1 may be a nonmonotonic function with respect to the pulse period T. This indicates that shortening the pulse period may lead to an increase of R1 and consequently may cause the outbreak of the pest, as shown in Figure 7. This result is not the same as the traditional results obtained in the previous works [3, 4, 9, 22, 23]. This may be because that frequent release leads to large number of current predators. As the number of predators released is guided by its current density, the large number of the current predators results in a small amount of releasing every time. So, under this circumstance, frequent release of the predator may be beneficial or detrimental to the eradication of pests.The releasing number and killing rate guided by both the pest and natural enemy populations have been introduced into the model; one of the main purposes is to analyze the effect of this kind of nonlinear control methods in terms of reducing the pest population. So, it is essential to suggest how to choose the shape parameterθ and the releasing rate τ and to determine how the releasing period affects the control output. Numerical results as shown in Figure 8 reveal that when releasing rate is relatively large, reducing the dependence of releasing amount and the predator population, namely, reducing the value of θ, has little effect on the reduction of the pest population, while it has large effect on the reduction of the population of the pest when the releasing rate is small. Moreover, frequent releasing of the predator can only lead to decreasing of R1 when the dependence of releasing amount and the predator population is small and frequent releasing of the predator may result in an increase of R1 when the dependence of releasing amount and the predator population is large.The results presented in this paper indicate that the successful biological control depends on the initial densities of the pest and predator and the reasonable implementing plan of control measures. The releasing of predators guided by the predator population is more sensitive to the pulse period and releasing amount every time. However, some simplifications of the real prey-predator interactions are neglected, such as predator hiding and effect of pesticide on the predator. In fact, the prey refuge can be broadly defined including spatial or temporal refuges, prey aggregations, or reduced prey activity [28–31]. It has been shown in many previous papers that the existence of refuges may have important effects on the coexistence of the prey and predator. Therefore, more realistic mathematical model and its dynamic analyses should be studied in the future which can help us to design control strategies and to make management decisions.



---

*Source: 8683207-2017-06-14.xml*