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.
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 aswhere 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 , and . The discrete-time prey-predator system can be obtained by the forward Euler scheme to system (2):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 based on (3) is as follows:where 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 . Furthermore, if , we only release sufficient natural enemies to prevent the pest population from exceeding ET. If and , then we adopt the chemical control strategy so that it can ensure pest population decrease when the pest population reaches the ET.
By , the nonnegative octant could be split in two parts:
In the region the more profitable prey density is below ET in system (3). In the region , 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 and system (4) in region so systems (6) in region and in region are divided into two different subsystems by utilizing threshold policy (ET), which are denoted as and , 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 of the plane. We first discuss the existence of fixed points for subsystems and , 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 Subsystems and
For the subsystem , we are interested in the positive point; let . By a simple computation, it is straightforward to obtain the following results.
If , the subsystem has a unique positive fixed point . The stability of the fixed point has been discussed by Liu and Xiao [4].
We have similar analysis for subsystem . It is clear that subsystem has unique positive fixed point , 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].
Lemma 1. Let . Suppose that and are two roots of . Then(i), and if and only if and ;(ii), and (or , and ) if and only if ;(iii), and if and only if and ;(iv), and if and only if and ;(v) and are complex and and if and only if and .
Note that the local stability of the fixed point is determined by the modules of eigenvalues of the characteristic equation. The Jacobian matrix of the subsystem at is given byand the characteristic equation of the Jacobian matrix can be written as
The characteristic polynomial is denoted as where
Then and .
Let and are two roots of the characteristic equation (8), which are called eigenvalues of the fixed point . The fixed point is called a sink if and ; the sink is locally asymptotically stable. is called a source if and ; the source is locally unstable. is called a saddle if and (or and ). is called nonhyperbolic if either or .
Then we can calculate the local dynamics of the fixed point by Lemma 1.
Proposition 2. Let be the positive fixed point of (8).(i)It is a sink if one of the following conditions holds:(i.1) and ;(i.2) and .(ii)It is a source if one of the following conditions holds:(ii.1) and ;(ii.2) and .(iii)It is a saddle if the following condition holds: and (iv)It is nonhyperbolic if the following condition holds: and and .
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. Let be the stable equilibrium point of the dynamics of region for . Then is called a real equilibrium if it belongs to subsystem and a virtual equilibrium if it belongs to , , for . The real equilibrium is denoted as and and the virtual equilibrium points are denoted as and , respectively.
For the subsystem with (i.e., it belongs to region ), if and , then it is a real equilibrium for the system , denoted by . If and , then it is a virtual equilibrium for the system , denoted by . About real equilibria for the subsystem with (i.e., it belongs to region ), on the basis of coordinate of equilibrium , we have the following results. If , and , then it is a real equilibrium for the system , denoted by . If , and , then it is a virtual equilibrium for the system , denoted by .
In order to consider the relationship of equilibria for two subsystems, we define five curves as follows: , and . Based on the above analysis and , we consider the following three cases.(i)If , the equilibria and can coexist as shown in Figure 1(a).(ii)If , the equilibria and can coexist (i.e., Figure 1(b)).(iii)If , then the equilibria and can coexist (shown in Figures 1(c) and 1(d)). Moreover, when , the attractor converges to the (see Figure 1(d)).
(a)
(b)
(c)
(d)
In the process of pest control, we should use IPM strategies to prevent pest outbreaks (i.e., ). According to the above discussions, it is obvious that , and are key factors in determining the existence of the above different types of equilibria of system (6). So we discuss the and 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 parameter space with and and Figure 2(b) shows parameter space with and In region I, and coexist, and and can coexist in region III. Moreover, It is obvious that the two virtual equilibria and 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., ) in the process of pest control. So, we should choose a set of parameters such that all the equilibria of subsystems and 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 and is virtual simultaneously; the region II (purple) is the ideal region to design optimal control strategy as shown in Figure 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 to , there exist different periodicity attractors. As further increases (i.e., ), the system enters chaotic state. There exists a wide multiattractors coexistence region when and especially; the details on multiattractors coexistence are discussed in the following section.
(a)
(b)
3.2. 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 , 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 ; 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 , 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.).
(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 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 and and Figure 5(a) is an enlargement of the Figure 5(a) with range and . 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 (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.
(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 and (or and ), 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 .
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 with the initial density of prey and predator , 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 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 which corresponds to Figure 4(c). However, if with the initial density of prey and predator , 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.
(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 between and generation; is a fixed time point.
The effects of key parameters , 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., ) 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 , 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 () is above the black curve (). In Figure 7(b), the mean switching times of the prey population concerning the parameter ET (i.e., ) are discussed and the other parameters are fixed. The maximum of mean switching times appear in the interval 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 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).
(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 (), 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 and 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 , and , and the baseline values of all parameters are given in the Figure 8.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
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 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.
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 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–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 generations, where 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.
(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 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.
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 and controlled subsystem . 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., ), the key parameters of equilibria of system (6) should be analysed. The parameter bifurcation diagram for the Equilibria in the and 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 parameters , 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 or 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.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Science Foundation of Hubei Province of China (2015CFB264), the Key Laboratory of Biologic Resource Protection and Utilization of Hubei Province (PKLHB1506, Changcheng Xiang), and the National Natural Science Foundation of China (no. 11601268, Wenjie Qin).