# Modelling and Analysis of a Host-Parasitoid Impulsive Ecosystem under Resource Limitation

Authors: Peipei Wang; Wenjie Qin; Guangyao Tang
Journal: Complexity
 (2019)

---

## Abstract

With a long history of theoretical development, biological model has focused on the interaction of a parasitoid and its host. In this paper, two Nicholson-Bailey models with a nonlinear pulse control strategy are proposed and analyzed to examine how limited resource affects the pest control. For a fixed-time discrete impulsive model, the existence and stability of the host-free periodic solution are derived. Threshold analysis suggests that it is critical to release parasitoid in an optimal number in case of the happening of the intra-specific competition, which will seriously affect the pest control. Bifurcation analysis reveals that the model exists complex dynamics including period doubling, chaotic solutions, coexistence of multiple attractors, and so on. For a state-dependent discrete impulsive model, the numerical simulations for bifurcation analysis are studied, the results show that how the key parameters and the initial densities of both populations affect the pest outbreaks, and consequently the relative biological implications with respect to pest control are discussed.


---

## Body

## 1. Introduction

Mathematical models can assist in the design and understanding of basic problems in biology, medicine, and life sciences, which can take many forms including dynamical system, differential equation, difference equation, and so on [1–3]. Different models can describe different biological phenomena, with each having different strengths. Some certain species, including many species of insect, have no overlap between successive generations and so their populations evolve in discrete time-steps. The discrete model governed by a difference equation is more appropriate than the continuous ones taking into account the fact that the population has a short life expectancy, non-overlapping generations in the real world [4–6], so it is reasonable to study biological models governed by difference equations [7–13].The insect pests’ outbreaks often cause serious ecological and economic problems, requiring complex control measures to reduce the harm brought by insect pests of agriculture and insect vectors of important plant, animal, and human diseases. In order to control pest more effective, at the same time, environmental influences and human interventions are taken into consideration, and the concept of Integrated Pest Management (IPM) strategy has been proposed [14–17], as a consequence it is challenging to evaluate the effectiveness of control measures such as biological (releasing predators), chemical (spraying pesticides), and physical control tactics. Based on the concept of IPM strategy, some mathematical models to describe the interaction between the pest and its natural enemy have been developed [18–21].One of the IPM strategies is to release natural enemies and spray chemical pesticides at a fixed time, since pesticides can quickly kill a significant portion of a pest population, and it is convenient to carry out. Meanwhile, the IPM strategy can be applied more effectively to achieve maximum effect at the minimum levels of pesticide if the life cycle of a pest is understood; i.e., it is used when the pest is at its most vulnerable. Pesticides can be relatively cheap and easy to apply and fast-acting and in most instances can be relied on to control the pests. Therefore, the periodic control strategy has been applied widely and extended in agriculture and ecology [18, 19].Another important concept in IPM strategy is the Threshold Policy Control (TPC) [22, 23], which maintains the pests’ density below the Economic Injury Level (EIL) [24–26] by releasing predators and spraying pesticides once the pests’ density reaches the Economic Threshold (ET) [27, 28]; it is also called by state-dependent feedback control. As we know, the main purpose of IPM is to maintain the density of the pests below the EIL rather than seeking to eradicate them, and the suitable tactic can be applied only when the density of pests reaches the given ET; it can minimize the damage of insecticides to non-target pests and preserve the quality of the environment. Therefore the state-dependent impulsive dynamical model could be appropriate to describe such control tactics [20, 21].However, the above IPM strategy has assumed that the killing rate of pesticide is a constant, which shows that the agricultural resources including chemical pesticides, labor forces, facilities, and costs are very effective and sufficient for pest control. In reality, every community or country has an appropriate or limited capacity for pesticides, costs, etc., especially for developing countries. Although it is critical to understand resource limitation to effective management and conservation of populations, limited resource is difficult to quantify partly because it is a dynamic process [29–32]. The main purpose of this paper is to describe the limited resources based on the classic Nicholson-Bailey model [33]; there are two different models with IPM strategies being derived. In order to characterize the saturation phenomenon of the limited resources, the constant fatality rate is replaced with a nonlinear saturation function which depends on its population density. Some theoretical, numerical, biological analyses are given to investigate how limited resources affect host outbreaks.The organization of the present paper is as follows: in the next section, a fixed-time impulsive control for a Nichloson-Bailey model with resource limitation is proposed; by employing the qualitative analysis and the discrete dynamical system determined by the stroboscopic map, the threshold conditions which guarantee the existence and stability of the host-free periodic solution are obtained. Meanwhile, 1-dimensional bifurcation analyses are investigated; the key parameters and initial values of both populations affect the pest outbreaks are given by taking advantage of numerical simulation. In Section3, we describe a Nichloson-Bailey ecosystem with resource limitation by state-dependent impulsive control; by employing numerical simulations, the model is analyzed to understand how resource limitation affects pest population outbreaks; the intersection between the initial densities and pest control will be discussed. The paper ends with some interesting biological conclusions, which complement the theoretical findings.


## 2. Fixed-Time Impulsive Control for a Nicholson-Bailey Model with Resource Limitation

### 2.1. Model Formulation

In 1935, Nicholson and Bailey [33] proposed a classic discrete host-parasitoid model, as an intergenerational survival rate for parasitoids was considered. In 2008, Tang [8] extended the basic Nicholson-Bailey model as follows: (1)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn.Here Hn and Pn represent the density of hosts(or pests) and parasitoids(or natural enemies) at n-th generation, respectively; r is the intrinsic growth rate for the host population; a denotes the searching efficiency of the parasitoid population; the terms exp⁡(-aPn) and [1-exp⁡(-aPn)] represent the probability that the host succeeds and fails in escaping from the parasitism, respectively; σ is the density-independent survival of the parasitoid at n-th generation, n∈Z≜{0,1,2,⋯}.In order to investigate the dynamics of system (1) with IPM strategies, by introducing periodic spraying pesticide only for hosts, and releasing parasitoids, Tang [8] proposed the following Nicholson-Bailey model: (2)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn,n∈ZHqk+=1-q1Hqk,Pqk+=1+q2Pqk+τ,q,k∈Nwith the initial densities (H0+,P0+)=(H0,P0) and N≜{1,2,⋯}. Meanwhile, q∈N is the pulse period, q1 is the killing rate for hosts, and q2 and τ are the release rate and constant for parasitoids at qk-th generation, respectively. Hqk and Pqk (Hqk+ and Pqk+) denote the densities of both populations at qk-th generation before (after) the impulsive perturbations, respectively.Now, taking into account the resource limitation and saturation effects, that is, a linear pulse in model (2) is replaced with a saturation phenomenon of limited resources; especially, a Hill function [34] is used to describe the resource limitation. For convenience and simplification, model (2) can be rewritten as (3)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn,n∈Z,Hqk+=1-q1maxHqkθ1+HqkHqk,Pqk+=1+q2maxPqkθ2+PqkPqk+τ,q,k∈N,where q1max and θ1 represent the maximal fatality rate and the half-saturation constant for hosts, respectively. q2max and θ2 represent the maximal release rate and the half-saturation constant for parasitoids, respectively. Assume that 0≤qimax<1 and θi>0 for i=1,2.


### 2.2. Mathematical Analysis of the Host-Free Periodic Solution and Threshold Conditions

With the main purpose of IPM strategy, it is necessary to illustrate the existence of the host-free periodic solution of model (3) and focus on determining its global attractivity. In order to do this, we first investigate the host-free set: A={(H,P)∈R+2,H=0}, which is clearly invariant by model (3); it follows from model (3) and its pest-free set A that (4)Pn+1=σPn,n∈Z,Pqk+=1+q2maxPqkθ2+PqkPqk+τ,k∈N,P0+=P0.Obviously, model (4) is a periodic system, so the periodic solution Pn can be defined at the subinterval [qs+,q(s+1)) with the initial value Pqs+, s∈Z. Noting that Pqs+ represents the parasitoid’s density after the impulsive perturbation.From the first equation of model (4), we can use this iteration equation at interval [qs+,q(s+1)) and obtain Pn+1=σn+1-qsPqs+. Especially, it gives that Pq(s+1)=σqPqs+ as n=q(s+1)-1.Now, combining the second equation of model (4) and the above one yields (5)Pqs+1+=q2max+1σ2qPqs+2+θ2+τσqPqs++θ2τθ2+σqPqs+.Denoting Pq(s+1)+ as Pn+1, then the above equation can be rewritten as the following difference equation: (6)Pn+1=FPn≜1+q2maxσ2qPn2+θ2+τσqPn+θ2τθ2+σqPn.Equation (6) is a stroboscopic map of model (4): it describes the relations of the parasitoids’ density between any two successive pulse points; that is, the existence of a positive steady state of (6) represents the existence of a positive periodic solution of model (4). Therefore, we first explore the positive steady state of the stroboscopic map (6); let us get the derivative of F(Pn) with Pn and yield (7)F′Pn=q2max+1σ3qPn2+2q2max+1σ2qθ2Pn+θ22σqσ2qPn2+2θ2σqPn+θ22.According to the stability theory of differential equation, it is necessary to ensure that the inequality |F′(Pn)|<1 holds true, and this condition is equivalent to (8)A0Pn2+B0Pn+C0<0,where A0=σ2q[(1+q2max)σq-1], B0=2θ2σq[(1+q2max)σq-1], C0=θ22(σq-1). By a simple calculation, one yields Δ0=B02-4A0C0=4θ22q2maxσqA0. It is easy to find that if A0<0 holds true, then Δ0<0, which implies |F′(Pn)|<1. For simplification, A0<0 can be rewritten as (9)R1≜1+q2maxσq<1.The above inequality can ensure the stability of a positive solution of the stroboscopic map (6). Meanwhile, it also plays an important role in discussing the existence of the positive fixed point for (6).Next, we continue to investigate the positive fixed point of the stroboscopic map (6). For convenience, we denote this fixed point by P~; then it follows from the stroboscopic map (6) that (10)P~=1+q2maxσ2qP~2+θ2+τσqP~+θ2τθ2+σqP~,which means (11)A1P~2+B1P~+C1=0,where A1=σq[1+q2maxσq-1], B1=(θ2+τ)σq-θ2, C1=θ2τ>0. Obviously, we can obtain that A1<0 and Δ1=B12-4A1C1>0 as R1<1 holds true. That is, the stroboscopic map (6) exists in only one positive fixed point P~=-(B1+Δ1)/(2A1). Therefore, we have the following result.Lemma 1.
Model (4) exists as a positive periodic solution Pn∗ provided that R1<1 holds true. Meanwhile, for every positive solution Pn of model (4), it is given that limn→+∞Pn-Pn∗=0. Here Pn∗=σn-qsP~ with n∈[qs+,q(s+1)),s∈Z.Based on Lemma1, we can obtain the general expression of the host-free periodic solution of model (3) in the interval [qs+,q(s+1)) for all s∈Z denoted by (12)0,Pn∗=0,-σn-qsB1+Δ12A1.In order to investigate the global stability of the host-free periodic solution (0,Pn∗) of model (3), the following theorem is derived.Theorem 2.
Let(Hn,Pn) be any positive solution of model (3); the host-free periodic solution (0,Pn∗) of model (3) is globally asymptotically stable in the first quadrant, provided that (13)R2≜2qrA11-σa1-σqB1+Δ1<1.Proof.
To prove Theorem2, the following condition (16) is considered first. It follows from R2<1, model (3), and Pn∗=σn-qsP~ that we have (14)∑n=qsqs+1-1r-aPn∗<0,which implies (15)∏n=qsqs+1-1exp⁡r-aPn∗<1.Therefore, from the above inequality, we can assume that there exists a sufficiently small positive ε such that (16)σ1≜∏n=qsqs+1-1exp⁡r-aPn∗-ε<1.
Next, we continue to prove Theorem2. Obviously, it is easy to show that Pn+1>σPn. Then the following comparison equation with pulse is considered. (17)Qn+1=σQn,n∈Z,Qqk+=1+q2maxQqkθ2+QqkQqk+τ,k∈N.By using Lemma 1, it gives that Pn≥Qn and limn→+∞Qn=Pn∗. Hence, for all n∈Z, the inequality Pn≥Qn>Pn∗-ε holds true. Then according to model (3), we have (18)Hn+1≤Hnexp⁡r-aPn∗-ε,n∈Z,Hqk+=1-q1maxHqkθ1+HqkHqk,k∈N.Associated with the condition (16), the following inequality is derived:(19)Hqs+1≤Hqs+∏n=qsqs+1-1exp⁡r-aPn∗-ε=Hqs1-q1maxHqsθ1+Hqs∏n=qsqs+1-1exp⁡r-aPn∗-ε<Hqs∏n=qsqs+1-1exp⁡r-aPn∗-ε=σ1Hqs,so it is easy to show that Hqs≤H0+σ1s and lims→+∞Hqs=0. It follows from model (3) that 0≤Hn≤Hn-1exp⁡r, and, for any n∈[qs+,q(s+1)), we have (20)0≤Hn≤exp⁡n-qsrHqs+≤1-q1maxHqsθ1+Hqsexp⁡qrHqs,in other words, limn→+∞Hn=0.
Thus, we have proved the global stability of the host-free periodic solution(0,Pn∗) of model (3). This completes the proof.Based on Theorem2, the global stability of the host-free periodic solution of model (3) is determined by R1<1 and R2<1; that is, it is crucial that those two threshold values make great effects on the dynamics of model (3). Therefore, in the following, we employ R1 and R2 to investigate the important factors which affect the threshold values most significantly by numerical simulations.As we know, resource limitation exists extensively in pest control, especially in some developing countries. There are so many factors affecting the limited resource; this study focuses on the thresholdsR1 and R2 which ensure the existence and stability of host-free periodic solution, and those thresholds will directly influence the success of pest control. According to the expressions of R1 and R2, the impulsive periodic q, the maximal release rate q2max, and the half saturation constant θ2 for parasitoid population are the key parameters that affect R1, R2 and the limited resource. In Figure 1(a), with the increasing of the pulse periodic q, the threshold R1 decreases, while the curve of R2-q is a typical single-peak curve. And only if the pulse periodic q>4, max⁡{R1,R2}<1 holds true, then the host-free periodic solution (0,Pn∗) can become stable. The higher the frequency of a control measure is, the lower the pulse periodic q is. Figure 1(a) indicates that the pest can be out of control by increasing the control frequency (i.e., decreasing the control periodic q), and it is not even realistic to do under the limited resource.Figure 1
Sensitivity analysis forR1,R2. (a) The curves of R1-q and R2-q, parameters are r=0.5, a=0.01, q2max=0.95, σ=0.6, θ2=0.2, τ=0.2; (b) the curve of R2-q2max with r=0.4, a=0.05, q2max=0.95, σ=0.6, θ2=4, τ=0.8, q=1; (c) the curve of R2-θ2 with r=0.2, a=0.05, q2max=0.3, σ=0.8, τ=1.2, q=2.
(a) (b) (c)Meanwhile, the optimal release rate of parasitoidq2max is between 0.34 and 0.90 which ensures R2<1, as shown in Figure 1(b). Based on this fact, it is crucial to choose a proper q2max to direct the adaptive resource allocation, especially under the limited control resource. To release the excessive parasitoid may not be good for pest control and will also bring about some other problems such as intra-specific competition for parasitoid population. Analogously, the curve of R2-θ2 shown in Figure 1(c) can also confirm this result, and Figure 1(c) shows that (1) the threshold value R2 is decreased gradually with increasing of the half-saturation constant θ2; (2) the critical value of θ2 is 2.16. From Figures 1(b) and 1(c), we should take appropriate measures to prevent the releasing of parasitoid which will be beneficial to pest control, which is the optimal strategy to wrestle with enormous challenges of the limited resources.This section mainly analyzes the host-free periodic solution; however, there are several other types of solutions for model (3), which will be useful for bifurcation analysis in the next subsection, so we list those solutions as follows.(i)
Host-free periodic solution, denoted by HF: the host population is eradicated, while the parasitoid population oscillates periodically, as shown in Figures2(a) and 2(b), respectively, and the periodic of Parasitoid is 2.(ii)
Host-parasitoid permanence periodic solution, denoted by HP: both individuals can oscillate with a periodic, as shown in Figures2(c) and 2(d), and the periodic of HP is 9.(iii)
Host-outbreak solution, denoted by HO: the host population has periodic or irregular outbreaks with large maximum amplitudes, as shown in Figures2(e) and 2(f).Figure 2
Three types of solutions for model (3). (a-b) HF periodic solution with 9=2; (c-d) HP periodic solution with q=9; (e-f) HO solution with q=20. Parameters are fixed as r=0.3, a=0.1, q1max=0.9, θ1=0.01, q2max=0.4, θ2=0.5, τ=1, and (H0,P0)=(1,0.2).
(a) (b) (c) (d) (e) (f)


### 2.3. Numerical Bifurcation Analysis

In order to get some dynamical behaviors of model (3), a single parameter bifurcation is adopted, which can give information about the dependence of the dynamics on a certain parameter. This subsection focuses on the complex dynamics including period doubling, chaotic solutions, coexistence of multiple attractors, the basin of attraction, and so on.Once the threshold value exceeds a critical level, that is,R2≥1, the host-free solution (0,Pn∗) becomes unstable; both host and parasitoid populations can oscillate periodically with a large amplitude; a series of bifurcation phenomena will be taken place at this moment. To investigate the complex dynamics that model (3) can have, we choose the release rate τ for parasitoids as a bifurcation parameter and fix all other parameters as those in Figure 3 for the different half-saturation constants θ1,θ2. Figures 3(a) and 3(b) are bifurcation diagrams with resource limitation, i.e., θ1=9,θ2=2; Figures 3(c) and 3(d) are bifurcation diagrams without resource limitation, i.e., θ1=θ2=0. Comparing those bifurcation diagrams, we conclude that the nonlinear pulse control measure can produce more complex dynamics than the linear ones, and the limited resource plays a key role.Figure 3
Bifurcation diagram of model (3) with τ. Here (a) & (b) with θ1=9,θ2=2; (c) & (d) with θ1=θ2=0, and all other parameters are a=0.39, r=3.8,σ=0.6,q1max=0.8, q2max=0.3, q=3.
(a) (b) (c) (d)To illustrate the complex and interesting dynamic behaviors of model (3) more clearly, the intrinsic growth rate r and half-saturation constant θ2 are chosen as a bifurcation parameter, respectively; numerical bifurcation analyses are derived to show some possible dynamics as shown in Figures 4 and 5. Specifically, Figure 4 shows that model (3) has periodic doubling bifurcation, periodic window, and chaos solution as the bifurcation parameters r∈(2,71,2.805), (2,819,2.825), and (2,83,2.86), respectively, while model (3) has periodic adding and halving phenomenon as the bifurcation parameters θ2∈(0,0.65) and (0.7,1.5), respectively. Meanwhile, model (3) has some other dynamic behaviors including quasiperiodic solutions, tangent bifurcation, multistability, and chaos crisis in more sensitive parameter space. The results indicate that the parameters of model (3) are highly sensitive, and the dynamical behaviors of model (3) will have a complex change along with parameter variations, even if there is a small perturbation in some key parameters.Figure 4
Bifurcation diagram of model (3) with r. Parameters are fixed as a=0.45, σ=0.4, q1max=0.9, θ1=0.01, q2max=0.1, θ2=0.5, τ=1.2, q=2.
(a) (b)Figure 5
Bifurcation diagram of model (3) with θ2. Parameters are fixed as r=0.8343, a=0.045, σ=0.5, q1max=0.5, θ1=3, q2max=0.5, τ=2, q=3.
(a) (b)From the above bifurcation diagrams (Figures3, 4, and 5), model (3) has several coexisting attractors in a special parameter space. Especially, all parameters are fixed as shown in Figure 6, model (3) has four attractors with different initial densities. Obviously, the cases in Figures 6(a) and 6(b) are a host-free and parasitoid-free solution, and the other cases shown in Figures 6(c)–6(h) are HP periodic solutions, and the maximum amplitudes of host population after the 100-th generation from top to bottom are 6.8,10,2, respectively. Similarly, the maximum amplitudes of parasitoid population after the 100-th generation from top to bottom are 26.5,30.5,23.4, respectively. Figure 6 indicates that those phenomena are caused by different initial densities.Figure 6
Four coexisting attractors of model (3). Parameters are r=2.63, a=0.15, q1max=0.2, σ=0.8, θ1=0.0001, q2max=0.3, θ2=0.00005, τ=10, q=5. The initial densities from top to bottom are (3.01,0.99),(0.93,20),(10.02,25.81),(1,25), respectively.In order to illustrate the initial sensitivities shown in Figure6 more specifically, the basins of attraction with respect to four different coexistence solutions are shown in Figure 7, the parameters are identical to Figure 6, and Figure 7(b) is the magnification part of Figure 7(a). The final stable states of host and parasitoid populations depend on their initial densities, the results indicate that the successful biological control depends on the initial densities of both populations. The proper initial densities can affect host outbreak and help us to design the proper control strategies as well as making management decisions.Figure 7
Basin of several attractions shown in Figure6. The magenta, red, blue, and green points are attracted to the attractor shown in Figure 6 from top to bottom, respectively. The black points represent the other attractors not shown in Figure 6.
(a) (b)


## 2.1. Model Formulation

In 1935, Nicholson and Bailey [33] proposed a classic discrete host-parasitoid model, as an intergenerational survival rate for parasitoids was considered. In 2008, Tang [8] extended the basic Nicholson-Bailey model as follows: (1)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn.Here Hn and Pn represent the density of hosts(or pests) and parasitoids(or natural enemies) at n-th generation, respectively; r is the intrinsic growth rate for the host population; a denotes the searching efficiency of the parasitoid population; the terms exp⁡(-aPn) and [1-exp⁡(-aPn)] represent the probability that the host succeeds and fails in escaping from the parasitism, respectively; σ is the density-independent survival of the parasitoid at n-th generation, n∈Z≜{0,1,2,⋯}.In order to investigate the dynamics of system (1) with IPM strategies, by introducing periodic spraying pesticide only for hosts, and releasing parasitoids, Tang [8] proposed the following Nicholson-Bailey model: (2)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn,n∈ZHqk+=1-q1Hqk,Pqk+=1+q2Pqk+τ,q,k∈Nwith the initial densities (H0+,P0+)=(H0,P0) and N≜{1,2,⋯}. Meanwhile, q∈N is the pulse period, q1 is the killing rate for hosts, and q2 and τ are the release rate and constant for parasitoids at qk-th generation, respectively. Hqk and Pqk (Hqk+ and Pqk+) denote the densities of both populations at qk-th generation before (after) the impulsive perturbations, respectively.Now, taking into account the resource limitation and saturation effects, that is, a linear pulse in model (2) is replaced with a saturation phenomenon of limited resources; especially, a Hill function [34] is used to describe the resource limitation. For convenience and simplification, model (2) can be rewritten as (3)Hn+1=Hnexp⁡r-aPn,Pn+1=Hn1-exp⁡-aPn+σPn,n∈Z,Hqk+=1-q1maxHqkθ1+HqkHqk,Pqk+=1+q2maxPqkθ2+PqkPqk+τ,q,k∈N,where q1max and θ1 represent the maximal fatality rate and the half-saturation constant for hosts, respectively. q2max and θ2 represent the maximal release rate and the half-saturation constant for parasitoids, respectively. Assume that 0≤qimax<1 and θi>0 for i=1,2.


## 2.2. Mathematical Analysis of the Host-Free Periodic Solution and Threshold Conditions

With the main purpose of IPM strategy, it is necessary to illustrate the existence of the host-free periodic solution of model (3) and focus on determining its global attractivity. In order to do this, we first investigate the host-free set: A={(H,P)∈R+2,H=0}, which is clearly invariant by model (3); it follows from model (3) and its pest-free set A that (4)Pn+1=σPn,n∈Z,Pqk+=1+q2maxPqkθ2+PqkPqk+τ,k∈N,P0+=P0.Obviously, model (4) is a periodic system, so the periodic solution Pn can be defined at the subinterval [qs+,q(s+1)) with the initial value Pqs+, s∈Z. Noting that Pqs+ represents the parasitoid’s density after the impulsive perturbation.From the first equation of model (4), we can use this iteration equation at interval [qs+,q(s+1)) and obtain Pn+1=σn+1-qsPqs+. Especially, it gives that Pq(s+1)=σqPqs+ as n=q(s+1)-1.Now, combining the second equation of model (4) and the above one yields (5)Pqs+1+=q2max+1σ2qPqs+2+θ2+τσqPqs++θ2τθ2+σqPqs+.Denoting Pq(s+1)+ as Pn+1, then the above equation can be rewritten as the following difference equation: (6)Pn+1=FPn≜1+q2maxσ2qPn2+θ2+τσqPn+θ2τθ2+σqPn.Equation (6) is a stroboscopic map of model (4): it describes the relations of the parasitoids’ density between any two successive pulse points; that is, the existence of a positive steady state of (6) represents the existence of a positive periodic solution of model (4). Therefore, we first explore the positive steady state of the stroboscopic map (6); let us get the derivative of F(Pn) with Pn and yield (7)F′Pn=q2max+1σ3qPn2+2q2max+1σ2qθ2Pn+θ22σqσ2qPn2+2θ2σqPn+θ22.According to the stability theory of differential equation, it is necessary to ensure that the inequality |F′(Pn)|<1 holds true, and this condition is equivalent to (8)A0Pn2+B0Pn+C0<0,where A0=σ2q[(1+q2max)σq-1], B0=2θ2σq[(1+q2max)σq-1], C0=θ22(σq-1). By a simple calculation, one yields Δ0=B02-4A0C0=4θ22q2maxσqA0. It is easy to find that if A0<0 holds true, then Δ0<0, which implies |F′(Pn)|<1. For simplification, A0<0 can be rewritten as (9)R1≜1+q2maxσq<1.The above inequality can ensure the stability of a positive solution of the stroboscopic map (6). Meanwhile, it also plays an important role in discussing the existence of the positive fixed point for (6).Next, we continue to investigate the positive fixed point of the stroboscopic map (6). For convenience, we denote this fixed point by P~; then it follows from the stroboscopic map (6) that (10)P~=1+q2maxσ2qP~2+θ2+τσqP~+θ2τθ2+σqP~,which means (11)A1P~2+B1P~+C1=0,where A1=σq[1+q2maxσq-1], B1=(θ2+τ)σq-θ2, C1=θ2τ>0. Obviously, we can obtain that A1<0 and Δ1=B12-4A1C1>0 as R1<1 holds true. That is, the stroboscopic map (6) exists in only one positive fixed point P~=-(B1+Δ1)/(2A1). Therefore, we have the following result.Lemma 1.
Model (4) exists as a positive periodic solution Pn∗ provided that R1<1 holds true. Meanwhile, for every positive solution Pn of model (4), it is given that limn→+∞Pn-Pn∗=0. Here Pn∗=σn-qsP~ with n∈[qs+,q(s+1)),s∈Z.Based on Lemma1, we can obtain the general expression of the host-free periodic solution of model (3) in the interval [qs+,q(s+1)) for all s∈Z denoted by (12)0,Pn∗=0,-σn-qsB1+Δ12A1.In order to investigate the global stability of the host-free periodic solution (0,Pn∗) of model (3), the following theorem is derived.Theorem 2.
Let(Hn,Pn) be any positive solution of model (3); the host-free periodic solution (0,Pn∗) of model (3) is globally asymptotically stable in the first quadrant, provided that (13)R2≜2qrA11-σa1-σqB1+Δ1<1.Proof.
To prove Theorem2, the following condition (16) is considered first. It follows from R2<1, model (3), and Pn∗=σn-qsP~ that we have (14)∑n=qsqs+1-1r-aPn∗<0,which implies (15)∏n=qsqs+1-1exp⁡r-aPn∗<1.Therefore, from the above inequality, we can assume that there exists a sufficiently small positive ε such that (16)σ1≜∏n=qsqs+1-1exp⁡r-aPn∗-ε<1.
Next, we continue to prove Theorem2. Obviously, it is easy to show that Pn+1>σPn. Then the following comparison equation with pulse is considered. (17)Qn+1=σQn,n∈Z,Qqk+=1+q2maxQqkθ2+QqkQqk+τ,k∈N.By using Lemma 1, it gives that Pn≥Qn and limn→+∞Qn=Pn∗. Hence, for all n∈Z, the inequality Pn≥Qn>Pn∗-ε holds true. Then according to model (3), we have (18)Hn+1≤Hnexp⁡r-aPn∗-ε,n∈Z,Hqk+=1-q1maxHqkθ1+HqkHqk,k∈N.Associated with the condition (16), the following inequality is derived:(19)Hqs+1≤Hqs+∏n=qsqs+1-1exp⁡r-aPn∗-ε=Hqs1-q1maxHqsθ1+Hqs∏n=qsqs+1-1exp⁡r-aPn∗-ε<Hqs∏n=qsqs+1-1exp⁡r-aPn∗-ε=σ1Hqs,so it is easy to show that Hqs≤H0+σ1s and lims→+∞Hqs=0. It follows from model (3) that 0≤Hn≤Hn-1exp⁡r, and, for any n∈[qs+,q(s+1)), we have (20)0≤Hn≤exp⁡n-qsrHqs+≤1-q1maxHqsθ1+Hqsexp⁡qrHqs,in other words, limn→+∞Hn=0.
Thus, we have proved the global stability of the host-free periodic solution(0,Pn∗) of model (3). This completes the proof.Based on Theorem2, the global stability of the host-free periodic solution of model (3) is determined by R1<1 and R2<1; that is, it is crucial that those two threshold values make great effects on the dynamics of model (3). Therefore, in the following, we employ R1 and R2 to investigate the important factors which affect the threshold values most significantly by numerical simulations.As we know, resource limitation exists extensively in pest control, especially in some developing countries. There are so many factors affecting the limited resource; this study focuses on the thresholdsR1 and R2 which ensure the existence and stability of host-free periodic solution, and those thresholds will directly influence the success of pest control. According to the expressions of R1 and R2, the impulsive periodic q, the maximal release rate q2max, and the half saturation constant θ2 for parasitoid population are the key parameters that affect R1, R2 and the limited resource. In Figure 1(a), with the increasing of the pulse periodic q, the threshold R1 decreases, while the curve of R2-q is a typical single-peak curve. And only if the pulse periodic q>4, max⁡{R1,R2}<1 holds true, then the host-free periodic solution (0,Pn∗) can become stable. The higher the frequency of a control measure is, the lower the pulse periodic q is. Figure 1(a) indicates that the pest can be out of control by increasing the control frequency (i.e., decreasing the control periodic q), and it is not even realistic to do under the limited resource.Figure 1
Sensitivity analysis forR1,R2. (a) The curves of R1-q and R2-q, parameters are r=0.5, a=0.01, q2max=0.95, σ=0.6, θ2=0.2, τ=0.2; (b) the curve of R2-q2max with r=0.4, a=0.05, q2max=0.95, σ=0.6, θ2=4, τ=0.8, q=1; (c) the curve of R2-θ2 with r=0.2, a=0.05, q2max=0.3, σ=0.8, τ=1.2, q=2.
(a) (b) (c)Meanwhile, the optimal release rate of parasitoidq2max is between 0.34 and 0.90 which ensures R2<1, as shown in Figure 1(b). Based on this fact, it is crucial to choose a proper q2max to direct the adaptive resource allocation, especially under the limited control resource. To release the excessive parasitoid may not be good for pest control and will also bring about some other problems such as intra-specific competition for parasitoid population. Analogously, the curve of R2-θ2 shown in Figure 1(c) can also confirm this result, and Figure 1(c) shows that (1) the threshold value R2 is decreased gradually with increasing of the half-saturation constant θ2; (2) the critical value of θ2 is 2.16. From Figures 1(b) and 1(c), we should take appropriate measures to prevent the releasing of parasitoid which will be beneficial to pest control, which is the optimal strategy to wrestle with enormous challenges of the limited resources.This section mainly analyzes the host-free periodic solution; however, there are several other types of solutions for model (3), which will be useful for bifurcation analysis in the next subsection, so we list those solutions as follows.(i)
Host-free periodic solution, denoted by HF: the host population is eradicated, while the parasitoid population oscillates periodically, as shown in Figures2(a) and 2(b), respectively, and the periodic of Parasitoid is 2.(ii)
Host-parasitoid permanence periodic solution, denoted by HP: both individuals can oscillate with a periodic, as shown in Figures2(c) and 2(d), and the periodic of HP is 9.(iii)
Host-outbreak solution, denoted by HO: the host population has periodic or irregular outbreaks with large maximum amplitudes, as shown in Figures2(e) and 2(f).Figure 2
Three types of solutions for model (3). (a-b) HF periodic solution with 9=2; (c-d) HP periodic solution with q=9; (e-f) HO solution with q=20. Parameters are fixed as r=0.3, a=0.1, q1max=0.9, θ1=0.01, q2max=0.4, θ2=0.5, τ=1, and (H0,P0)=(1,0.2).
(a) (b) (c) (d) (e) (f)


## 2.3. Numerical Bifurcation Analysis

In order to get some dynamical behaviors of model (3), a single parameter bifurcation is adopted, which can give information about the dependence of the dynamics on a certain parameter. This subsection focuses on the complex dynamics including period doubling, chaotic solutions, coexistence of multiple attractors, the basin of attraction, and so on.Once the threshold value exceeds a critical level, that is,R2≥1, the host-free solution (0,Pn∗) becomes unstable; both host and parasitoid populations can oscillate periodically with a large amplitude; a series of bifurcation phenomena will be taken place at this moment. To investigate the complex dynamics that model (3) can have, we choose the release rate τ for parasitoids as a bifurcation parameter and fix all other parameters as those in Figure 3 for the different half-saturation constants θ1,θ2. Figures 3(a) and 3(b) are bifurcation diagrams with resource limitation, i.e., θ1=9,θ2=2; Figures 3(c) and 3(d) are bifurcation diagrams without resource limitation, i.e., θ1=θ2=0. Comparing those bifurcation diagrams, we conclude that the nonlinear pulse control measure can produce more complex dynamics than the linear ones, and the limited resource plays a key role.Figure 3
Bifurcation diagram of model (3) with τ. Here (a) & (b) with θ1=9,θ2=2; (c) & (d) with θ1=θ2=0, and all other parameters are a=0.39, r=3.8,σ=0.6,q1max=0.8, q2max=0.3, q=3.
(a) (b) (c) (d)To illustrate the complex and interesting dynamic behaviors of model (3) more clearly, the intrinsic growth rate r and half-saturation constant θ2 are chosen as a bifurcation parameter, respectively; numerical bifurcation analyses are derived to show some possible dynamics as shown in Figures 4 and 5. Specifically, Figure 4 shows that model (3) has periodic doubling bifurcation, periodic window, and chaos solution as the bifurcation parameters r∈(2,71,2.805), (2,819,2.825), and (2,83,2.86), respectively, while model (3) has periodic adding and halving phenomenon as the bifurcation parameters θ2∈(0,0.65) and (0.7,1.5), respectively. Meanwhile, model (3) has some other dynamic behaviors including quasiperiodic solutions, tangent bifurcation, multistability, and chaos crisis in more sensitive parameter space. The results indicate that the parameters of model (3) are highly sensitive, and the dynamical behaviors of model (3) will have a complex change along with parameter variations, even if there is a small perturbation in some key parameters.Figure 4
Bifurcation diagram of model (3) with r. Parameters are fixed as a=0.45, σ=0.4, q1max=0.9, θ1=0.01, q2max=0.1, θ2=0.5, τ=1.2, q=2.
(a) (b)Figure 5
Bifurcation diagram of model (3) with θ2. Parameters are fixed as r=0.8343, a=0.045, σ=0.5, q1max=0.5, θ1=3, q2max=0.5, τ=2, q=3.
(a) (b)From the above bifurcation diagrams (Figures3, 4, and 5), model (3) has several coexisting attractors in a special parameter space. Especially, all parameters are fixed as shown in Figure 6, model (3) has four attractors with different initial densities. Obviously, the cases in Figures 6(a) and 6(b) are a host-free and parasitoid-free solution, and the other cases shown in Figures 6(c)–6(h) are HP periodic solutions, and the maximum amplitudes of host population after the 100-th generation from top to bottom are 6.8,10,2, respectively. Similarly, the maximum amplitudes of parasitoid population after the 100-th generation from top to bottom are 26.5,30.5,23.4, respectively. Figure 6 indicates that those phenomena are caused by different initial densities.Figure 6
Four coexisting attractors of model (3). Parameters are r=2.63, a=0.15, q1max=0.2, σ=0.8, θ1=0.0001, q2max=0.3, θ2=0.00005, τ=10, q=5. The initial densities from top to bottom are (3.01,0.99),(0.93,20),(10.02,25.81),(1,25), respectively.In order to illustrate the initial sensitivities shown in Figure6 more specifically, the basins of attraction with respect to four different coexistence solutions are shown in Figure 7, the parameters are identical to Figure 6, and Figure 7(b) is the magnification part of Figure 7(a). The final stable states of host and parasitoid populations depend on their initial densities, the results indicate that the successful biological control depends on the initial densities of both populations. The proper initial densities can affect host outbreak and help us to design the proper control strategies as well as making management decisions.Figure 7
Basin of several attractions shown in Figure6. The magenta, red, blue, and green points are attracted to the attractor shown in Figure 6 from top to bottom, respectively. The black points represent the other attractors not shown in Figure 6.
(a) (b)


## 3. State-Dependent Impulsive Control for a Nicholson-Bailey Model with Resource Limitation

### 3.1. Model Formulation

In Section2, model (3) with a fixed-time impulsive control strategy has been investigated. One of the greatest advantages of this periodic control strategy is the purpose of eradicating pests rapidly and implementing easily, so it is widely applied and extended in agriculture and ecology. However, the most fatal disadvantage is that it can cause serious environmental pollution and resource waste, which defeats the purpose of IPM strategy. Especially in the resource limitation, it is always questioned, so it is necessary to look for a new measure to control pest.It is well known that the main purpose of IPM is to maintain the density of the pests below the EIL rather than seeking to eradicate them, and the suitable tactic will be applied only when the density of host reaches the given ET. That is the threshold control strategy. Therefore, in model (1), only when the density of pests exceedes the given threshold, the pesticide will be applied to host population and decrease its density to ET, that is (21)1-q1nHnexp⁡r-aPn=ET,where q1n is the killing rate for host, and the survival rate of host after a pesticide application is 1-q1n. With the idea of the threshold control strategy, we can redefine the killing rate as follows: (22)q1n=1-ETHnexp⁡r-aPn,ifHn+1>ET,0,ifHn+1≤ET.By employing threshold policy control and IPM strategy, we can use the combination of the biological (releasing the parasitoid) and chemical (spraying pesticide) tactics to suppress the pest to a controllable level ET. Meanwhile, the resource limitations and saturation effects are taken into account; for convenience and simplification, then model (3) can be rewritten as (23)Hn+1=min⁡Hnexp⁡r-aPn,ET,Pn+1=Hn1-exp⁡-aPn+σPn,Hn+1≠ET,n∈ZPn+1+=1+q2maxPnθ2+PnPn+1+τ,Hn+1=ETwith (H0+,P0+)=(H0,P0) and H0<ET. The term q2maxPn/(θ2+Pn) represents the half-saturation function, and the parameters are the same as model (3).


### 3.2. Numerical Bifurcation Analysis

Compared with the fixed-time impulsive model (3), the state-dependent impulsive model (23) describes the threshold control policy and IPM strategy more complexly, and it is difficult to solve model (23) and especially in resource limitation model (23) becomes a complicated nonlinear dynamic system; it is difficult to investigate model (23) by theoretically, so we will investigate the dynamical behavior of model (23) by numerical simulations in the subsection and the analysis focus on the effects of some key parameters on dynamics of model (23) with limited resource.Firstly, we choose the intrinsic growth rater of hosts as a bifurcation parameter; the bifurcation diagrams of model (23) with the killing rate q1n are derived. As shown in Figure 8, both populations Hn,Pn and the killing rate q1n have quite interesting and complicated dynamics, including periodic doubling bifurcation, periodic window, chaotic solutions, and so on, when the intrinsic growth rate r varies from 1 to 2.8. The killing rate q1n always oscillates between 0 and 1, which can pose a major challenge to be estimated and predicted.Figure 8
Bifurcation diagram of model (23) with r. Parameters are fixed as a=0.25, σ=0.6, q2max=0.1, θ2=1, τ=2, ET=12.
(a) (b) (c)Meanwhile, there are several types of solution including HF, HP periodic solutions, and HO solution; i.e., model (23) has several attractors with different initial densities. Figure 8 shows that model (23) has at least four different HP solutions or attractors. The maximum amplitudes of parasitoid population from top to bottom are 21.316,23.996,19.2084,26.045, respectively; see Figures 9(a)–9(d) for details.Figure 9
Four coexisting attractors of model (23). Parameters are r=3, a=0.2, σ=0.7, q2max=0.5, θ2=0.05, τ=8,ET=2. The initial values from top to bottom are (1,1),(0.2,4),(1,9),(0.2,20), respectively, and q1(0)=0.Furthermore, in order to illustrate the different initial densities that can cause the different dynamics, the basins of attractors are given in Figure10. It shows that the dynamical behavior of model (23) and the host-killing rate q1n are quite sensitive with varies of their initial value (H0,P0,q10). Therefore, to successfully carry out a control pest plan, it is necessary to fully understand the initial state and analyze the changes in dynamics.Figure 10
Basin of several attractions shown in Figure9. The red, magenta, blue, and green points are attracted to the attractor shown in Figure 9 from top to bottom, respectively. The yellow points represent the other attractors not shown in Figure 9. All the parameters are identical to Figure 9. (b) is the magnification part of (a).
(a) (b)


## 3.1. Model Formulation

In Section2, model (3) with a fixed-time impulsive control strategy has been investigated. One of the greatest advantages of this periodic control strategy is the purpose of eradicating pests rapidly and implementing easily, so it is widely applied and extended in agriculture and ecology. However, the most fatal disadvantage is that it can cause serious environmental pollution and resource waste, which defeats the purpose of IPM strategy. Especially in the resource limitation, it is always questioned, so it is necessary to look for a new measure to control pest.It is well known that the main purpose of IPM is to maintain the density of the pests below the EIL rather than seeking to eradicate them, and the suitable tactic will be applied only when the density of host reaches the given ET. That is the threshold control strategy. Therefore, in model (1), only when the density of pests exceedes the given threshold, the pesticide will be applied to host population and decrease its density to ET, that is (21)1-q1nHnexp⁡r-aPn=ET,where q1n is the killing rate for host, and the survival rate of host after a pesticide application is 1-q1n. With the idea of the threshold control strategy, we can redefine the killing rate as follows: (22)q1n=1-ETHnexp⁡r-aPn,ifHn+1>ET,0,ifHn+1≤ET.By employing threshold policy control and IPM strategy, we can use the combination of the biological (releasing the parasitoid) and chemical (spraying pesticide) tactics to suppress the pest to a controllable level ET. Meanwhile, the resource limitations and saturation effects are taken into account; for convenience and simplification, then model (3) can be rewritten as (23)Hn+1=min⁡Hnexp⁡r-aPn,ET,Pn+1=Hn1-exp⁡-aPn+σPn,Hn+1≠ET,n∈ZPn+1+=1+q2maxPnθ2+PnPn+1+τ,Hn+1=ETwith (H0+,P0+)=(H0,P0) and H0<ET. The term q2maxPn/(θ2+Pn) represents the half-saturation function, and the parameters are the same as model (3).


## 3.2. Numerical Bifurcation Analysis

Compared with the fixed-time impulsive model (3), the state-dependent impulsive model (23) describes the threshold control policy and IPM strategy more complexly, and it is difficult to solve model (23) and especially in resource limitation model (23) becomes a complicated nonlinear dynamic system; it is difficult to investigate model (23) by theoretically, so we will investigate the dynamical behavior of model (23) by numerical simulations in the subsection and the analysis focus on the effects of some key parameters on dynamics of model (23) with limited resource.Firstly, we choose the intrinsic growth rater of hosts as a bifurcation parameter; the bifurcation diagrams of model (23) with the killing rate q1n are derived. As shown in Figure 8, both populations Hn,Pn and the killing rate q1n have quite interesting and complicated dynamics, including periodic doubling bifurcation, periodic window, chaotic solutions, and so on, when the intrinsic growth rate r varies from 1 to 2.8. The killing rate q1n always oscillates between 0 and 1, which can pose a major challenge to be estimated and predicted.Figure 8
Bifurcation diagram of model (23) with r. Parameters are fixed as a=0.25, σ=0.6, q2max=0.1, θ2=1, τ=2, ET=12.
(a) (b) (c)Meanwhile, there are several types of solution including HF, HP periodic solutions, and HO solution; i.e., model (23) has several attractors with different initial densities. Figure 8 shows that model (23) has at least four different HP solutions or attractors. The maximum amplitudes of parasitoid population from top to bottom are 21.316,23.996,19.2084,26.045, respectively; see Figures 9(a)–9(d) for details.Figure 9
Four coexisting attractors of model (23). Parameters are r=3, a=0.2, σ=0.7, q2max=0.5, θ2=0.05, τ=8,ET=2. The initial values from top to bottom are (1,1),(0.2,4),(1,9),(0.2,20), respectively, and q1(0)=0.Furthermore, in order to illustrate the different initial densities that can cause the different dynamics, the basins of attractors are given in Figure10. It shows that the dynamical behavior of model (23) and the host-killing rate q1n are quite sensitive with varies of their initial value (H0,P0,q10). Therefore, to successfully carry out a control pest plan, it is necessary to fully understand the initial state and analyze the changes in dynamics.Figure 10
Basin of several attractions shown in Figure9. The red, magenta, blue, and green points are attracted to the attractor shown in Figure 9 from top to bottom, respectively. The yellow points represent the other attractors not shown in Figure 9. All the parameters are identical to Figure 9. (b) is the magnification part of (a).
(a) (b)


## 4. Discussion and Conclusion

In order to investigate the effect of resource limitation on pest control, this paper develops two Nicholso-Bailey IPM models with nonlinear impulsive control strategies. Theoretical, numerical, and biological analyses are given; those results show that it is critical to release parasitoid with an optimal number as the excessive releasing may cause intra-specific competition of parasitoid, which will seriously affect the pest control, especially for the limited resource. Meanwhile, the initial densities of both host and parasitoid populations can affect the pest outbreaks.Model (3) describes the periodic control strategy, that is, a fixed-time impulsive control for a host-parasitoid model with resource limitation. According to the theoretical analysis, the existence and stability of host-free periodic solution of model (3) are derived; numerical threshold results show that we should take appropriate measures to prevent the releasing of parasitoid which will be beneficial to pest control, which is the optimal strategy to wrestle with enormous challenges of the limited resources. Meanwhile, numerical bifurcation results indicate that the nonlinear impulsive can make the dynamics of model (3) become more complex and interesting. There are several different bifurcation phenomena; we can find that the routes to chaos are very complicated; that is, there are several hidden factors that can adversely affect our control strategy. Those present a major challenge for controlling the pest populations in practice.Model (23) describes the threshold control strategy in which measures could overcome the cost ineffectiveness or potentially damage to the environment in model (3). The threshold control strategy is to maintain the density of the hosts below the EIL rather than seeking to eradicate them, and the suitable tactic will be applied only as the density of hosts reaches the given ET. Numerical simulations clarify that the initial densities of both populations can affect the outcome of classical biological control, and the final stable states of host and parasitoid populations depend on their initial densities. Those results are further confirmed by basins of attraction of initial densities.This work only focuses on the limited resource; there are several other factors that can influence the pest control in reality, including Allee effect, delayed responses, and residual effects on hosts of pesticides. If those factors are considered, it will be of great theoretical and practical significance to investigate the effective prevention and control strategy for agricultural pests. Addressing these issues needs more future work.



---

*Source: 9365293-2019-01-14.xml*