Analyzing a generalized pest-natural enemy model with nonlinear impulsive control

1 Introduction

Over the past decade, controlling insect pests and other arthropods in agriculture became an increasing important issue. How to reduce losses due to insect pests becomes a great concern for entomologists and the society. To realize this purpose, a wide range of pest control strategies are available to farmers [1,2]. In particular, integrated pest management (IPM) has been proposed and designed which is a long term management strategy with aims to minimize economic, health and environmental risks by combining biological, chemical, cultural and physical tools [3,4,5,6].

With the development of the theory and application of impulsive differential equations [7,8], it is possible to depict the control strategies involved in IPM by establishing mathematical models. In particular, impulsive differential equations can accurately depict the dynamic process of spraying insecticides and releasing natural enemies [8,9,10,11,12,13,14]. For example, we can assume that the pesticide is applied at each fixed period, and a certain proportion of pests will be killed instantly after each spray. Similarly, the natural enemy is released simultaneously, where a constant amount of a natural enemy is administered at each fixed period. Note that the impulsive differential equations with fixed moments can provide a natural description for above assumptions, which can assist in the design and understanding of the ecological systems including the pests, their natural enemies, the surrounding environment and their inter-relationships. Moreover, theoretical analyses can provide valuable information about how to determine the optimal times or application frequencies of spraying pesticides, releasing natural enemies and infected pests [10,11,12,13,14,15].

Recently, many mathematical models concerning IPM have been developed and investigated [16,17,18,19,20,21,22,23]. In particular, the special prey-predator models with various functional response functions have been employed, and the main assumptions for those models are as follows: a proportion of pest population is killed after a pesticide is applied and simultaneously a natural enemy is released [24,25,26,27,28]. Based on those assumptions, the pest-natural enemy ecological systems with linear IPM measures have been extensively investigated, and almost all of those works focused on the existence and stability of pest free periodic solution. In particular, the threshold conditions under which the pest free periodic solution is locally or globally stable were provided, and this could help us to evaluate the effectiveness of pesticide and its application period on the successful pest control. Moreover, extensively numerical investigations revealed that those models could involve very complex dynamics including the coexistence of multiple attractors, chaotic solutions and period-doubling bifurcations.

The existence and stability of periodic solution for some generalized prey-predator model with linear or constant pulse actions studied so far [29,30,31,32] have provided some analytical techniques to deal with the generalized models. However, all of the previous works are basically assuming that the control strategy is linear or constant, which is not consistent with the actual situation. In fact, due to the resources limitation and the saturation effect of pesticide efficiency, the instant killing rate is a monotonic increasing function of pest population which should be a saturation function with a maximal killing rate. Similarly, the number of natural enemies is released depending on current pest density, which means that the larger the natural enemy, the fewer natural enemy is released and vice versa.

Therefore, the main purpose of this paper is to construct and investigate the dynamics of a quite generalized predator-prey model with nonlinear impulsive control due to resource limitation. By using Floquet theorem and qualitative techniques, it is proved that there exists a globally stable pest-free periodic solution under certain threshold conditions. By employing an operator theoretic approach which reduces the existence of the nontrivial periodic solutions to a fixed point and bifurcation problem, we show the existence and stability of positive periodic solution once the pest-free periodic solution loses its stability. In order to apply the main results, we choose the classical pest-natural enemy model with Holling type II functional response function and nonlinear impulsive control, then the exact thresholds for the local and global stability of pest free periodic solution are obtained. The results show that local stability does not imply the global stability which is confirmed by the bi-stability, and this is a novel result comparing with the model under the linear pulse perturbations [19,20,21,22,23,24,25].

2 The pest-natural enemy model with nonlinear pulse control

The generalized prey-predator model or pest-natural enemy model employed in the present paper is as follows:

where x(t), y(t) represent the densities of prey and predator populations, respectively. The function g(x) represents the intrinsic growth rate of the pest in the absence of natural enemy, p(x) denotes the predator response function, c is the efficiency rate, and D is the death rate of the predator population. In order to use our main results for a wide range of biological systems which have been investigated in the literature, we made the following assumptions related to the function p(x) and g(x). Let g(x) and p(x) be locally Lipschitz functions on R⁺ such that:

1. There exists a positive constant K > 0 such that g(x) > 0 for 0 ≤ x < K, g(K) = 0 and g(x) < 0 for K < x.

2. The functional response function satisfies p(0) = 0, p’(0) > 0 and p(x) > 0 for all x > 0.

3. The function xg(x)p(x)is upper bounded for all x > 0.

The first condition means that the pest population follows the density dependent growth in the absence of the natural enemy, and the second condition indicates that the functional response function is positive and monotonically increasing for small pest populations. The last one shows that the pest population can not increase infinitely once the biological control is introduced, and, on the other hand, if the density of pest population is too large then the biological control is impossible.

We assume that the IPM strategy is applied at every time point nT at which the natural enemies are released and pesticides are applied simultaneously, where T denotes the period of control actions and n ∈ 𝓝, which denotes the positive integer set. Moreover, the nonlinear saturation functions or density dependent functions are employed to depict the effects of resource limitation on the pest control, i.e. we choose

where δ ≥ 0 and h ≥ 0 represent the maximal fatality rate and the half saturation constant for the pest with δ < 1, λ ≥ 0 is the release amount of the natural enemy, and θ ≥ 0 denotes the shape parameter. We assume that the densities of both the pest and natural enemy populations are updated to (1−δx(t)x(t)+h)x(t)andy(t)+λ1+θy(t) at every discrete time point nT and n ∈ 𝓝, respectively.

Taking the control measures shown as the above and model (1) into account, one yields the following differential equation model with IPM strategies:

The positivity and boundedness of solutions of model (2) are useful for the coming analyses, and we have:

One of the main purposes of implementation IPM is to eradicate the pest population when the fixed moments are applied at every period T. To address this, 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 with nonlinear control, which can be analytically solved, and we have the following main result.

In particular, if θ = 0 then model (3) has a globally stable positive periodic solution

where y∗=λ1−exp⁡(−DT).

Therefore, we obtain the general expression of the unique pest-free periodic solution of model (2), i.e.

where y∗=−1+1+4λθexp⁡(−DT)(1−exp⁡(−DT))−12θexp⁡(−DT) when θ ≠ 0 or y∗=λ1−exp⁡(−DT) when θ = 0.

As mentioned before, one of the main purposes of IPM is to design suitable control measures such that the pest population dies out eventually, i.e. the pest free periodic solution (x_p(t), y_p(t)) is globally stable. Thus, the threshold conditions under which the pest free periodic solution is globally stable are crucial in this work. To do this, we first show the local stability, and consider the behavior of small amplitude perturbations of the solution.

where x̃(t) and ỹ(t) are small perturbations, then model (1) becomes

The impulsive effects on x̃ are unchanged because of x_p(t) = 0, so we have

The impulsive effects on ỹ are defined as follows:

The linear approximation of the deviation system of model (5) around the periodic solution (x_p(t), y_p(t)) is as follows:

In the following, we will show the sufficient condition for the global stability of pest-free periodic solution (x_p(t), y_p(t)) of model (2), and we have the following main results for model (2).

Comparing the formula of both R₀ and R₁ shown in equations (7) and (8) we can see that the conditions for the local and global stability of the pest free periodic solution are different, which depends on the relations between Ms=supxg(x)p(x) and g(0)f′(0). Note that limx→0xg(x)p(x)=g(0)f′(0) due to l’Hospital rule, and Ms=g(0)f′(0)⇔R0 = R₁ as x → 0. In general, the globally attractive condition is stronger than the local stability condition due to Ms≥g(0)f′(0)⇔R0≤R1. The question is whether the local stability implies the global stability of the pest free periodic solution, which will be discussed in more detail in the application section.

3 Threshold condition of bifurcation

For the existence of interior periodic solutions of model (2), we can investigate the bifurcation near the pest free periodic solution, i.e. (x_p(t), y_p(t)). To do this, for computation convenience we first exchange the variables x(t) and y(t), and denote u(t) = y(t) and v(t) = x(t), then system (2) becomes as follows:

We let Φ be the flow associated to the first two equations of (10), and the fundamental solution matrix of (10) is X(t) = Φ(t, X₀) with X₀ = X(0) = (u(0⁺), v(0⁺)) and Φ = (Φ₁, Φ₂). We employ the notations used in this section as those in [33], then we can define the mapping Θ₁, Θ₂: R² → R² as follows

and the mapping F₁, F₂: R² → R² by

Furthermore, we define Ψ : [0, +∞) × R² → R² by

Based on the above notations we can see that Ψ is determined by the values of solutions at impulsive points 0 and T, which is called as stroboscopic map of model (10) and T is the stroboscopic time snapshot. We know that X = (u, v) is a periodic solution of (10) with period T if and only if its initial value X₀ = X(0) is a fixed point for map Ψ(T, ⋅). Therefore, in order to establish the existence of nontrivial periodic solutions of (10), we should prove the existence of the nontrivial fixed points of Ψ.

For model (10), it follows from the discussion in the previous section that model (10) has a stable boundary T periodic solution, denoted by

where y∗=−1+1+4λθexp⁡(−DT)(1−exp⁡(−DT))−12θexp⁡(−DT) is a positive constant. In order to employ the analytical methods developed in [33, 34], we now consider the bifurcation of nontrivial periodic solutions near (u_s(t), 0) with initial value X(0) = (u_s(0⁺), 0).

In order to obtain a nontrivial periodic solution of period τ with initial value X(0), we have only to find the fixed point problem X = Ψ(τ, u). Denoting τ = T + τ̃, X = X₀ + X̃, and the fixed point problem X = Ψ(τ, u) equals to

Defining

so X₀ + X̃ is a fixed point of Ψ(T, ⋅) if N(τ̃, X̃) = 0.

According to the variational equations of the first two equations of (10), we have

which relates to the dynamics of the first two equations in (10). So we obtain that

with the condition D_X(Φ(0, X₀)) = I₂, which is the identity matrix in M₂(R). Thus, it follows from equation (10) that we have the particular form

with initial value D_X(Φ(0, X₀)) = I₂.

According to the initial value ∂Φ2(0,X0)∂u=0, we obtain the following

i.e. ∂Φ2(t,u0)∂u=0 for all t > 0. Further, we obtain

According to the initial condition D_X(Φ(0, X₀)) = I₂, we obtain that

for all 0 ≤ t ≤ T, where ρ(t) = exp ∫0t(g(0)−p′(0)us(v))dv.

We then compute the derivation of N according to (11), and observe the following matrix

Letting a′=a0′,b′=b0′,c′=c0′ and d′=d0′ when (τ̃, X̃) = (0, (0, 0)), at which ∂Θ1∂v=∂Θ2∂u=∂Φ2∂u=0, then we have

Thus, by simple calculations we have

and

Based on the above equations, we can see that D_XN(0, (0, 0)) is an upper triangular matrix with a0′>0. Thus, the necessary condition for the bifurcation of nontrivial solution is

which reduces to d0′=0. By simple calculation, we can see that d0′=0 is equivalent to R₀ = 1. Therefore, in the following we focus on d0′=0 and address the sufficient condition for the bifurcation of nontrivial solution.

Note that dim (ker(D_XN(0, (0, 0)))) = 1 and a basis of ker(D_XN(0, (0, 0))) is (−b0′a0′,1). Thus, the equation N(τ̃, X̃) = 0 is equivalent to

where E₀ = (1, 0), Y₀ = (−b0′a0′,1) and X̃ = aY₀ + zE₀ represents the direct summation decomposition of X̃ using the projections onto ker(D_XN(0, (0, 0))) (i.e. the central manifold) and Im(D_XN(0, (0, 0))) (i.e. the stable manifold).

Now, we define

and consequently we have