Abstract
We study the role of predation dynamics in oscillation of pest population in insect ecology. A two-dimensional pest control model (under the use of insecticides) with time delay in predation is considered in this paper. By the Hopf bifurcation theory, we prove the existence of the stable oscillation of the system. We also consider the economic viability of the control process. First we improve the Pontryagin maximum principle (PMP) where the delay in the system is sufficiently small and control function is linear, and then we apply the improved version of PMP to perform the optimal analysis of the pest control model as a special case.
1. Introduction
Long-term forecasts of pest pressure are central to the effective management of many agricultural and forest insect pests, because prediction of abundance and distribution (pest pressure) of a pest species is crucial for both strategic planning and tactical decision-making [1]. For example, timely forecasts would be useful for determining insecticide budgets. However, forecasting pest pressure is more problematic, because many factors influence the abundance of pests. For instance, abundance of predators interacting in the food web associated with the pest species or the effects of weather, rainfall, and temperature (e.g., [2]) may lead to a sustained oscillation in the pest population. There are a number of insects (mostly insects attacking forest trees) which exhibit fluctuating population density [3–5], although usually regular and occurring repeatedly between well-defined lower and upper boundaries. As these insects cause widespread economic damage, the reasons of their outbreaks have been a focus of intensive research.
There are various factors which cause the population changes in insect ecology such as density-dependent birth-death mechanism [5] or external factors like weather. For example, Stenseth et al. [6] found that the dynamics of Canadian lynx (Lynx canadensis) populations is mostly driven by weather condition. Since 1950s, ecologists have been proposing that “populations are limited either by extrinsic factors—such as weather, especially extremes of cold, drought, or rainfall—or by intrinsic factors—such as birth and death rates, or interactions with other species” [3, 5, 7, 8]. But the biological mechanisms that drive the oscillations are not yet well investigated.
1.1. A Case Study
Turchin and his group [9] designed a long-term experiment to test the hypothesis in favor of the existence of cycles in one forest insect, the southern pine beetle (SPB) Dendroctonus frontalis. They suspected that the mechanism behind this oscillation was beetle’s population interaction with its predators. Time-series analysis indicated that SPB fluctuations were driven primarily by endogenous (density-dependent) factors: approximately 80% of the variance in the rate of population change was explained by a joint action of current and lagged population densities. The evidence for second-order dynamics, that is, delayed density dependence [5], was strong, because regression of the rate of population change on lagged density was highly significant and it alone explained 55% of the variance [9].
However, they have demonstrated [10] and drawn three broad outcomes corresponding to the hypotheses that predation was (i) an exogenous, (ii) a first-order endogenous, or (iii) a second-order endogenous factor (see Figure 1 in Supplementary Material: SI_figures available online at http://dx.doi.org/10.1155/2013/653080). In the first case, there is no dynamical feedback between prey density and the predation impact. The average predator-induced mortality may be very high, and still predators would have no dynamical impact, simply reducing the intrinsic rate of population increase to a lower value. Fluctuations in predator-imposed mortality affect prey density in a stochastic manner, but they cannot drive a regular oscillation. In the second case, predators respond to changes in prey population without a significant lag time. Therefore the dynamical role of predators, in this case, was stabilizing rather than causing oscillations. Generalist predators may act in this manner, reducing the amplitude of oscillations or preventing diverging oscillations. Only in the third case, when acting in a delayed density-dependent manner, the predators are actually causing the oscillation. Thus their experimental study suggests that predation with time delay may be one of the fundamental causes of population oscillation in insect ecology.
1.2. Our Motivation
The experimental study by Turchin and his colleagues motivated us to formulate a mathematical framework to describe this oscillatory phenomenon in insect ecosystem. The goal of the paper is twofold. The first goal is to characterize the delay involved in predation that leads to oscillation in population and also what makes the oscillation stable! To mention, our study on describing oscillation for such kind of biological mechanism is not only to support Turchin’s experiments on SPB, but also important in other aspects of biological questions such as whether any control measure can stabilize the oscillation and if any, how it happens.
The second goal relates to the economic viability of the control process using optimization theory. We develop Pontryagin’s maximum principle (PMP) for delay differential equation where the delay is sufficiently small and underlined control function in the system is linear, and then we apply this theorem in our model. This is however the simplest condition necessary to develop PMP under delay. Similar things may be generalized with measurable functions under suitable restrictions, but this is outside the scope of current study. However, the second part on optimal control theory has its own importance independent of the first part.
2. Mathematical Model and Preliminaries
The mathematical formulation of pest control problem and its management (i.e., bionomic aspects of the pest control model under various control measure, like spraying of insecticide, release of additional predator or sterile males) have been taken up very recently by Bhattacharya and Karan [11], Bhattacharyya and Bhattacharya [12, 13], Ghosh et al. [14, 15]. But ecological aspect of the pest control problem has not been emphasized yet very much [16]. We consider a simple pest control model of a density dependent pest population under application of insecticide at the constant rate “” given by Figure 1 where and represent the density of the population of pest and its natural enemy, respectively. denotes the positive intrinsic birth rate of pest, is its carrying capacity, is the natural mortality of predators, and is the density dependence factor. The terms in the model “” and “” define the pest and predator mortality due to application of insecticide, where and denote proportional damage of species due to per unit application of pesticide. We consider the Holling-type II predational form for describing the grazing phenomenon [17, 18].
The initial condition for model (1) may be taken as any point in the region , where by we mean and denotes the set of real numbers.
Equilibria. The system has three equilibria, namely, the trivial equilibrium , axial equilibrium , and interior equilibrium , where and is given by the following equation: where , and .
When carrying capacity “” is high, the coefficient is negative and is positive, when . Hence, by Decarte’s rule of sign, (5) has either three or one positive root depending on values of other parameters. Moreover, if we notice the expression of the coefficients, it reveals that application of insecticide decreases the pest population in the system.
Now, the constant term , which shows that if we apply the insecticide of amount , pest will disappear from the system and the system reaches the trivial equilibrium . We will also show later that under this condition the trivial equilibrium becomes a stable one. Again, for the persistence of the predator in the system, it is necessary that and moreover Equation (6) indicates that the minimum number of pest is required for persistence of predator in the system. However, , and so it would be enough to keep the balance of the ecosystem, if we maintain
3. Model with Distributed Delay
We assume that the predation activity of predators involves some time lag. We also assume that this predation activity follows a gamma distribution. This form of distributed delay kernels are widely studied in many biological modelings ([19, 20] and references therein). Under the use of such kernels, our system (1) can be expressed as where , a nonnegative integer, is the order of the delay kernel and is a nonnegative real number. Both and are related with the mean time lag by . However we shall study the system for . With the transformation , the above system reduces to the following equivalent coupled differential equations:
The system (7) possesses the same equilibria as obtained for the system (1). It is easy to see that trivial equilibrium becomes asymptotically stable if , otherwise it becomes unstable for . Again, the axial equilibrium is locally asymptotically stable if . But this violates the positivity of in (4). Thus, existence of positive interior equilibrium implies that the axial equilibrium is an unstable saddle in character. The characteristic equation around is given by where Using Routh-Hurewitz criterion, we see that real parts of all roots are negative if holds, and thus is locally asymptotically stable. Hence, in this case, there is no possibility of exchange of stability.
Next we show by applying Benedixson-Dulac criteria [21] that the system (8) around has no periodic solutions.
Let , where is to be chosen suitably. Obviously, , whenever . We define
Then So, if we take , , and , then for all , , . However, these imply together . Therefore, there will be no periodic solution of the system (8).
Hence we find that the cyclic nature of the pest population in such forestry or agriculture ecosystem cannot be explained by this type of distribution of predation activity of natural enemies. It is also worth noticing that if the order of the delay kernel goes to infinity, while keeping the mean delay fixed, then the distributed delay can be viewed as a discrete delay [22]. Hence to explain the periodic nature of pest population, we shall assume that the process of predation activity follows a discrete time variation.
4. Model with Discrete Delay
Under this assumption that predation activity follows a discrete time variation, the system (1) takes the following form: where denotes the discrete time delay.
There are series of papers published by Chaplain and his colleagues on delayed predator-prey system. For example, Xu and Chaplain [23] studied the following delayed Gause-type predator-prey system without dominating instantaneous negative feedbacks mechanism: They investigated the uniform persistence of the system under some suitable parametric conditions. In another paper, Xu et al. [24] considered a ratio-dependent predator-prey model with distributed time delay. There is another recent paper by Lin and Yuan [25] on prey-predator system with distributed time delay, where the delay kernel was general gamma distribution: where the convolution is defined by . However, both of the studies do not consider the intraspecific competition of the predator species.
One special case of system (14) is the following model studied by Freedman and Ruan [26]: Zhao et al. [27] establish the existence of periodic solution of system by constructing an appropriate map and showing that the map has a nontrivial fixed point. Faria [28] considered the prey-predator system with delay in predation function: where and are constants. There is another most recent work on delayed system by Yan and Li [29] as a modification of the above system. This paper also discusses the global existence of periodic solution in the system. However, results in these works may be found as a special case of our system.
Another work that studied the effect of delay in prey-predator system is done by Chattopadhyay et al. [30]: They considered the phytoplankton-zooplankton system to understand the mechanism of harmful algal blooms in presence of toxic substance and incorporated the delay in the zooplankton mortality term due to liberation of toxic substance by phytoplankton.
The novelty of our paper lies in the fact that the results build on the current literature offering insight into population control strategies; and to best of our knowledge, all the studies mentioned above do not reflect this scenario. Our paper also gives important insight into the effects of control measures showing that a control agent can stabilize a population provided the delay in predator response function is sufficiently small. Thus, our paper discusses the delay problem in mathematical ecology in the light of control theory. This aspect indeed makes our paper important as well as different from earlier works on delayed problem.
However, for our system (14), as in the previous case, it has three equilibria, namely, a trivial equilibrium , an axial equilibrium , and a positive interior equilibrium . The trivial equilibrium is an unstable saddle point, and existence of makes the the axial equilibrium also an unstable saddle.
To investigate the nature of the system around , we perturb system (14) around and apply Nyquist theorem. A detailed discussion is given in Section 1 of Supplementary Material.
Bifurcation of the Solution. We consider the parameter time delay as bifurcation parameter of the system (14). Assuming as a solution of the characteristic equation of the perturbed system, we derive the following theorem.
Theorem 1. Suppose that (11) holds and . Then the real parts of the solutions of characteristic equation (see Supplementary Material) are negative for , where is the smallest value for which there is a solution of equation with real part zero. For , is unstable. Further increases through , and bifurcates into a small amplitude of periodic oscillation.
A detailed discussion on proof of the above theorem is given in Section 2 of Supplementary Material.
Hence we can observe that introduction of a delay into the predation activity of the natural enemy brings an oscillatory nature in the system (1). Thus our analysis supports the experimental findings of Turchin et al. [10].
5. Stability of the Bifurcation
In the previous section we obtained the condition that guarantees that the system (1) undergoes the Hopf bifurcation at the critical values given by Theorem 1. In this section we determine the formula that establishes the stability of bifurcating periodic orbits using the approach adapted in Hassard et al. [31]. The detailed calculation of the normal form and the specific condition on delay are derived in Supplementary Material. We see that stability of the bifurcation of the periodic solution depends on the system parameters. For example, the bifurcation is supercritical if and subcritical if . Also, if , the period of the solution increases with increase of , where where The terms and signs are defined in the Supplementary Material.
6. Optimization of Net Profit: The Second Goal
So far we have seen that under suitable threshold limits of model parameters and control parameters “,” the system has unique asymptotically stable equilibrium when delay “” is small, and for further increase of delay, the system enters into a Hopf bifurcation and gives a stable oscillation. But during the process considering spraying of insecticides as a control measure, there is an obvious question of incurring some cost and allied profit in the whole process. Thus the objective is to quantify the units to express the net profit during the given time of experiment. In other words, this precisely means that we are to construct an economic model out of the given dynamic model of pest population control. In this case the problem reduces to an optimal control problem. Our task is then to formulate an optimal policy when forms of control measure in the system are defined in a mathematical form and finally to find out the restrictions on the economic parameters of the model for making the policy viable. In fact, we perform the optimal analysis of the model (14) when delay is small enough. In this section we first improve the Pontryagin maximum principle (PMP) for linear control functions and delay is sufficiently small. However, we apply PMP to our model (14) of the system with constant control only to find out the optimal biomass of the species conveniently.
The most general form of the control problem involving small delay of discrete type is given by where is an -dimensional state vector and is an -dimensional control vector. Moreover, and . We assume that the vector valued function has continuous partial derivatives of all orders with respect to all its arguments. The set of all permissible control is denoted by . Moreover, it is assumed for definiteness that the control in questions is continuous from the left at their discontinuities, that is, We assume that the system (22) has a unique solution with initial condition and under the permissible control , for .
Hence under the above conditions, we now formulate the standard optimal control problem with delay of discrete type as follows:
We may mention here that the above optimal control problem is a relaxed optimal control problem, in the sense that we have not included the target in the decision problem. In ecological problem much work remains to be done on how we should specify the target in an optimal control problem which arises in ecology. In our case, we assume that the control measure, that is, insecticide is applied for certain period of time and the rate of control measure is linear function of time. The objective function consists of some desirable criterion at the end of the planning period plus a criterion involving the profit function which describes some desirable performance index during the planning period. Necessary conditions can be stated compactly in terms of Hamiltonian function. By definition Hamiltonian function is The functions are called costate variables. The set of necessary conditions that we have deduced for the above optimal control problem is as follows.
Theorem 2. Whenever delay in the system is sufficiently small and control function is linear in time , a necessary condition for an admissible set to be optimal is that there exists such that where denotes the transpose of the vector .
The details of the proof are given in Appendix of Supplementary Material.
Applying the necessary conditions (26) in Theorem 2, we will now continue the optimal analysis for our system. It may be mentioned once again that the control in questions of our particular system is constant. The profit function for our system is given by where is the projected profit per unit killing of pest and is the projected loss per unit killing of predator. The problem is then to optimize the integral over the control parameter , where . The detailed calculation of the optimal biomass of prey and predator and also optimal value of the control measure are derived in Supplementary Material.
7. Discussion and Conclusion
This paper considers the following two major issues:(a)the cyclic nature of pest population in an agroecosystem or a forestry ecosystem,(b)the improvement of Pontryagin’s maximum principle (PMP) for system of differential equations with small delay of discrete type and its application to the pest management problem.The experimental study of Turchin and his colleague [10] suggests that predation by natural enemy is one fundamental cause driving oscillation of pest population in the insect ecology. As the process behind the oscillation is still not clear, we have investigated the phenomenon through a simple model consisting of insect pest and its predator under two types of distribution of predation activity of the predator. We have observed that the cyclic nature of pest population, which is very common in agro- or forestry ecosystem, cannot be explained if the predation activity follows Holling-type II rule with gamma distribution. However, if it follows a delay of discrete type, we have observed that the system around the positive equilibrium enters a Hopf bifurcation and exhibits the cyclic nature for certain amount of time delay. To ascertain the local behavior, we have performed the stability analysis of bifurcating periodic solutions and have obtained the conditions for supercritical or subcritical bifurcations.
We have employed center manifold theory to show the stability of bifurcation of the oscillation. The approach of center manifold theory developed by Hassard and colleagues [31] shows that characteristic exponent determining the stability of the bifurcated periodic orbit is the same whether computed for original system or the system restricted on the associated center manifold. However, there is another existing approach showing that stability of bifurcating cycle is application of Poore’s condition [32]. The stability information of the bifurcating periodic orbits is contained in the following theorem of Poore, coupled with an algebraic expression, which completely reduces the determination of stability to an algebraic problem. In Poore’s theorem, the restricted original differential equation on center manifold reduced to a perturbed differential equation of a small parameter. Thus stability of the original periodic solution depends on characteristic multiplier of the variational matrix of perturbed system. This is a more sophisticated way of showing periodic stability compared to technique developed by Hassard and colleagues.
However, we have also performed a numerical experiment to substantiate the analytical findings with the following set of values (Table 1). Numerical solutions of equations were carried out using the modified fourth-order Runga-Kutta method. The parameter values are taken in such a way that it reflects the most realistic scenario in context of pest control problem.
For convenience of the simulation we have used the logarithmic scale. For those sets of values, we see that, under the condition in the Theorem 1, we have only one positive root and this implies . In fact, it is observed from Figure that for the system creates an stable limit cycle bifurcating from the positive asymptotically stable equilibrium. In fact, it is seen that, for , the system remains stable and approaches an asymptotically stable equilibrium (2.32034, 1.70924) (Figure 2). Further we obtained the value of in the eigenvector in (see Supplementary Material) and and in eigenvector . Moreover, it was found that , , , and . With these values of 's, we obtained the real part . This implies that and further . Thus we see that for this set of parameters, , that is, the bifurcation is supercritical, and the system exhibits stable limit cycle (Figures 3 and 4).
(a)
(b)
(a)
(b)
(c)
Further since , the period of the oscillations increases with increase of (Figures 5 and 6). The results indicate that the system has asymptotic equilibrium for and becomes unstable (by growing oscillations) for . The system exhibits a stable periodic solution for (Figures 7, 8, and 9). These observations indicate that there is threshold limit of the delay in the predation activity of the natural enemy, below which the system shows no excitability and above which the system enters into excitable range. These findings demonstrate the delay effect of predation activity and cyclic nature of pest population, which was suspected through experimental observations by Turchin and his colleagues [10]. We may mention here that delayed density-dependent predation is not only the single reason for oscillation of pest population in the insect ecology, but there are other factors like weather, especially extremes of cold, drought, or rainfall [33]. In spite of that our mathematical findings cannot be ignored, as it supports the experimental conclusion of Turchin and his colleagues. In fact, this will motivate the biologist to perform more explicit study in the laboratory in this direction.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Another important fact that we investigate in our study is control of this oscillation. A careful observation of the bifurcation analysis indicates that the coefficients which are responsible for rise and fall of pest population, its upper and lower boundaries, time period, even the value of bifurcation parameter (see Supplementary Material) are highly dependent on the control parameter . In fact, we perform some numerical experiments whenever , , and . It is seen that with these parameters estimation, and . This means that system will remain stable for (Figure ). Thus one can control the oscillation of pest population up to certain extent by applying control measures such as use of insecticide and release of additional predators. However, for larger delay, the system again enters into Hopf bifurcation and exhibits stable periodic solution (Figure ).
An allied problem with this control policy in pest control program is its economic viability. Some of the recent works on control theory in the respect of pest management problem could be found in the references like Bhattacharya and Karan [11, 34], and Bhattacharyya and Bhattacharya [35]. But all the earlier works involve application of Pontryagin’s maximum principle (PMP) on system of ordinary differential equations (under single control parameter or multicontrol parameters). In fact, control theory for ordinary differential equation is well known. But it is less developed for delay differential equations. The well-known PMP, which is an essential tool for solving the corresponding optimal control problem in ODE, has no exactly convenient form for delay differential equations. This part of the paper studies the corresponding problem with delay differential equations with small delay. It is better to mention here that system remains stable and solutions approach an asymptotically stable equilibrium for small delay. However, we have been able to obtain a set of necessary conditions for optimality whenever the delay in system is sufficiently small and control is linear function in time and then apply it suitably in our ecosystem model to get the optimal biomass under optimal control.
Finally, we would like to mention that dynamics of insect ecology is very complex one. There are indeed many factors, as we stated earlier, extrinsic factors such as weather, especially extremes of cold, drought, or rainfall; or by intrinsic factors such as birth and death rates, or interactions with other species may cause oscillation. Recent theoretical and empirical studies have started to combine the extrinsic and intrinsic hypotheses to explain the fundamental problem in insect ecology [36]. However, the present work on one hand gives some insight on this important ecological phenomenon and on the other hand throws some light on the control of the situation.
Acknowledgment
The authors are thankful to the anonymous reviewers for their comments and suggestions that greatly improved the paper.
Supplementary Materials
Supplementary Material: The supplementary materials contains figure from experimental study of Turchin and colleagues on population cycles (Turchin et al. 1999) and mathematical derivation and proofs of all underlined theorems described in the paper.