# Numerical Analysis of Discrete Switching Prey-Predator Model for Integrated Pest Management

Authors: Changcheng Xiang; Yi Yang; Zhongyi Xiang; Wenjie Qin
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

The switching discrete prey-predator model concerning integrated pest management has been proposed, and the switches are guided by the economic threshold (ET). To begin with, the regular and virtual equilibria of switching system have been discussed and the key parameter bifurcation diagrams for the existence of equilibria have been proposed, which reveal the three different regions of equilibria. Besides, numerical bifurcation analyses show that the switching discrete system may have complicated dynamics behavior including chaos and the coexistence of multiple attractors. Finally, the effects of key parameters on the switching frequencies and switching times are discussed and the sensitivity analysis of varying parameter values for mean switching times has also been given. The results proved that economic threshold (ET) and the growth rate (α) were the key parameters for pest control.


---

## Body

## 1. Introduction

The integrated pest management (IPM) is applied [1–3]. IPM is a long-term control strategy that combines biological, cultural, and chemical tactics to reduce the density of pest populations to economic injury level (EIL) when the pests reach an economic threshold (ET) [1, 2]. The classic Lotka-Volterra type prey-predator system for pest control presents as(1)dxdt=rx1-xK-bxy,dydt=cxy-dy,where x(t) and y(t) denote prey and predator densities, respectively; all the parameters in system (1) are positive constants.For the sake of simplicity, we rewrite system (1) as follows:(2)dxdt=αx1-x-βxy,dydt=βx-δy,where α=rK,β=K2c, and δ=dK. The discrete-time prey-predator system can be obtained by the forward Euler scheme to system (2):(3)xt+1=xt+σαxt1-xt-βxtyt,yt+1=yt+σ-δ+βxtyt,where σ is the step size. In model (3), if the prey stands for the pest population and the predator stands for the natural enemies, then what we want to know is that how many key parameters affect the pest control when integrated pest management (IPM) is applied. An ET is usually defined as the density of pest population in the field when control actions must be taken to prevent the EIL from being reached and exceeded. This indicates that the main objective of IPM is to maintain the density of pest population below the EIL rather than eradicate it.In this study, The prey population and the predator population are regarded as the pest and the natural enemy, respectively. we want to discuss discrete-time prey-predator models by utilizing a threshold policy (TP) to control the pest population (prey). When the pest population density is above the economic threshold (ET), the pesticides are applied for pest population and the natural enemies (predator) population is released simultaneously. The IPM subsystem withx(t)≥ET based on (3) is as follows:(4)xt+1=1-pxt+σαxt1-xt-βxtyt,yt+1=yt+σ-δ+βxtyt+στyt,xt≥ET,where 0<p<1 is the reduction proportion of the pest population density by killing or trapping when the number of pests reaches ET and τ is the increasing proportion of natural enemies when pest population exceeds ET. In order to control the pest, in this paper we assume τ>0. Furthermore, if p=0, we only release sufficient natural enemies to prevent the pest population from exceeding ET. If p>0 and τ=0, then we adopt the chemical control strategy so that it can ensure pest population decrease when the pest population reaches the ET.ByR+2={(x(t),y(t))∣x(t)≥0,y(t)≥0}, the nonnegative octant could be split in two parts: (5)G1=xt,yt∈R+2∣xt<ET,G2=xt,yt∈R+2∣xt≥ET.In the regionG1 the more profitable prey density is below ET in system (3). In the region G2, if the prey density is above ET in system (4), the IPM strategy is applied, which includes spraying the pesticides and releasing the natural enemies.We combine system (3) in region G1 and system (4) in region G2(6)xt+1=xt+σαxt1-xt-βxtyt,yt+1=yt+σ-δ+βxtyt,xt<ET,xt+1=1-σpxt+σαxt1-xt-βxtyt,yt+1=yt+σ-δ+βxtyt+στyt,xt≥ET,so systems (6) in region G1 and in region G2 are divided into two different subsystems by utilizing threshold policy (ET), which are denoted as SG1 and SG2, respectively.


## 2. Fixed Point in Discrete Switching System

In this section, from the point of pest control, we consider the discrete-time switching system (6) in the closed first quadrant R+2 of the (x,y) plane. We first discuss the existence of fixed points for subsystems SG1 and SG2, respectively, and then study the switching behavior of fixed points by utilizing a threshold policy.


### 2.1. Existence and Locally Asymptotic Stability of Fixed Points for the SubsystemsSG1 and SG2

For the subsystemSG1, we are interested in the positive point; let x(t+1)=x(t)=x∗,y(t+1)=y(t)=y∗. By a simple computation, it is straightforward to obtain the following results.Ifβ>δ, the subsystem SG1 has a unique positive fixed point (x1∗,y1∗)=δ/β,αβ-δ/β2. The stability of the fixed point has been discussed by Liu and Xiao [4].We have similar analysis for subsystemSG2. It is clear that subsystem SG2 has unique positive fixed point (x2∗,y2∗)=(δ-τ)/β,(α(β-δ)+ατ-βp)/β2, if δ>τ and α(β-δ)>βp-ατ. Before we investigate the stability of the fixed point, we need the following lemma, which can be easily proved by the relations between roots and coefficients of a quadratic equation [5].Lemma 1.
LetF(λ)=λ2+Pλ+Q. Suppose that F(1)>0 and λ1,λ2 are two roots of F(λ)=0. Then(i)
λ1<1, and λ2<1 if and only if F(-1)>0 and Q>1;(ii)
λ1<1, and λ2>1 (or λ1>1, and λ2<1) if and only if F(-1)<0;(iii)
λ1>1, and λ2>1 if and only if F(-1)>0 and Q>1;(iv)
λ1=-1, and λ2≠1 if and only if F(-1)=0 and P≠0,2;(v)
λ1 and λ2 are complex and λ1=1 and λ2=1 if and only if P2-4Q<0 and Q=1.Note that the local stability of the fixed point(x2∗,y2∗) is determined by the modules of eigenvalues of the characteristic equation. The Jacobian matrix J of the subsystem SG2 at (x2∗,y2∗) is given by(7)Jx2∗,y2∗=-2ασx2∗+ασ-pσ+1-σβy2∗-σβx2∗σβy2∗1-σδ+σβx2∗+στand the characteristic equation of the Jacobian matrix J(x2∗,y2∗) can be written as(8)λ2+-2+Hσλ+1-Hσ+Gσ2=0.The characteristic polynomial is denoted as(9)Fλ=λ2+-2+Hσλ+1-Hσ+Gσ2,where(10)G=-αδ2β+αδ-ατ+pτ+2αδτβ-ατ2β-δp=δ-ταβ-αδ+ατ-pββ>0,H=αδβ-ατβ=αδ-τβ>0.ThenF(1)=Gσ2>0 and F(-1)=4-2Hσ2+Gσ2.Letλ1 and λ2 are two roots of the characteristic equation (8), which are called eigenvalues of the fixed point (x2∗,y2∗). The fixed point (x2∗,y2∗) is called a sink if λ1<1 and λ2<1; the sink is locally asymptotically stable. (x2∗,y2∗) is called a source if λ1>1 and λ2>1; the source is locally unstable. (x2∗,y2∗) is called a saddle if λ1<1 and λ2>1 (or λ1>1 and λ2<1). (x2∗,y2∗) is called nonhyperbolic if either λ1=1 or λ2=1.Then we can calculate the local dynamics of the fixed point(x2∗,y2∗) by Lemma 1.Proposition 2.
Let(x2∗,y2∗) be the positive fixed point of (8).(i)
It is a sink if one of the following conditions holds:
(i.1)
G≥H2/4 and σ>(H-H2-4G)/G;(i.2)
G<H2/4 and H/G>σ>0.(ii)
It is a source if one of the following conditions holds:
(ii.1)
G≥H2/4 and σ>(H+H2-4G)/G;(ii.2)
G<H2/4 and σ>H/G.(iii)
It is a saddle if the following condition holds:G>H2/4 and (H+H2-4G)/G>σ>(H-H2-4G)/G.(iv)
It is nonhyperbolic if the following condition holds:G>H2/4 and σ=(H±H2-4G)/G and σ≠2/H,4/H.


### 2.2. Equilibria for Switching Dynamics of System (6)

In the following section, the definition of the real and virtual equilibria concepts is briefly summarized [6–8], and several types of equilibrium are discussed.Definition 3.
LetZi=(xi∗,yi∗) be the stable equilibrium point of the dynamics of region Gi for i=1,2. Then Zi is called a real equilibrium if it belongs to subsystem SGi and a virtual equilibrium if it belongs to SGj, j≠i, for i,j=1,2. The real equilibrium is denoted as Zr1 and Zr2 and the virtual equilibrium points are denoted as Zv1 and Zv2, respectively.For the subsystemSG1 with x<ET (i.e., it belongs to region G1), if β>δ and δ/β<ET, then it is a real equilibrium for the system SG1, denoted by Zr1. If β>δ and δ/β≥ET, then it is a virtual equilibrium for the system SG1, denoted by Zv1. About real equilibria for the subsystem SG2 with x≥ET (i.e., it belongs to region G2), on the basis of coordinate of equilibrium Z2=(x2∗,y2∗), we have the following results. If δ>τ,α(β-δ)>βp-ατ, and (δ-τ)/β≥ET, then it is a real equilibrium for the system SG2, denoted by Zr2. If δ>τ,α(β-δ)>βp-ατ, and (δ-τ)/β<ET, then it is a virtual equilibrium for the system SG2, denoted by Zv2.In order to consider the relationship of equilibria for two subsystems, we define five curves as follows:L1=δ/β,L2=(δ-τ)/β,L3=αβ-δ/β2,L4=(α(β-δ)+ατ-βp)/β2, and L5=ET. Based on the above analysis and (δ-τ)/β<δ/β, we consider the following three cases.(i)
IfET≤(δ-τ)/β<δ/β, the equilibria Zv1 and Zr2 can coexist as shown in Figure 1(a).(ii)
If(δ-τ)/β<ET≤δ/β, the equilibria Zv1 and Zv2 can coexist (i.e., Figure 1(b)).(iii)
If(δ-τ)/β<δ/β<ET, then the equilibria Zr1 and Zv2 can coexist (shown in Figures 1(c) and 1(d)). Moreover, when ET=0.8, the attractor converges to the Zr1 (see Figure 1(d)).Figure 1
Real/virtual equilibria. The parameter values are fixed as follows:p=0.1,τ=1,σ=0.35,β=4,α=3, and δ=2.0. (a) Zv1 and Zr2 coexist with ET=0.2. (b) Zv1 and Zv2 coexist with ET=0.45. (c) Zr1 and Zv2 coexist with ET=0.52. (d) Zr1 and Zv2 coexist with ET=0.8.
(a) (b) (c) (d)In the process of pest control, we should use IPM strategies to prevent pest outbreaks (i.e.,x(t)<ET). According to the above discussions, it is obvious that β,τ, and ET are key factors in determining the existence of the above different types of equilibria of system (6). So we discuss the β and ET parameter space and β and τ parameter space, respectively, which can be divided into three regions (see Figure 2), and coexistence of real or virtual equilibria is indicated in each region. Figure 2(a) shows the β and ET parameter space with δ=2.0 and τ=1 and Figure 2(b) shows parameter space with ET=0.8 and δ=2. In region I, Zr1 and Zv2 coexist, and Zv1 and Zr2 can coexist in region III. Moreover, It is obvious that the two virtual equilibria Zv1 and Zv2 can coexist in region II. The region II is very important for pest control. In practice, we should apply IPM strategies to prevent multiple pest outbreaks (i.e., x(t)≥ET) in the process of pest control. So, we should choose a set of parameters such that all the equilibria of subsystems SG1 and SG2 are virtual equilibria, which have been widely used in pest control and plant disease [1, 3, 9]. Therefore, from the perspective of pest management, the appropriate control strategies and ET should be designed such that the interior equilibrium of SG1 and SG2 is virtual simultaneously; the region II (purple) is the ideal region to design optimal control strategy as shown in Figure 2.Figure 2
Parameter bifurcation diagram for the existence of equilibria of system (6). Region I (red): Zr1 and Zv2 coexist. Region II (purple): Zv1 and Zv2 coexist. Region III (green): Zv1 and Zr2 coexist. (a) β and ET parameter space with δ=2.0 and τ=1. (b) β and τ parameter space with ET=0.8 and δ=2.
(a) (b)


## 2.1. Existence and Locally Asymptotic Stability of Fixed Points for the SubsystemsSG1 and SG2

For the subsystemSG1, we are interested in the positive point; let x(t+1)=x(t)=x∗,y(t+1)=y(t)=y∗. By a simple computation, it is straightforward to obtain the following results.Ifβ>δ, the subsystem SG1 has a unique positive fixed point (x1∗,y1∗)=δ/β,αβ-δ/β2. The stability of the fixed point has been discussed by Liu and Xiao [4].We have similar analysis for subsystemSG2. It is clear that subsystem SG2 has unique positive fixed point (x2∗,y2∗)=(δ-τ)/β,(α(β-δ)+ατ-βp)/β2, if δ>τ and α(β-δ)>βp-ατ. Before we investigate the stability of the fixed point, we need the following lemma, which can be easily proved by the relations between roots and coefficients of a quadratic equation [5].Lemma 1.
LetF(λ)=λ2+Pλ+Q. Suppose that F(1)>0 and λ1,λ2 are two roots of F(λ)=0. Then(i)
λ1<1, and λ2<1 if and only if F(-1)>0 and Q>1;(ii)
λ1<1, and λ2>1 (or λ1>1, and λ2<1) if and only if F(-1)<0;(iii)
λ1>1, and λ2>1 if and only if F(-1)>0 and Q>1;(iv)
λ1=-1, and λ2≠1 if and only if F(-1)=0 and P≠0,2;(v)
λ1 and λ2 are complex and λ1=1 and λ2=1 if and only if P2-4Q<0 and Q=1.Note that the local stability of the fixed point(x2∗,y2∗) is determined by the modules of eigenvalues of the characteristic equation. The Jacobian matrix J of the subsystem SG2 at (x2∗,y2∗) is given by(7)Jx2∗,y2∗=-2ασx2∗+ασ-pσ+1-σβy2∗-σβx2∗σβy2∗1-σδ+σβx2∗+στand the characteristic equation of the Jacobian matrix J(x2∗,y2∗) can be written as(8)λ2+-2+Hσλ+1-Hσ+Gσ2=0.The characteristic polynomial is denoted as(9)Fλ=λ2+-2+Hσλ+1-Hσ+Gσ2,where(10)G=-αδ2β+αδ-ατ+pτ+2αδτβ-ατ2β-δp=δ-ταβ-αδ+ατ-pββ>0,H=αδβ-ατβ=αδ-τβ>0.ThenF(1)=Gσ2>0 and F(-1)=4-2Hσ2+Gσ2.Letλ1 and λ2 are two roots of the characteristic equation (8), which are called eigenvalues of the fixed point (x2∗,y2∗). The fixed point (x2∗,y2∗) is called a sink if λ1<1 and λ2<1; the sink is locally asymptotically stable. (x2∗,y2∗) is called a source if λ1>1 and λ2>1; the source is locally unstable. (x2∗,y2∗) is called a saddle if λ1<1 and λ2>1 (or λ1>1 and λ2<1). (x2∗,y2∗) is called nonhyperbolic if either λ1=1 or λ2=1.Then we can calculate the local dynamics of the fixed point(x2∗,y2∗) by Lemma 1.Proposition 2.
Let(x2∗,y2∗) be the positive fixed point of (8).(i)
It is a sink if one of the following conditions holds:
(i.1)
G≥H2/4 and σ>(H-H2-4G)/G;(i.2)
G<H2/4 and H/G>σ>0.(ii)
It is a source if one of the following conditions holds:
(ii.1)
G≥H2/4 and σ>(H+H2-4G)/G;(ii.2)
G<H2/4 and σ>H/G.(iii)
It is a saddle if the following condition holds:G>H2/4 and (H+H2-4G)/G>σ>(H-H2-4G)/G.(iv)
It is nonhyperbolic if the following condition holds:G>H2/4 and σ=(H±H2-4G)/G and σ≠2/H,4/H.


## 2.2. Equilibria for Switching Dynamics of System (6)

In the following section, the definition of the real and virtual equilibria concepts is briefly summarized [6–8], and several types of equilibrium are discussed.Definition 3.
LetZi=(xi∗,yi∗) be the stable equilibrium point of the dynamics of region Gi for i=1,2. Then Zi is called a real equilibrium if it belongs to subsystem SGi and a virtual equilibrium if it belongs to SGj, j≠i, for i,j=1,2. The real equilibrium is denoted as Zr1 and Zr2 and the virtual equilibrium points are denoted as Zv1 and Zv2, respectively.For the subsystemSG1 with x<ET (i.e., it belongs to region G1), if β>δ and δ/β<ET, then it is a real equilibrium for the system SG1, denoted by Zr1. If β>δ and δ/β≥ET, then it is a virtual equilibrium for the system SG1, denoted by Zv1. About real equilibria for the subsystem SG2 with x≥ET (i.e., it belongs to region G2), on the basis of coordinate of equilibrium Z2=(x2∗,y2∗), we have the following results. If δ>τ,α(β-δ)>βp-ατ, and (δ-τ)/β≥ET, then it is a real equilibrium for the system SG2, denoted by Zr2. If δ>τ,α(β-δ)>βp-ατ, and (δ-τ)/β<ET, then it is a virtual equilibrium for the system SG2, denoted by Zv2.In order to consider the relationship of equilibria for two subsystems, we define five curves as follows:L1=δ/β,L2=(δ-τ)/β,L3=αβ-δ/β2,L4=(α(β-δ)+ατ-βp)/β2, and L5=ET. Based on the above analysis and (δ-τ)/β<δ/β, we consider the following three cases.(i)
IfET≤(δ-τ)/β<δ/β, the equilibria Zv1 and Zr2 can coexist as shown in Figure 1(a).(ii)
If(δ-τ)/β<ET≤δ/β, the equilibria Zv1 and Zv2 can coexist (i.e., Figure 1(b)).(iii)
If(δ-τ)/β<δ/β<ET, then the equilibria Zr1 and Zv2 can coexist (shown in Figures 1(c) and 1(d)). Moreover, when ET=0.8, the attractor converges to the Zr1 (see Figure 1(d)).Figure 1
Real/virtual equilibria. The parameter values are fixed as follows:p=0.1,τ=1,σ=0.35,β=4,α=3, and δ=2.0. (a) Zv1 and Zr2 coexist with ET=0.2. (b) Zv1 and Zv2 coexist with ET=0.45. (c) Zr1 and Zv2 coexist with ET=0.52. (d) Zr1 and Zv2 coexist with ET=0.8.
(a) (b) (c) (d)In the process of pest control, we should use IPM strategies to prevent pest outbreaks (i.e.,x(t)<ET). According to the above discussions, it is obvious that β,τ, and ET are key factors in determining the existence of the above different types of equilibria of system (6). So we discuss the β and ET parameter space and β and τ parameter space, respectively, which can be divided into three regions (see Figure 2), and coexistence of real or virtual equilibria is indicated in each region. Figure 2(a) shows the β and ET parameter space with δ=2.0 and τ=1 and Figure 2(b) shows parameter space with ET=0.8 and δ=2. In region I, Zr1 and Zv2 coexist, and Zv1 and Zr2 can coexist in region III. Moreover, It is obvious that the two virtual equilibria Zv1 and Zv2 can coexist in region II. The region II is very important for pest control. In practice, we should apply IPM strategies to prevent multiple pest outbreaks (i.e., x(t)≥ET) in the process of pest control. So, we should choose a set of parameters such that all the equilibria of subsystems SG1 and SG2 are virtual equilibria, which have been widely used in pest control and plant disease [1, 3, 9]. Therefore, from the perspective of pest management, the appropriate control strategies and ET should be designed such that the interior equilibrium of SG1 and SG2 is virtual simultaneously; the region II (purple) is the ideal region to design optimal control strategy as shown in Figure 2.Figure 2
Parameter bifurcation diagram for the existence of equilibria of system (6). Region I (red): Zr1 and Zv2 coexist. Region II (purple): Zv1 and Zv2 coexist. Region III (green): Zv1 and Zr2 coexist. (a) β and ET parameter space with δ=2.0 and τ=1. (b) β and τ parameter space with ET=0.8 and δ=2.
(a) (b)


## 3. Numerical Bifurcation Analysis

In order to show the richer dynamic behavior of discrete switching system (6) aiming to address the implications in biological pest control, we carry out numerical bifurcation investigations in this section.


### 3.1. Bifurcation Diagrams with Different Key Parameters

To discuss the complex dynamics of model (6), we first choose the step size σ as a bifurcation parameter and fix all other parameter values as those in Figure 3, where Figures 3(a) and 3(b) are bifurcation diagrams for prey and predator populations, respectively. They indicate that if the amplitude of the step size σ is increasing from 0.45 to 0.79, the system goes through very complicated behavior, including chaos bands, narrow and wide period-windows, and chaos crises. As the parameter σ increases from 0.45 to 0.68, there exist different periodicity attractors. As σ further increases (i.e., σ∈(0.68,0.79)), the system enters chaotic state. There exists a wide multiattractors coexistence region when σ∈(0.46,0.53) and (0.645,0.67) especially; the details on multiattractors coexistence are discussed in the following section.Figure 3
Bifurcation diagram for system (6):The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.43,β=4,α=3,δ=2.0, and ET=0.6.
(a) (b)


### 3.2. Multiple Attractors, Coexistence, and Initial Sensitivity

From Figure3, we know that several types of attractors coexist. In view of this, we will focus on how the initial densities affect the final states, pest outbreaks, and successful pest control. To confirm this, we fix all parameters as those in Figure 4 and choose different initial densities. For example, two different prey-outbreak attractors coexist at σ=0.66, from which we can find that the two prey-outbreak attractors have different amplitudes and frequencies. Let the initial density of prey and predator be (x0,y0)=(0.4,0.15); then the outbreak patterns for prey population are quite complex, as shown in Figure 4(a), the maximal prey-outbreak amplitude (density) is 0.802, and the system experiences 6 generation and reaches the next maximal outbreak (e.g., the first maximal outbreak lies in 1403 generation, and the next maximal outbreak lies in 1409 generation). If we let the initial values be (x0,y0)=(0.52,0.2), then the maximal prey-outbreak density of system (6) is 0.956 (see Figure 4(c)), and the attractor oscillates with period 7 (e.g., the first maximal outbreak lies in 1402 generation, and the next maximal outbreak lies in 1409 generation). Note that these attractors show that they will go through three different smaller amplitude outbreaks between 1402 generation and 1409 generation (which lies in 1403, 1404, and 1408 generation, resp.).Figure 4
Several different types of typical solutions of system (6). The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4,α=3,δ=2.0, and ET=0.6. (a-b) The initial density of prey and predator (x0,y0)=(0.4,0.15) and (c-d) (x0,y0)=(0.52,0.2). Notice that there exists another attractor that tends to infinity (not shown here).
(a) (b) (c) (d)Furthermore, there exists another attractor that tends to infinity (i.e., the prey and predator densities eventually tend to infinity). This attractor is very harmful for pest control. Moreover, in this case, we can not spray insecticides and release natural enemies to keep prey density below the ET. The above results suggest that the control of insect pests may depend on the initial densities of prey and predator populations. Thus, the initial attraction regions of these stable attractors play a pivotal role in pest management. To show this, the basins of attraction with respect to Figure4 coexistence are given in Figure 5. It indicates that the three attractors can coexist in various prey-predator initial densities. The horizontal axis and the vertical axis are the prey and predator initial values, respectively. In Figure 5(b), the initial value ranges are 0≤x0≤2 and 0≤y0≤2 and Figure 5(a) is an enlargement of the Figure 5(a) with range 0≤x0≤1 and 0≤y0≤1. There exists the basin of attraction for three attractors: the magenta, yellow, and green areas. The magenta and yellow attraction regions are related to the periodic solutions which are shown in Figures 4(a), 4(b), 4(c), and 4(d), respectively. Moreover, the green areas correspond to attractor tending to infinity. Thus, we conclude that the magenta and yellow area are easy to take the integrated pest management strategies for preventing pest (prey) outbreak according to analysis of coexistence attractors shown in Figure 4, whereas the green areas are very harmful for pest control. If the initial densities of prey and predator populations fall into the green regions, the prey and predator populations tend to infinity and can not be controlled by the integrated pest management strategies. Furthermore, we can find the fractal properties of the self-similarity and fractal basin boundaries. Note that the basin of attraction separates the attraction regions into two parts by the line x0=ET (here ET = 0.6), which reveals that different prey-outbreak patterns and control strategies exist there. From pest control point of view, the integrated control strategies may strictly depend on the initial densities of both populations, because the different attractors have different outbreak amplitudes and outbreak frequencies.Figure 5
Basins of attraction of four coexisting attractors of system (6); the horizontal axis and the vertical axis are the prey and predator initial values x0 and y0, respectively. (a) Basin of attraction of two periodic solutions shown in Figure 4 and one infinity solution, where the intervals of initial values are 0≤x0≤1 and 0≤y0≤1; (b) 0≤x0≤2 and 0≤y0≤2. The other parameters are fixed the same as Figure 4. The magenta regions represent Figures 4(a) and 4(b) attractor, the yellow regions are related to Figures 4(c) and 4(d) attractor, and the green region represents infinity attractor.
(a) (b)


## 3.1. Bifurcation Diagrams with Different Key Parameters

To discuss the complex dynamics of model (6), we first choose the step size σ as a bifurcation parameter and fix all other parameter values as those in Figure 3, where Figures 3(a) and 3(b) are bifurcation diagrams for prey and predator populations, respectively. They indicate that if the amplitude of the step size σ is increasing from 0.45 to 0.79, the system goes through very complicated behavior, including chaos bands, narrow and wide period-windows, and chaos crises. As the parameter σ increases from 0.45 to 0.68, there exist different periodicity attractors. As σ further increases (i.e., σ∈(0.68,0.79)), the system enters chaotic state. There exists a wide multiattractors coexistence region when σ∈(0.46,0.53) and (0.645,0.67) especially; the details on multiattractors coexistence are discussed in the following section.Figure 3
Bifurcation diagram for system (6):The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.43,β=4,α=3,δ=2.0, and ET=0.6.
(a) (b)


## 3.2. Multiple Attractors, Coexistence, and Initial Sensitivity

From Figure3, we know that several types of attractors coexist. In view of this, we will focus on how the initial densities affect the final states, pest outbreaks, and successful pest control. To confirm this, we fix all parameters as those in Figure 4 and choose different initial densities. For example, two different prey-outbreak attractors coexist at σ=0.66, from which we can find that the two prey-outbreak attractors have different amplitudes and frequencies. Let the initial density of prey and predator be (x0,y0)=(0.4,0.15); then the outbreak patterns for prey population are quite complex, as shown in Figure 4(a), the maximal prey-outbreak amplitude (density) is 0.802, and the system experiences 6 generation and reaches the next maximal outbreak (e.g., the first maximal outbreak lies in 1403 generation, and the next maximal outbreak lies in 1409 generation). If we let the initial values be (x0,y0)=(0.52,0.2), then the maximal prey-outbreak density of system (6) is 0.956 (see Figure 4(c)), and the attractor oscillates with period 7 (e.g., the first maximal outbreak lies in 1402 generation, and the next maximal outbreak lies in 1409 generation). Note that these attractors show that they will go through three different smaller amplitude outbreaks between 1402 generation and 1409 generation (which lies in 1403, 1404, and 1408 generation, resp.).Figure 4
Several different types of typical solutions of system (6). The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4,α=3,δ=2.0, and ET=0.6. (a-b) The initial density of prey and predator (x0,y0)=(0.4,0.15) and (c-d) (x0,y0)=(0.52,0.2). Notice that there exists another attractor that tends to infinity (not shown here).
(a) (b) (c) (d)Furthermore, there exists another attractor that tends to infinity (i.e., the prey and predator densities eventually tend to infinity). This attractor is very harmful for pest control. Moreover, in this case, we can not spray insecticides and release natural enemies to keep prey density below the ET. The above results suggest that the control of insect pests may depend on the initial densities of prey and predator populations. Thus, the initial attraction regions of these stable attractors play a pivotal role in pest management. To show this, the basins of attraction with respect to Figure4 coexistence are given in Figure 5. It indicates that the three attractors can coexist in various prey-predator initial densities. The horizontal axis and the vertical axis are the prey and predator initial values, respectively. In Figure 5(b), the initial value ranges are 0≤x0≤2 and 0≤y0≤2 and Figure 5(a) is an enlargement of the Figure 5(a) with range 0≤x0≤1 and 0≤y0≤1. There exists the basin of attraction for three attractors: the magenta, yellow, and green areas. The magenta and yellow attraction regions are related to the periodic solutions which are shown in Figures 4(a), 4(b), 4(c), and 4(d), respectively. Moreover, the green areas correspond to attractor tending to infinity. Thus, we conclude that the magenta and yellow area are easy to take the integrated pest management strategies for preventing pest (prey) outbreak according to analysis of coexistence attractors shown in Figure 4, whereas the green areas are very harmful for pest control. If the initial densities of prey and predator populations fall into the green regions, the prey and predator populations tend to infinity and can not be controlled by the integrated pest management strategies. Furthermore, we can find the fractal properties of the self-similarity and fractal basin boundaries. Note that the basin of attraction separates the attraction regions into two parts by the line x0=ET (here ET = 0.6), which reveals that different prey-outbreak patterns and control strategies exist there. From pest control point of view, the integrated control strategies may strictly depend on the initial densities of both populations, because the different attractors have different outbreak amplitudes and outbreak frequencies.Figure 5
Basins of attraction of four coexisting attractors of system (6); the horizontal axis and the vertical axis are the prey and predator initial values x0 and y0, respectively. (a) Basin of attraction of two periodic solutions shown in Figure 4 and one infinity solution, where the intervals of initial values are 0≤x0≤1 and 0≤y0≤1; (b) 0≤x0≤2 and 0≤y0≤2. The other parameters are fixed the same as Figure 4. The magenta regions represent Figures 4(a) and 4(b) attractor, the yellow regions are related to Figures 4(c) and 4(d) attractor, and the green region represents infinity attractor.
(a) (b)


## 4. Switching Times, Switching Frequencies, and Parameter Sensitivity

In this section, we will discuss the effects of key parameters on the switching frequencies of system (6).


### 4.1. Switching Times and Switching Frequencies

For convenience, we provide the definition of switching times and switching frequencies as follows.Definition 4.
For time series of the prey population in system (6), if x(t)≥ET and x(t+1)<ET (or x(t)≤ET and x(t+1)>ET), then one says the system experiences one time switch and t is called switched point. The interval between two switched points is defined as switching times N. The switching frequencies are defined as 1/N.The switching times and switching frequencies are critical factors for successful pest control. If switching times (or switching frequencies) of prey population are a constant number with the time passing, then the system is stable, and we can easily design the control strategy. Moreover, the smaller the switching times (or the larger the switching frequencies) are, the more frequently the control tactics should be applied; that is, the pesticide applications and releasing strategies should be implemented frequently, and this is not effective and may result in adverse effects. On the contrary, if the switching times and switching frequencies are dynamic, the system is unstable or has entered a chaotic state. What is more, it is difficult to design suitable control measures for pest control because we do not know when and how control strategies should be adopted.After that, the switching frequencies are analysed with the spectrogram, which is a visual representation of the spectrum of frequencies in the time series of the prey population as they vary with time. To address the effects of parameters on switching times and switching frequencies, we let the initial values and intrinsic growth ratesα vary, as shown in Figure 6, from which we can see that the switching frequencies have different patterns with different initial values and intrinsic growth rates. For example, the switching times of the prey population are always stable when α=3 with the initial density of prey and predator (x0,y0)=(0.4,0.15), as shown in Figure 6(a). The corresponding spectrogram (see Figure 6(c)) confirms that the time series is periodic. Moreover, the main frequency is about 0.165 Hz, so the switching times are about 1/0.1666=6 and have consistency with Figure 4(a). Similarly, Figure 6(b) shows that the prey population is eventually stable, and switching frequencies are about 0.1428 (see Figure 6(e)). So the switching times is about 1/0.1428=7 which corresponds to Figure 4(c). However, if α=3.5 with the initial density of prey and predator (x0,y0)=(0.8,0.2), the systems is always unstable (see Figure 6(c)), and the spectrogram (Figure 6(f)) shows the disorder frequencies. What is more, it is difficult to design suitable control measures for pest control because we do not know the switching times or switching frequencies.Figure 6
Switching times and switching frequencies; (a–c) switching times and (d–f) switching frequencies. The parameter values are fixed as follows:p=0.1,τ=0.3,σ=0.66,β=4,δ=2.0, and ET=0.6. (a, d) α=3; the initial density of prey and predator (x0,y0)=(0.4,0.15). (b, e) α=3; the initial density of prey and predator (x0,y0)=(0.52,0.2). (c, f) α=3.5; the initial density of prey and predator (x0,y0)=(0.8,0.2); the color bar at the right tells us that the deeper red represents the frequencies where the power spectrum is highest.
(a) (b) (c) (d) (e) (f)In order to show more details of the effects of key parameters on switching frequencies, we first introduce the following definition.Definition 5.
Mean switching times are the mean of all switching times betweeni and i+1999 generation; i is a fixed time point.The effects of key parametersp,ET, and α on the mean switching times are shown in Figure 7. In Figure 7(a), the mean switching times are analysed by varying the parameter α (i.e., α=2,2.2,2.5,2.8,3.2) and keeping the other parameters constant, and five curves for different α can be got. It shows that there is a great oscillating feature if p lies in the interval (0.25,0.6), and the maximum of mean switching times also appears in this interval. It implies that the selection of suitable p may be crucial in prolonging the pest outbreak period. If p increases slightly from 0.6 to 1, the mean switching times of all five curves are gradually smaller and became stable. Moreover, the bigger α is, the smaller the mean switching times are. For example, any point of the blue curve (α=2) is above the black curve (α=3.2). In Figure 7(b), the mean switching times of the prey population concerning the parameter ET (i.e., ET=0.6,0.7,0.75,0.8,0.85) are discussed and the other parameters are fixed. The maximum of mean switching times appear in the interval (0.3,0.6) and the ET become larger; the mean switching times are bigger. Therefore, if we want to prolong mean switching times (pest outbreak period), we should let p be in the interval (0.3,0.6) and select appropriate ET. Most importantly, the higher the killing rate of the insecticide is, the smaller the mean switching times frequently are (i.e., the outbreak of pest should be broken off frequently).Figure 7
Mean switching times versusp. The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4, and δ=2.0; the initial density of prey and predator (x0,y0)=(0.2,0.5). (a) ET=0.6 and α=2,2.2,2.5,2.8,3.2, respectively. (b) α=2 and ET=0.85,0.80,0.75,0.70,0.60, respectively. Mean switching times are calculated between 50 and 2050 generation.
(a) (b)


### 4.2. Sensitivity Analysis of Varying Parameter Values for Mean Switching Times

To examine the sensitivity of mean switching times to parameter variation, we use Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCs). The LHS method is a type of stratified monte carlo sampling, which was first proposed by McKay et al. [10], and was applied to disease transmission [11], integrated pest management [2], HIV/AIDS epidemic [12], cellular network dynamics [13], and systems biology [14]. We investigate the important parameters which affect mean switching times most significantly by using LHS and PRCCs, including the reduction proportion of the pest population (p), the economic threshold (ET), the increasing proportion of natural enemies (τ), the conversion ratio (β), the growth rate (α), and death rate (δ). Here, the mean switching times are calculated between 50 and 2050 generation. We performed uncertainty and sensitivity analyses for all parameters using LHS with 3,000 samples and assumed that all parameters with wide ranges were submitted to uniform distribution, such as p~U(0.1,0.3),ET~U(0.5,0.8),τ~U(0.1,0.3),β~U(3.6,4.1),α~U(1,3), and δ~U(1.8,2.2), and the baseline values of all parameters are given in the Figure 8.Figure 8
PRCC results and PRCC scatter plots. The baseline parameters are fixed as follows:p=0.2,τ=0.2,β=4, and δ=2.0, Here we fix σ=0.25; initial density of prey and predator (x0,y0)=(0.2,0.2). (a) PRCC results. (b–g) PRCC scatter plots. The title of each plot represents the PRCC value with the corresponding p value.
(a) (b) (c) (d) (e) (f) (g)Figure8(a) shows the PRCC results, which illustrates the dependence of mean switching times at each parameter, and PRCC scatter plots for each parameter are given in Figures 8(b)–8(g), respectively. Based on the magnitude of the absolute value of the PRCC, we ranked the relative importance of the 3 input parameters (i.e., the mean switching times are highly dependent on changes in the economic threshold (ET), conversion ratio (β), and the growth rate (α)). Moreover, PRCC results for the mean switching times show that the growth rate (α) has the highest influence on the results and the death rate (δ) is strongly negatively correlated with the mean switching times. On the contrary, the economic threshold (ET) is strongly positively correlated to the mean switching times. The positive sign of their PRCCs indicates that if the parameters increase, the mean switching times increase accordingly (and vice versa). That is, the negative sign suggests that if increased, the mean switching times decrease (and vice versa). Therefore, the parameters ET and β are responsible for the increasing of the mean switching times, so increasing all those parameters is beneficial for prolonging pest outbreak period. This fact is proved numerically that when ET become larger, the mean switching times are bigger in Figure 7(b). On the other side, increasing parameters α, p, τ, and δ will lead to a more serious pest outbreak; we should reduce negative parameters to extend survival pest outbreak period. Figure 7(b) shows that the bigger α is, the smaller the mean switching times are.Figures8(c) and 8(f) show that the ranking of relative importance of parameters could affect the mean switching times. For example, the economic threshold (ET) and the growth rate (α) have relatively high positive and negative influence for pest control, respectively. Moreover, the PRCC values for the reduction proportion of the pest population (p) and the increasing proportion of natural enemies (τ) are small (see Figures 8(b) and 8(g)), which indicate that the killing efficiency rate and the decay rate with respect to natural enemies do not significantly affect the mean switching times.


### 4.3. Global Sensitivity Analysis of Varying Parameter by Heat Map

PRCCs tell us how the mean switching times are affected if we increase (or decrease) a specific parameter at a time point (see Figure7). Here, we focus primarily on time-dependent global sensitivity analysis; we examine the response of mean switching times to parameter variation within a period of time rather than at a time point. We now discuss how to analyse the global sensitivity by using the heat map, which is a graphical representation of data where the individual values contained in a matrix are represented as colors. The developed sensitivity heat map is capable of exploring the sensitivity to many parameters simultaneously and over time [13, 15–17]. To do this, we firstly analyse the PRCCs of mean switching times to identify all those generations (i.e., i is the ith generation; the ith means switching times ranges between i and i+1999 generations, where i=1,2,…,150 generation); then PRCCs of all generations are plotted in the heat map (see Figure 9(b)). On the examination of the heat maps for the mean switching times, we noted that ET was also strongly positive parameter and α was the essentially negative parameter over all generations.Figure 9
Mean switching times versusp; the parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4, and δ=2.0; the initial density of prey and predator (x0,y0)=(0.2,0.5). (a) ET=0.6 and α=2,2.2,2.5,2.8,3.2, respectively. (b) α=2 and ET=0.85,0.80,0.75,0.70,0.60, respectively.
(a) (b)Specifically, in order to show the relative change in the mean switching times over all generations, we performed a box plot to all PRCCs for comparison (see Figure9(a)). The box plot is a convenient way of graphically depicting groups of numerical data and indicated the degree of dispersion (spread) and skewness in the data [18] and was used to analyse the variation of the reproduction numbers to predict the HIV/AIDS epidemic [12]. It shows that the economic threshold (ET) and the growth rate (α) have very little variation over all generations. Therefore, the economic threshold (ET) and the growth rate (α) are very more important parameters for pest control.


## 4.1. Switching Times and Switching Frequencies

For convenience, we provide the definition of switching times and switching frequencies as follows.Definition 4.
For time series of the prey population in system (6), if x(t)≥ET and x(t+1)<ET (or x(t)≤ET and x(t+1)>ET), then one says the system experiences one time switch and t is called switched point. The interval between two switched points is defined as switching times N. The switching frequencies are defined as 1/N.The switching times and switching frequencies are critical factors for successful pest control. If switching times (or switching frequencies) of prey population are a constant number with the time passing, then the system is stable, and we can easily design the control strategy. Moreover, the smaller the switching times (or the larger the switching frequencies) are, the more frequently the control tactics should be applied; that is, the pesticide applications and releasing strategies should be implemented frequently, and this is not effective and may result in adverse effects. On the contrary, if the switching times and switching frequencies are dynamic, the system is unstable or has entered a chaotic state. What is more, it is difficult to design suitable control measures for pest control because we do not know when and how control strategies should be adopted.After that, the switching frequencies are analysed with the spectrogram, which is a visual representation of the spectrum of frequencies in the time series of the prey population as they vary with time. To address the effects of parameters on switching times and switching frequencies, we let the initial values and intrinsic growth ratesα vary, as shown in Figure 6, from which we can see that the switching frequencies have different patterns with different initial values and intrinsic growth rates. For example, the switching times of the prey population are always stable when α=3 with the initial density of prey and predator (x0,y0)=(0.4,0.15), as shown in Figure 6(a). The corresponding spectrogram (see Figure 6(c)) confirms that the time series is periodic. Moreover, the main frequency is about 0.165 Hz, so the switching times are about 1/0.1666=6 and have consistency with Figure 4(a). Similarly, Figure 6(b) shows that the prey population is eventually stable, and switching frequencies are about 0.1428 (see Figure 6(e)). So the switching times is about 1/0.1428=7 which corresponds to Figure 4(c). However, if α=3.5 with the initial density of prey and predator (x0,y0)=(0.8,0.2), the systems is always unstable (see Figure 6(c)), and the spectrogram (Figure 6(f)) shows the disorder frequencies. What is more, it is difficult to design suitable control measures for pest control because we do not know the switching times or switching frequencies.Figure 6
Switching times and switching frequencies; (a–c) switching times and (d–f) switching frequencies. The parameter values are fixed as follows:p=0.1,τ=0.3,σ=0.66,β=4,δ=2.0, and ET=0.6. (a, d) α=3; the initial density of prey and predator (x0,y0)=(0.4,0.15). (b, e) α=3; the initial density of prey and predator (x0,y0)=(0.52,0.2). (c, f) α=3.5; the initial density of prey and predator (x0,y0)=(0.8,0.2); the color bar at the right tells us that the deeper red represents the frequencies where the power spectrum is highest.
(a) (b) (c) (d) (e) (f)In order to show more details of the effects of key parameters on switching frequencies, we first introduce the following definition.Definition 5.
Mean switching times are the mean of all switching times betweeni and i+1999 generation; i is a fixed time point.The effects of key parametersp,ET, and α on the mean switching times are shown in Figure 7. In Figure 7(a), the mean switching times are analysed by varying the parameter α (i.e., α=2,2.2,2.5,2.8,3.2) and keeping the other parameters constant, and five curves for different α can be got. It shows that there is a great oscillating feature if p lies in the interval (0.25,0.6), and the maximum of mean switching times also appears in this interval. It implies that the selection of suitable p may be crucial in prolonging the pest outbreak period. If p increases slightly from 0.6 to 1, the mean switching times of all five curves are gradually smaller and became stable. Moreover, the bigger α is, the smaller the mean switching times are. For example, any point of the blue curve (α=2) is above the black curve (α=3.2). In Figure 7(b), the mean switching times of the prey population concerning the parameter ET (i.e., ET=0.6,0.7,0.75,0.8,0.85) are discussed and the other parameters are fixed. The maximum of mean switching times appear in the interval (0.3,0.6) and the ET become larger; the mean switching times are bigger. Therefore, if we want to prolong mean switching times (pest outbreak period), we should let p be in the interval (0.3,0.6) and select appropriate ET. Most importantly, the higher the killing rate of the insecticide is, the smaller the mean switching times frequently are (i.e., the outbreak of pest should be broken off frequently).Figure 7
Mean switching times versusp. The parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4, and δ=2.0; the initial density of prey and predator (x0,y0)=(0.2,0.5). (a) ET=0.6 and α=2,2.2,2.5,2.8,3.2, respectively. (b) α=2 and ET=0.85,0.80,0.75,0.70,0.60, respectively. Mean switching times are calculated between 50 and 2050 generation.
(a) (b)


## 4.2. Sensitivity Analysis of Varying Parameter Values for Mean Switching Times

To examine the sensitivity of mean switching times to parameter variation, we use Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficients (PRCCs). The LHS method is a type of stratified monte carlo sampling, which was first proposed by McKay et al. [10], and was applied to disease transmission [11], integrated pest management [2], HIV/AIDS epidemic [12], cellular network dynamics [13], and systems biology [14]. We investigate the important parameters which affect mean switching times most significantly by using LHS and PRCCs, including the reduction proportion of the pest population (p), the economic threshold (ET), the increasing proportion of natural enemies (τ), the conversion ratio (β), the growth rate (α), and death rate (δ). Here, the mean switching times are calculated between 50 and 2050 generation. We performed uncertainty and sensitivity analyses for all parameters using LHS with 3,000 samples and assumed that all parameters with wide ranges were submitted to uniform distribution, such as p~U(0.1,0.3),ET~U(0.5,0.8),τ~U(0.1,0.3),β~U(3.6,4.1),α~U(1,3), and δ~U(1.8,2.2), and the baseline values of all parameters are given in the Figure 8.Figure 8
PRCC results and PRCC scatter plots. The baseline parameters are fixed as follows:p=0.2,τ=0.2,β=4, and δ=2.0, Here we fix σ=0.25; initial density of prey and predator (x0,y0)=(0.2,0.2). (a) PRCC results. (b–g) PRCC scatter plots. The title of each plot represents the PRCC value with the corresponding p value.
(a) (b) (c) (d) (e) (f) (g)Figure8(a) shows the PRCC results, which illustrates the dependence of mean switching times at each parameter, and PRCC scatter plots for each parameter are given in Figures 8(b)–8(g), respectively. Based on the magnitude of the absolute value of the PRCC, we ranked the relative importance of the 3 input parameters (i.e., the mean switching times are highly dependent on changes in the economic threshold (ET), conversion ratio (β), and the growth rate (α)). Moreover, PRCC results for the mean switching times show that the growth rate (α) has the highest influence on the results and the death rate (δ) is strongly negatively correlated with the mean switching times. On the contrary, the economic threshold (ET) is strongly positively correlated to the mean switching times. The positive sign of their PRCCs indicates that if the parameters increase, the mean switching times increase accordingly (and vice versa). That is, the negative sign suggests that if increased, the mean switching times decrease (and vice versa). Therefore, the parameters ET and β are responsible for the increasing of the mean switching times, so increasing all those parameters is beneficial for prolonging pest outbreak period. This fact is proved numerically that when ET become larger, the mean switching times are bigger in Figure 7(b). On the other side, increasing parameters α, p, τ, and δ will lead to a more serious pest outbreak; we should reduce negative parameters to extend survival pest outbreak period. Figure 7(b) shows that the bigger α is, the smaller the mean switching times are.Figures8(c) and 8(f) show that the ranking of relative importance of parameters could affect the mean switching times. For example, the economic threshold (ET) and the growth rate (α) have relatively high positive and negative influence for pest control, respectively. Moreover, the PRCC values for the reduction proportion of the pest population (p) and the increasing proportion of natural enemies (τ) are small (see Figures 8(b) and 8(g)), which indicate that the killing efficiency rate and the decay rate with respect to natural enemies do not significantly affect the mean switching times.


## 4.3. Global Sensitivity Analysis of Varying Parameter by Heat Map

PRCCs tell us how the mean switching times are affected if we increase (or decrease) a specific parameter at a time point (see Figure7). Here, we focus primarily on time-dependent global sensitivity analysis; we examine the response of mean switching times to parameter variation within a period of time rather than at a time point. We now discuss how to analyse the global sensitivity by using the heat map, which is a graphical representation of data where the individual values contained in a matrix are represented as colors. The developed sensitivity heat map is capable of exploring the sensitivity to many parameters simultaneously and over time [13, 15–17]. To do this, we firstly analyse the PRCCs of mean switching times to identify all those generations (i.e., i is the ith generation; the ith means switching times ranges between i and i+1999 generations, where i=1,2,…,150 generation); then PRCCs of all generations are plotted in the heat map (see Figure 9(b)). On the examination of the heat maps for the mean switching times, we noted that ET was also strongly positive parameter and α was the essentially negative parameter over all generations.Figure 9
Mean switching times versusp; the parameter values are fixed as follows: p=0.1,τ=0.3,σ=0.66,β=4, and δ=2.0; the initial density of prey and predator (x0,y0)=(0.2,0.5). (a) ET=0.6 and α=2,2.2,2.5,2.8,3.2, respectively. (b) α=2 and ET=0.85,0.80,0.75,0.70,0.60, respectively.
(a) (b)Specifically, in order to show the relative change in the mean switching times over all generations, we performed a box plot to all PRCCs for comparison (see Figure9(a)). The box plot is a convenient way of graphically depicting groups of numerical data and indicated the degree of dispersion (spread) and skewness in the data [18] and was used to analyse the variation of the reproduction numbers to predict the HIV/AIDS epidemic [12]. It shows that the economic threshold (ET) and the growth rate (α) have very little variation over all generations. Therefore, the economic threshold (ET) and the growth rate (α) are very more important parameters for pest control.


## 5. Discussion

It is well documented that the threshold policy (or on-off policy, switching policy) is commonly used in biological system, which is to allow removal of the prey (pest) or increase the predator (natural enemies), such as on-off policy in the herbivore-vegetation dynamics [19], endogenous switching of harvesting strategy in prey-predator fishery model [20], gause prey-predator model with a refuge [21], and two stage pest control models with economic thresholds [1]. In this study, we discuss discrete-time prey-predator models by utilizing a threshold policy (TP) to control the pest population (prey). This switching system is divided into two subsystems, namely, the control-free subsystem SG1 and controlled subsystem SG2. The existence and stability of several types of equilibria of the full switching system are discussed; the relationship between the real and virtual equilibrium for two subsystems is showed as in Figure 1. In order to prevent multiple pest outbreaks (i.e., x(t)<ET), the key parameters of equilibria of system (6) should be analysed. The parameter bifurcation diagram for the Equilibria in the β and ET parameter space and β and τ parameter space is shown as in Figure 2 respectively, which are divided into three regions. Furthermore, the region II is very important for pest control.The bifurcation diagram of parameterσ provided the evidence that the switching discrete system may have very complex dynamics including coexistence of multiple attractors among attractors (see Figure 3). To confirm the coexistence of multiple attractors and discuss their biological implications, Figure 4 showed that the control of insect pests may depend on initial densities of prey and predator populations. Thus, the initial attraction regions of these stable attractors play a pivotal role in pest management as shown in Figure 5. It indicates that the three attractors can coexist in various prey-predator initial densities.The key parameters concerning switching frequencies or mean switching times have been analysed, and consequently the relative biological implications with respect to pest control are discussed. If the switching frequencies and switching times are always unstable, it is difficult to design suitable control measures for pest control because we do not know when and how control strategies should be adopted. In order to show more details of the effects of key parameters on unstable states, we first introduce the mean switching times. The effects of key parametersp,ET, and α on the mean switching times are shown in Figure 7. Most importantly, (i) the bigger α is, the smaller the mean switching times are, (ii) the higher the killing rate of the insecticide is, and the smaller the mean switching times frequently are (i.e., the outbreak of pest should be broken off frequently).To examine the sensitivity of mean switching times to parameter variation, we used Latin Hypercube Sampling (LHS) and partial rank correlation coefficients (PRCCs). The results show that the economic threshold (ET) is strongly positively correlated to the mean switching times and the growth rate (α) is strongly negatively correlated with the mean switching times (Figures 8(b)–8(g)). The global sensitivity using the heat map and PRCCs proved that economic threshold (ET) and the growth rate (α) were the key parameters for pest control.The numerical bifurcation analyses clarify that the switching system could have very complex dynamics including multiple attractor coexistence and chaotic solutions. From the pest control point of view, we have carried out extensively numerical investigations on the switching times and their biological implications have been discussed in more detail. The theoretical/mathematical analysis about subsystemSG1 or SG2 has been studied extensively. Thus, it is very important to select the corresponding parameters for numerical investigation. In order to consider the effects of parameter selection, the global sensitivity analysis of all parameters has been discussed.



---

*Source: 8627613-2016-12-08.xml*