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

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 ofmultiple 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.


Introduction
The integrated pest management (IPM) is applied [1][2][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 where () and () 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: where  = ,  =  2 , and  = .The discrete-time prey-predator system can be obtained by the forward Euler scheme to system (2):

𝑥 (𝑡 + 1) = 𝑥 (𝑡) + 𝜎 [𝛼𝑥 (𝑡) (1 − 𝑥 (𝑡)) − 𝛽𝑥 (𝑡) 𝑦 (𝑡)] ,
( + 1) =  () +  (− +  ())  () , 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 with () ≥ ET based on (3) is as follows: ( + 1) = (1 − ) ( () +  [ () (1 −  ()) −  ()  ()]) ,  ( + 1) =  () +  (− +  ())  () +  () ,  () ≥ ET, (4) where 0 <  < 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  = 0, we only release sufficient natural enemies to prevent the pest population from exceeding ET.If  > 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.
By  2 + = {((), ()) | () ≥ 0, () ≥ 0}, the nonnegative octant could be split in two parts: In the region  1 the more profitable prey density is below ET in system (3).In the region  2 , 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.

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  2 + of the (, ) plane.We first discuss the existence of fixed points for subsystems   1 and   2 , respectively, and then study the switching behavior of fixed points by utilizing a threshold policy.

Existence and Locally Asymptotic Stability of Fixed Points
for the Subsystems   1 and   2 .For the subsystem   1 , we are interested in the positive point; let ( + 1) = () =  * , ( + 1) = () =  * .By a simple computation, it is straightforward to obtain the following results.
We have similar analysis for subsystem   2 .It is clear that subsystem   2 has unique positive fixed point ( * 2 ,  * 2 ) = (( − )/, (( − ) +  − )/ 2 ), if  >  and ( − ) >  − .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].
In the process of pest control, we should use IPM strategies to prevent pest outbreaks (i.e., () < 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 V and  2 V 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., () ≥ ET) in the process of pest control.So, we should choose a set of parameters such that all the equilibria of subsystems   1 and   2 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   1 and   2 is virtual simultaneously; the region II (purple) is the ideal region to design optimal control strategy as shown in Figure 2.

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.

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.

Multiple Attractors, Coexistence, and Initial Sensitivity.
From Figure 3, 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 ( 0 ,  0 ) = (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 ( 0 ,  0 ) = (0.52, 0.2), then the maximal preyoutbreak 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.).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 Figure 4 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 ≤  0 ≤ 2 and 0 ≤  0 ≤ 2 and Figure 5(a) is an enlargement of the Figure 5(a) with range 0 ≤  0 ≤ 1 and 0 ≤  0 ≤ 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

Coexist Coexist
Coexist 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 selfsimilarity and fractal basin boundaries.Note that the basin of attraction separates the attraction regions into two parts by the line  0 = 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.

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).

Switching Times and Switching
Frequencies.For convenience, we provide the definition of switching times and switching frequencies as follows.Figure 4: Several different types of typical solutions of system (6).The parameter values are fixed as follows:  = 0.1,  = 0.3,  = 0.66,  = 4,  = 3,  = 2.0, and ET = 0.6.(a-b) The initial density of prey and predator ( 0 ,  0 ) = (0.4,0.15) and (c-d) ( 0 ,  0 ) = (0.52, 0.2).Notice that there exists another attractor that tends to infinity (not shown here).Definition 4. For time series of the prey population in system (6), if () ≥ ET and ( + 1) < ET (or () ≤ ET and (+1) > ET), then one says the system experiences one time switch and  is called switched point.The interval between two switched points is defined as switching times .The switching frequencies are defined as 1/.
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 ( 0 ,  0 ) = (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 ( 0 ,  0 ) = (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.
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 between  and  + 1999 generation;  is a fixed time point.
The effects of key parameters , 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  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  may be crucial in prolonging the pest outbreak period.If  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  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).

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 (), 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  ∼ (0.1, 0.3), ET ∼ (0.5, 0.8),  ∼ (0.1, 0.3),  ∼ (3.6,4.1),  ∼ (1, 3), and  ∼ (1.8,2.2), and the baseline values of all parameters are given in the Figure 8. Figure 8(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 , , , 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.
Figures 8(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 () 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.

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 Figure 7).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][16][17].To do this, we firstly analyse the PRCCs of mean switching times to identify all those generations (i.e.,  is the th generation; the th means switching times ranges between  and  + 1999 generations, where  = 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.
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 Figure 9(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.

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   1 and controlled subsystem   2 .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., () < 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 preypredator 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 parameters , 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 subsystem   1 or   2 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.

2 V
(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,  1  and  coexist, and  1 V and  2  can coexist in region III.Moreover, It is obvious that the two virtual equilibria  1

Figure 5 :Figure 6 :
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  0 and  0 , 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 ≤  0 ≤ 1 and 0 ≤  0 ≤ 1; (b) 0 ≤  0 ≤ 2 and 0 ≤  0 ≤ 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.