Complex dynamics and coexistence of period-doubling and period-halving bifurcations in an integrated pest management model with nonlinear impulsive control

An expectation for optimal integrated pest management is that the instantaneous numbers of natural enemies released should depend on the densities of both pest and natural enemy in the field. For this, a generalised predator–prey model with nonlinear impulsive control tactics is proposed and its dynamics is investigated. The threshold conditions for the global stability of the pest-free periodic solution are obtained based on the Floquet theorem and analytic methods. Also, the sufficient conditions for permanence are given. Additionally, the problem of finding a nontrivial periodic solution is confirmed by showing the existence of a nontrivial fixed point of the model’s stroboscopic map determined by a time snapshot equal to the common impulsive period. In order to address the effects of nonlinear pulse control on the dynamics and success of pest control, a predator–prey model incorporating the Holling type II functional response function as an example is investigated. Finally, numerical simulations show that the proposed model has very complex dynamical behaviour, including period-doubling bifurcation, chaotic solutions, chaos crisis, period-halving bifurcations and periodic windows. Moreover, there exists an interesting phenomenon whereby period-doubling bifurcation and period-halving bifurcation always coexist when nonlinear impulsive controls are adopted, which makes the dynamical behaviour of the model more complicated, resulting in difficulties when designing successful pest control strategies.


Introduction
In order to reduce losses due to insect pests, a wide range of pest control strategies are available to farmers [1,2]. Integrated pest management (IPM) is a long-term management strategy involving a combination of biological, cultural and chemical tactics aiming to minimise economic, health and environmental risks [3,4].
An important component of an IPM strategy is biological control, which is defined as using natural enemies such as predators and parasites to suppress pest populations and typically involves an active human role [2]. For example, one known approach of biological control is augmentation, i.e. releasing beneficial natural enemies to control insects to tolerable levels [5][6][7]. Chemical control is also an important method for IPM, being very useful because it can kill a significant proportion of an insect population quickly, and sometimes it provides the only feasible method for preventing economic loss. In practical application, in order to avoid or delay the development of resistance, a variety of different types of pesticides are used instead of a single one. Moreover, different pest control techniques should function together rather than be antagonistic to each other. Experience with use of IPM has shown it to be more effective than classical control methods such as biological control or chemical control applied alone [2].
With the development of the theory and application of impulsive differential equations [8,9], mathematical models can assist in the design of IPM strategies and the study of how to control pest populations at tolerable levels. Recently, many researchers have constructed impulsive differential equations to describe IPM strategies [10][11][12][13][14][15]. However, the major assumption in previous studies is that a proportion of the pest population will be killed instantly after spraying pesticide, while simultaneously releasing a constant number of natural enemies [16][17][18][19][20][21][22]. In reality, every community or country has an appropriate or limited capacity for agricultural resources such as pesticides, labour force, equipment and excessive use of the pesticides will potentially cause damage to the environment. So, the methods used for releasing natural enemies and the ratios of numbers of natural enemies to be released to their current density in the field could be significantly affected by limitations to agricultural resources. In order to understand how such a resource limitation could affect an IPM strategy, some mathematical models incorporating saturation of the limited resources have been proposed and analysed recently. These have mainly focused on investigating the global dynamic and addressing how nonlinear impulsive control actions affect pest control tactics [23][24][25][26]. However, the nonlinear impulsive functions in these papers only depend on the density of the natural enemy population. A more realistic case modelled in this paper is that the instantaneous number of natural enemies to be released should depend on the densities of both pest and natural enemy densities in the field, i.e. the density dependent releasing function is a saturation function of the pest population. In this case, the lower the pest population in the field, the lower the number of the predators that should be released and vice versa.
The main purpose of this paper is to construct a generalised predator-prey model with nonlinear impulsive control to investigate the effect of limited predator releases on the outbreak of pest populations. Other generalised predator-prey models with linear or constant pulse actions have already provided some analytical techniques to deal with the global stability of periodic solutions [27][28][29]. In Sect. 2, by using the Floquet theorem and analytical techniques, we show that there exists a globally stable pest-free periodic solution when the period of impulsive predator releases is less than some threshold, and the condition for the permanence of the system is given in Sect. 3. Next, by employing an operator theoretic approach which reduces the existence of the nontrivial periodic solutions to a fixed point and bifurcation problem [30][31][32][33], we show that a nontrivial periodic solution appears via a supercritical bifurcation once the pest-free periodic solution loses its stability. In order to apply the main results, we chose a classic pest-natural enemy model with a Holling type II functional response function and nonlinear impulsive control to investigate how the nonlinear pulse perturbations affect success or otherwise of the pest control, which is presented in Sect. 5. Further, numerical simulations show that the model with nonlinear impulsive actions has more complex and interesting dynamic behaviour than models without them. Finally, related biological implications are discussed.

The predator-prey model with nonlinear pulse control
The generalised predator-prey model we consider in this paper is as follows: where x(t), y(t) represent the densities of prey and predator populations, respectively. r is the intrinsic growth rate of the prey, K denotes the carrying capacity of the prey, c represents the efficiency rate for predator conversion of its prey, D is the death rate of the predator. The function p(x) denotes the predator response function, and let F(x) = rx(1 -x K )/p(x). In order to use our main results for as many models as possible, we made the following assumptions related to the function p(x). Let p(x) be locally Lipschitz functions such that: (i) The functional response function satisfies p(x) ∈ C 2 (0, +∞), p(0) = 0, p (0) > 0 and p (x) < 0 for all x > 0, (ii) The function F(x) and p(x) x is upper bounded for all x > 0. The first condition indicates that the functional response function is positive and monotonically increasing for small pest populations. The p(x) is bounded or linear, so p(x) x is upper bounded for all x > 0. The pest population cannot increase infinitely once the biological control is introduced, so F(x) is upper bounded. It is assumed that the IPM strategy is applied at every time point nT . Moreover, the nonlinear saturation functions or density dependent functions are employed, i.e. we choose where 0 ≤ δ ≤ 1 and h ≥ 0 represent the maximal mortality rate and the half saturation constant for the pest due to chemical control, respectively. We assume that the pesticide is sprayed immediately before releases of natural enemies to minimise the chance of the pesticide killing newly released predators. λ 1 ≥ 0 is the number of the natural enemies released, and θ ≥ 0 denotes the shape parameter. The pest and natural enemy populations are updated to (1 -δx(t) x(t)+h )x(t) and y(t) + λ 1 x(t) 1+θx(t) + λ 2 at every discrete time point nT and n ∈ N , respectively, which is a more realistic consideration compared with previous studies [25,26].
On the basis of the above assumptions and taking model (1) into account, the following differential equation model was formulated: Before our main results, we show the positivity and boundedness of solutions of model (2), and we have the following.

Lemma 1 The solutions of model (2) are positive and bounded.
Proof It is easy to show that the positive initial conditions indicate the positivity of the solutions, and x(t) < max{K, x(0 + )} according to x(nT + ) < x(nT).
For the boundedness of the y-component, we denote V (t) = cx(t)+y(t), then when t = nT we obtain . When t = nT we have Therefore, for t ∈ (nT, (n + 1)T], we have Thus, V (t) is uniformly ultimately bounded. According to the definition of V (t), we can see that y(t) is bounded.
One of the main purposes of implementing IPM is to find suitable control measures to eradicate the pest population, i.e. if there exists a critical value such that pest-free periodic solution is globally stable when the period T is less than the critical value. The existence and stability of the pest-free periodic solution are crucial for this point. Thus, we first consider the following subsystem: Subsystem (3) is a simple linear growth model, so we can easily obtain the following main result.
Therefore, we obtain the general expression of the unique pest-free periodic solution of model (2), i.e.
x p (t), y p (t) = 0, y * exp -D(t -nT) , t ∈ nT, (n + 1)T , In the following, we show the sufficient condition for the global stability of pest-free periodic solution (x p (t), y p (t)) of model (2) and the main results for model (2).
and it is globally attractive if where M s = sup x≥0 F(x).
Proof Firstly, we prove the local stability of the solution (x p (t), y p (t)) of model (2). Let (t) be the fundamental matrix of (2), and then (t) satisfies , where (0) = I represents the identity matrix and the term * is not calculated in the exact form for the next analyses. The linearization of the third and fourth equation of model (2) becomes x(nT + ) Therefore, the local stability of the solution (x p (t), y p (t)) of model (2) is determined by the absolute values of both eigenvalues of the following matrix: The two Floquet multipliers are as follows: All these results confirm that the pest-free solution (x p (t), y p (t)) is locally stable if and only if In the following, we prove the global attractivity of the periodic solution (x p (t), y p (t)), for which we only need to show that any solution (x(t), y(t)) of model (2) tends to (x p (t), y p (t)) as t goes to infinity. If T < λ 2 M s D , then we can choose ε > 0 sufficiently small such that According to the theory of differential equations, we have the following comparison equation: for t large enough. Without loss of generality, we assume that y(t) ≥ y p (t)ε holds true for all t ≥ 0. From the first equation of model (2) Integrating any interval (nT, (n + 1)T], we have For any t ∈ (lT, (l + 1)T], l ∈ Z + , we get For the boundedness of the second term of the right-hand side, it is clear that x(t) converges to 0 as t → ∞.
Next, we prove that y(t) → y p (t) as t → ∞. Since x(t) goes to zero as t → ∞ and assumption (i), there exists a finite time t 1 such that p(x) < ε and λ 1 x(t) 1+θx(t) < ελ 1 for all t > t 1 . Therefore, for all t > t 1 , we deduce that Thus, we consider the following equation: From the comparison theorem of impulsive differential equations, we have i.e. the pest-free periodic solution of model (2) is globally asymptotically stable.

Permanence
Now, we investigate the permanence of model (2).
Proof By Lemma 1, for simplicity, we may assume that (6), we know that y(t) ≥ y p (t)ε for t large enough and for some ε > 0, so we have for t large enough. Thus, we only need to find m 1 such that x(t) ≥ m 1 for t large enough. We will do this in the following two steps.
x . Suppose that x(t) < m 3 for all t > 0, then we have By Lemma 2, we have y(t) ≤ z(t) and and z p (t) = According to assumption (i), we have p(x) < p (0)x, further, for t > T 1 . Let N ∈ Z + and NT > T 1 . Integrating equation (10) on (nT, (n + 1)T], n ≥ N , we have for t ∈ (nT, (n + 1)T] and n 1 ≤ n ≤ n 1 + n 2 + n 3 . Then Furthermore, we have Thus, we have which is a contradiction. (11) we have . For t >t, the same arguments can be continued since x(t) ≥ m 3 .
Remark Define T * = λ 2 p (0) rD , the pest-free periodic solution is unstable if T > T * . Therefore, T * plays a role in a bifurcation threshold value. The interesting question is how to determine the dynamics of model (2) once the pest-free periodic solution loses its stability.

Threshold condition of bifurcation
In the following, we deal with problem of the bifurcation of a nontrivial periodic solution of model (2) near the pest-free periodic solution, i.e. (x p (t), y p (t)). We have used the bifurcation theory in earlier publications [34,35]. We denote to be the solution of the first two equations of model (2), and we have X(t) = (t, X 0 ) where X 0 = X(0) = (y + 0 , 0) and = ( 1 , 2 ). We define the mapping 1 , 2 : R 2 → R 2 by 1 (x, y) = x + λ 1 y 1 + θ y + λ 2 , 2 (x, y) = 1 -δy y + h y and the mapping F 1 , F 2 : R 2 → R 2 as follows: Furthermore, let us define : It is clear that is the stroboscopic map associated with model (2), which is determined by the values of solutions at impulsive points 0 and T, and T is the stroboscopic time snapshot.
We reduce the problem of finding a periodic solution to a fixed problem. Hence, X = (y p (t), 0) is a periodic solution of period T for model (2) if and only if its initial value X 0 = X(0) is a fixed point for map (T, ·). Consequently, in order to establish the existence of nontrivial periodic solutions of model (2), we need to prove the existence of the nontrivial fixed points of .
According to the variational equations which characterise the dynamics of the first two equations in model (2), we get Then we deduce that model (2) has the particular form d dt with initial value D X ( (0, X 0 )) = I 2 , which is the identity matrix.
According to the initial value, it follows that where ρ(t) = exp t 0 (rp (0)y p (v)) dv. We denote the derivation of N at the fixed point (τ ,X) = (0, (0, 0)) as So, we obtain It is easy to see that the necessary condition for the bifurcation of the nontrivial solution is which can reduce to d 0 = 0 because D X N(0, (0, 0)) is an upper triangular matrix with a 0 > 0. Further, we deduce that d 0 = 0 is equivalent to T = T * . Therefore, it remains to prove the sufficient condition for the bifurcation of nontrivial solution. We have the following main result.
Solving the above equation, we have ã τ ≈ -B C > 0, which implies that there is a supercritical bifurcation to a nontrivial periodic solution.
In order to verify our main results, we choose the Holling type II functional response function as an example in the coming section.

Application of the main results
In order to show the applications of the main results and discuss the biological implications of the threshold conditions, we choose the Holling type II functional response function for p(x) = αx 1+ωx . Thus model (2) becomes the following special system: where α and ω are positive constants, respectively. It follows from Lemma 2 that we obtain the pest-free periodic solution of model (17) as follows: x p (t), y p (t) = 0, y * exp -D(t -nT) , t ∈ nt, (n + 1T) with y * = λ 2 1-exp(-DT) . For the stability of the pest-free periodic solution, we employ the main results shown in Theorem 1, so we have the following theorem. (17) is locally stable provided that

Theorem 4 The pest-free periodic solution of model
and is globally attractive if Based on the main results shown in Theorem 3, we can theoretically address the bifurcations for model (17), and we have the following main results.
Theorem 5 If T = T * , then the pest-free periodic solution (x p (t), y p (t)) can bifurcate to another periodic solution, which is supercritical if Kω ≤ 1.
Based on the above discussion, the dynamics of model (17) could be very complex when T > T * . To show this, we carry out a one-dimensional bifurcation analysis numerically. Figure 1 shows that if the impulsive period T exceeds some threshold levels, then both prey and predator populations can oscillate periodically with quite different amplitudes. As parameter T increases, the T-periodic solution loses its stability, then a period-doubling bifurcation occurs at T ≈ 8. As T further increases, the period-doubling bifurcations lead model (17) to chaos. After that, chaos disappears and period-halving bifurcations occur.   (17) with the parameters are identical to those in Fig. 1. The blue and green points are attracted to the attractors shown in Fig. 2

from top to bottom, respectively
When T = 14.5, a non-unique attractor appears, i.e. two different types of T-periodic solutions coexist (see Fig. 2). In addition, an interesting phenomenon is that 1-2-4-2-1 Tperiodic solutions occur periodically when T > 14.5, i.e. the dynamics of model (17) periodically recur as T increases, which means that only decreasing the period of predator releases could bring the pest population down to small amplitudes of periodic solutions. The above results reveal that the impulsive control period T plays a crucial role in dynamics and consequently in pest control.
The results shown in Fig. 3 reveal the basins of attraction for two alternative attractors, the blue and green areas are the two different types of T-periodic solutions, respectively. These results clarify that the final stable states of model (17) depend on their initial den- sities. This reveals the initial sensitivity of the model and the initial value dependence on the pest control. It also shows the importance of frequent monitoring of pest population in the process of pest control.
To further investigate the effect of the impulsive period T on the dynamical behaviour of system (17), consider the results shown in Fig. 4 which reveal that model (17) has more complex and interesting dynamic behaviour including period-doubling bifurcation, chaotic solutions, chaos crisis, period-halving bifurcations and periodic windows with increase of pulse period T. Meanwhile, bifurcation analyses also indicate that multiple attractors can coexist for a wide range of parameters. For example, Fig. 5 shows that a T-periodic solution coexists with a 6T-periodic solution when T = 21.1. If we choose different initial conditions (x 0 , y 0 ) = (1, 0.5) and (x 0 , y 0 ) = (2.5, 1.5), the stable attractors are the T-periodic solution and 6T periodic solution, respectively. If we fix all other parameter values as those in Fig. 4 and Fig. 6, it further shows that the solutions with different initial values will approach different amplitudes of attractors. In order to show how the nonlinear impulsive interventions affect the dynamics of model (17), we choose λ 1 as a bifurcation parameter and fix all other parameters as shown in Fig. 7(A) and (B). A period-doubling bifurcation occurs when λ 1 approaches 0.2, which means that a T-periodic solution disappears suddenly and a 2T-periodic solution occurs at this point. As λ 1 further increases, the 2T-periodic solution becomes unstable and the period-doubling bifurcations lead model (17) to chaotic dynamics. Following these, period-halving bifurcations result in various periodic solutions with different periods as λ 1 increases. When λ 1 > 5, there is a transition from periodicity to period-doubling bifurcation to a chaotic region with some narrow periodic windows. When λ 1 ≈ 13.1, the chaos bands disappear suddenly and a 2T-periodic solution occurs. With λ 1 in the range [13.1, 20], there are period-doubling bifurcation and period-halving bifurcations again. Meanwhile, we choose λ 2 as a bifurcation parameter as shown in Fig. 7(C) and (D). With λ 2 increasing, bifurcation analyses indicate that model (17) also exhibits very complex dynamics including period-doubling bifurcation, period-halving bifurcation, multi-stability and chaos. Numerical simulations show that the nonlinear impulsive interventions have a substantial influence on the dynamics of model (17) and cause the pest and natural enemy populations to oscillate periodically with different periods and amplitudes, which makes it more difficult for pest management decisions.
The results shown in Fig. 8 reveal how the maximal mortality rate due to the pesticide δ affects the dynamics of model (17). There exists a very interesting phenomenon of multiple attractors with a chaotic attractor and a periodic attractor always coexisting as parameter δ increases and in relation to changes of initial densities of pest and natural enemy populations. The red branches are generated from initial values (x 0 , y 0 ) = (2.7, 2.5), with the parameter δ increasing, the model presents transitions from 2T-periodic solutions to chaos by period-doubling bifurcation. However, the blue branches are generated from initial values (x 0 , y 0 ) = (2.5, 2.7), the model displays transitions from chaotic solutions to 2T-periodic solutions by period-halving bifurcation as the parameter δ increases. Note that period-doubling bifurcation and period-halving bifurcation occur in a dualism as parameter δ increases. In particular, if we choose initial values (x 0 , y 0 ) = (2.7, 2.5) and (x 0 , y 0 ) = (2.5, 2.7) respectively, 2T-periodic solutions via period-doubling bifurcation and 4T-periodic solutions via period-halving bifurcation coexist when δ = 0.31 (see Fig. 9). Further, Fig. 8 reveals that the model could approach chaos or stable attractors with different initial values when parameter δ is fixed, which further proves that initial densities of pest and natural enemy populations are crucial for pest control. Similarly, if we choose h as a bifurcation parameter, Fig. 10 also shows that different amplitudes and period of attractors always coexist as parameter h increases and in relation to changes of initial densities of pest and natural enemy populations. When we choose initial values (x 0 , y 0 ) = (2.7, 2.5), a series of period-halving bifurcation lead model to go from chaos to stable attractors as parameter h increases. If we choose initial values (x 0 , y 0 ) = (2.5, 2.7), model (17) undergoes a series of period-doubling bifurcation ending in a chaotic attractor as parameter h increases. We can see that period-doubling bifurcation and period-halving bifurcation occur in a dualism as parameter δ increases as well. Further, when fixing initial values and varying insecticide dosages (the maximal mortality rate δ or half saturation constant h), the model may have very complicated and interesting dynamic behaviour, such as the same types of periodic solutions that coexist, different types of periodic solutions coexisting, periodic solutions and chaotic attractors coexisting, which implies that the low initial densities of natural enemy will lead pest populations to chaos or to pest outbreaks or to small amplitude attractors with the parameters δ and h varying. Obviously, solutions with small amplitudes are biologically desirable. All these results have shown that nonlinear impulsive interventions have a substantial influence on the dynamics of model (17), which implies that the dynamics of the model are not only dependent on the selection of the initial densities of pest and natural enemy populations, but are also closely related to the system's parameters. Choosing appropriate parameters and initial densities can help us to design control strategies and to make appropriate management decisions.

Conclusion
The numbers of natural enemies to be released must be dependent on the densities of pest populations in the field. To show this, we focus on a generalised predator-prey model with a nonlinear impulsive control strategy. The existence and global stability of the pest-free periodic solution were addressed in detail, and the condition for permanence of the system was given. When T = T * , the trivial periodic solution loses its stability and a stable nontrivial periodic solution emerges via a supercritical bifurcation once a threshold condition is reached. To verify the main results obtained and to study the effects of parameter space on the threshold conditions, we further considered the Holling type II functional response function as an example.
One of the main purposes of this paper is to understand how nonlinear impulsive control affects the dynamics of the system, which will affect the success or failure of pest control strategies. To achieve these aims, we chose the nonlinear impulsive control parameters as bifurcation parameters, and our results indicated that nonlinear impulsive control makes the dynamical behaviour of systems change more dramatically and become more complicated. For example, bifurcation diagrams in Figs. 1 and 4 indicate that system (17) has more complex and interesting dynamic behaviour, including period-doubling bifurcation, chaotic solutions, chaos crisis, period-halving bifurcations, periodic windows and multistability (see Figs. 2 and 5). In particular, bifurcation diagrams reveal that the dynamic behaviour of the system depends on the initial densities of pest and natural enemy populations, which can help us to design strategies for controlling the pests by varying amounts and frequencies of insecticide applications.
Furthermore, the bifurcation diagram with respect to the maximal mortality δ and half saturation constant h reveals that the models with nonlinear impulsive interventions have very complicated and interesting dynamic behaviour, such as period-doubling bifurcation and period-halving bifurcation co-occurring, which implies that the model can be dramatically affected by small changes in the initial values or in the value of a relevant control parameter. In particular, there exists an interesting phenomenon whereby period-doubling bifurcation and period-halving bifurcation coexist as parameters δ and h increase, which implies that it is a very difficult task to control a pest's population as the pest population could have outbreaks in a quite complex way. So it is crucial for pest control to choose appropriate control parameters and initial values. Successful pest control depends not only on the initial densities of the pest and natural enemy populations, but also on the parameters related to the nonlinear impulsive interventions.
Our model with nonlinear impulsive interventions has more interesting dynamics and is more realistic than models with linear impulsive control, and nonlinear regulatory factors should be taken into account when considering IPM. Furthermore, it is believed that the methods developed in the generalised systems could be widely used to population dynamics in relation to glucose-insulin monitor model, tumour-immune dynamic model and epidemic model. However, it will also be more realistic if we consider delayed responses and residual effects of pesticides when these are also used as part of IPM on nonlinear impulsive control. Such future research will involve more complicated analyses of model dynamics.