Abstract
We propose a model which describes the interaction between pest and its natural predator. We assume that pest can be infected with diseases or pathogens such as bacteria, fungi, and viruses. The model is constructed by combining the Leslie-Gower model and S-I epidemic model. It is also considered the effects of pest harvesting. Harvesting in this case is intended to take a number of pests as one of the pest population control strategies. The proposed model will be analyzed dynamically to study its qualitative behaviour. The dynamical analysis includes the determination of all possible equilibrium points and their stability properties. Furthermore we also discuss the implementation of pesticide control where its optimal strategy is determined by Pontryagin’s maximum principle. To support our analytical studies, we perform some numerical simulations and their interpretation.
1. Introduction
Pest is organism that can be very damaging to agricultural production, which is a serious problem throughout the world [1]. Therefore, many researchers are interested in developing effective pest control methods. The classical method for eradicating pests is to spray chemical pesticides, which can quickly kill a significantly number of pests. However, the use of chemical pesticides is known to have side effects such as pests becomes resistant to pesticides or accumulation of chemical toxic residues in human food chain [2, 3]. To overcome the side effects of the use of chemical pesticides, pest control is carried out biologically which is expected to maintain the ecological environment and provide pest control solutions in a sustainable manner. An example of biological control is by releasing natural enemies to kill pest or maintain pest density to remain below the threshold of economic or ecological damage. An alternative biological control is by providing pathogens (bacteria, fungi, and viruses) to infect or to kill the pest [4–6]. Recently Meng et al. [7] also studied the combination effects of harvesting and pathogens in prey on the dynamics of prey and predator. To implement an integrated pest control we need to know what kind of natural enemy should be released, how many pest should be killed or trapped, and how much biopesticides (pathogens) or chemical pesticides should be sprayed. To regulate the pest control strategy, it is very useful to study the population changes of both pest and predator in a long time.
The interaction of pest with its natural predator with disease in pest is theoretically described by ecoepidemic model, which is a combination of the ecological model and the epidemiological model. The ecoepidemic model has been widely proposed and investigated [8–13]. Anderson and May [8] have investigated the interaction of pest (susceptible and infected) with its predator and observed that the infectious disease within animal and plant communities may destabilize the system. Recently, Ghosh et al. [10], Jana and Kar [11], and Kar et al. [12] proposed mathematical models which describe the interaction of pest and its predator where pest population is controlled by releasing infectious disease, natural enemy, and the use of chemical methods (pesticides). Here, the use of pesticide is optimized by Pontryagin’s maximum principle to minimize the implementation cost for the pesticide as well as its side effects, while effectively reducing the pest population. The models proposed in [10–12] are based on Lotka-Volterra model where the growth rate of predator is proportional to the rate of predation. If the predator is assumed to grow logistically where its carrying capacity depends on the size of prey, then we can apply a Leslie-Gower model. Ecoepidemic models based on the Leslie-Gower equation have been proposed by a number of researchers [14–19]. Zhou et al. [19] and Sharma and Samanta [16] considered an ecoepidemic model based on Leslie-Gower model with infectious disease in prey where the functional response follows the Holling type II. They considered that predator only consumes infected preys and the transmission of disease follows a bilinear incidence rate. The prey harvesting effect on model proposed in [19] has been investigated by Purnomo et al. [15]. Recently, Suryanto [17] has also modified model in [19] by assuming the saturated incidence rate. Moreover, Suryanto et al. [18] also studied the Leslie-Gower ecoepidemic model where predator only eats susceptible prey. In this paper, we develop model in [18] by including the effects of pest harvesting. The pest harvesting is intended to take a number of pests. Hence, the pest harvesting is considered as one of the pest population control strategies. The ecoepidemic model in this article is then given bywith initial conditions Here, and are the susceptible and infected pest, respectively, while represents the predator. All parameters in (1a), (1b), and (1c) are assumed to be positive and their biological interpretation are given in Table 1. Notice that model (1a), (1b), and (1c) is a combination of Leslie-Gower model and epidemic model with bilinear incidence rate. Furthermore, the susceptible pest is continuously harvested with rate .
2. Dynamical Properties
2.1. Positivity and Boundedness
It is obvious that the population is impossible to have negative value, while the boundedness of solutions can be understood as a natural limitation for growth as a consequence of limited resources. Positivity and boundedness of solutions of system (1a), (1b), and (1c) will be presented as follows.
Lemma 1. Every solution of system (1a), (1b), and (1c) with initial condition (2) will be positive .
Proof. Based on system (1a), (1b), and (1c), we directly getIt is clear that every solution of initial value problem (1a), (1b), and (1c)-(2) will always be positive
Lemma 2. (i) The number of preys is bounded above.
(ii) The number of predators is bounded above.
(iii) System (1a), (1b), and (1c) is uniformly bounded.
Proof. (i) For , we directly get the result. For , the positivity solution gives . Furthermore, from (1a) and (1b) we getwhich imply that and .
(ii) The result is directly obtained for . Similarly, for we have . From (1c) and previous result, we obtainTherefore, we have that . Since is arbitrary, we can conclude that .
(iii) We set and calculate the derivative of with respect to to get Now, for any we have thatFor a given , we have such that for any and . Applying Grönwall’s differential inequality, we getand for and taking , we getHence, system (1a), (1b), and (1c) is uniformly bounded.
2.2. Equilibria of the System
System (1a), (1b), and (1c) has six biologically feasible equilibria:(i), where all populations are extinct. The extinction of all populations equilibrium () always exists.(ii), where only susceptible prey survives. The survival of susceptible prey equilibrium () is feasible if , i.e., when the intrinsic growth rate of prey is larger than the harvesting rate of susceptible prey.(iii), where prey is extinct and predator survives. The prey-free equilibrium () is always feasible.(iv), where the predator is extinct with and . The predator-free equilibrium () is feasible if .(v), where there is not infected prey. Here , , and . The free of infected prey equilibrium () exists if (vi), where prey and predator coexist. Here, and . The coexistence equilibrium is feasible if .
2.3. Stability of Equilibria
The local stability of equilibrium point of system (1a), (1b), and (1c) is determined by the eigenvalues of Jacobian matrix of system (1a), (1b), and (1c). Here, the Jacobian matrix at an equilibrium point is The direct evaluation of the Jacobian matrix at , and is, respectively, and All eigenvalues of , and are very easy to be determined and analyzed. An equilibrium is asymptotically (locally) stable if the real parts of all eigenvalues of its Jacobian matrix are negative, and consequently we have the following result.
Theorem 3. Equilibrium points , and are always unstable. If , then is locally asymptotically stable.
The evaluation of Jacobian matrix (11) at gives The characteristics equation of is given by the following cubic equation: where and
. It is obvious that if and , then , and also
From the Routh-Hurwitz criterion, we can conclude the following result.
Theorem 4. Equilibrium point is asymptotically stable (locally) if and .
Using similar arguments, the stability of coexistence equilibrium can be obtained as follows.
Theorem 5. If , then the coexistence point () is asymptotically stable (locally).
2.4. Global Stability
The main goal of pest control is to eradicate pest. Mathematically, this can be achieved whenever the pest-free equilibrium is stable. The sufficient condition for the pest-free equilibrium to be globally asymptotically stable is given by the following theorem.
Theorem 6. If , then the pest-free equilibrium () is globally asymptotically stable.
Proof. From Theorem 3, we know that if then the equilibrium is locally asymptotically stable. To show the global stability of , we have to prove that , and . From (1c), we notice that which gives that . Based on this result, Lemmas 1 and 2(i), and (1a), we have that . Then, if . For any , there exists some sufficiently large such that for . Using (1b), we obtain By taking , we get the following comparison equation: Obviously is the only equilibrium of (19). Furthermore, is globally asymptotically stable equilibrium. Based on the comparison theorem, we get . Using similar argument, we can also show that .
2.5. Numerical Simulations and Discussion
To illustrate the dynamical behaviour of system (1a), (1b), and (1c), we perform some numerical simulations using hypothetical value of parameters. First, we set parameters , and . Using these parameters, we show in Figure 1 the stability areas of , and with respect to and . It is shown that, for a relatively small intrinsic predator growth rate () and a relatively low harvesting rate of susceptible prey (), the system will be convergent to the coexistence equilibrium (). Increasing the value of or may cause the instability of and simultaneously may induce the stability of (). If the value of or is further increased, then equilibrium () may be stable. Thus, there exists a transcritical bifurcation phenomenon driven by or . The dependence of stability properties of equilibrium can also be seen in Figure 2. It is seen that, for , the coexistence point () is stable if . Furthermore, the free of infected prey equilibrium () is stable for , while the prey free equilibrium () is stable if . If we take , then the coexistence equilibrium, the free of infected prey, and the prey free equilibrium will be stable if and , respectively.
(a)
(b)
(c)
(d)
(e)
(f)
Next we plot in Figure 3 the stability regions of , and using the same parameter values as before, except and varying the death rate of prey induced by the disease (). It is observed that, for relatively small rate of harvesting and , the prey-free equilibrium is always unstable. However, it can be seen that there is a minimal threshold of harvesting rate such that the prey-free state is always stable while other equilibrium points are unstable. This threshold should be satisfied if we want to eradicate the prey. In more detail, we also plot the dependence of equilibrium points stability for and ; see Figure 4. It can be seen from Figures 4(a)–4(c), for , the coexistence point and the free of infected prey point will be stable if and , respectively. On the other hand, if we take , then the prey free point will be stable for any ; see Figures 4(d)–4(f). Hence system (1a), (1b), and (1c) also exhibits transcritical bifurcation driven by .
(a)
(b)
(c)
(d)
(e)
(f)
3. Optimal Control Problem
We have discussed the dynamics of pest-predator model in which the pest population is controlled biologically by its natural enemy (predator) or the pest infection as well as mechanically by trapping (harvesting) the susceptible pest. In many cases, such control is not enough to eradicate the pest and therefore we need other types of controls. One of the well-known alternative controls is the use of pesticides. The use of pesticide controls to reduce the population of pest has been studied theoretically in [10–12]. Here, we will also implement the pesticide as a pest control in addition to the previous biological and mechanical controls. For this aim, we modify system (1a), (1b), and (1c) by assuming that the spraying of pesticide with rate reduces the growth rate of susceptible and infected prey, respectively, by the amount of and , where . Therefore, system (1a), (1b), and (1c) is now modified as subject to the initial conditions (2).
Our main goal is to reduce the number of pest (susceptible and infected), by applying pesticide. However, we have to consider the cost for applying pesticide which may be very high as well as its side effects which may be harmful for the environment or may cause the crops poisonous. To minimize the cost for the pesticide control and also its side effects, Kar et al. [12] and Joshi [20] suggested to minimize the square of the cost of applying the pesticide. Hence, our optimal control problem is to minimize the following objective functional: subject to system of differential equations (20a), (20b), and (20c) and initial condition (2). Here, and are, respectively, the weight factors which correspond to susceptible and infected pest, while is the weight factor related to the use of pesticide. The square of the control parameter is chosen to eliminate the side effects of the pesticide; see [12, 20]. Due to the practical limitations on the maximum rate of spraying pesticide, we consider the control set We solve our optimal control problem using Pontryagin’s maximum principle. For that aim, we first define the Hamiltonian where are the adjoint (or costate) variables. The adjoint equations are with transversality conditions Using the optimality condition, the optimal control is characterized by Since the control variable has to be in , the characterization is then given by
We now have our optimality system consisting of the state equations (20a), (20b), and (20c) with initial conditions (2), the adjoint (costate) equations (24a), (24b), and (24c) subject to the transversality conditions (25), and the optimal control (28). This optimality system will be solved numerically using the forward-backward sweep method [21]. We start by taking an initial guess of . Using the forward fourth-order Runge-Kutta method, we solve the state-equations with its initial conditions. Using the solutions of these state variables, we successively solve the adjoint equations with the transversality conditions by the backward fourth-order Runge-Kutta method. The control variable is then updated by substituting the values of the state and adjoint solutions into (28). This process is repeated until a convergent solution is achieved.
Since we are not considering a specific case of quantitative nature, we use the hypothetical parameter values as , and . In this case, system (1a), (1b), and (1c) is convergent to the free of infected pest equilibrium because , meaning that the pest will always be endemic. To reduce the number of pest we apply pesticide control. To minimize the total population of pest, we set the weight factors . By keeping in mind that the cost to kill a single pest is quite low, we take . Using those values of parameter we perform a simulation with initial values and , where the pesticide control is assumed to be applied in 100 units of time; i.e., . The number of susceptible pest, infected pest, and predator as function of time, both without and with pesticide control, are plotted in Figures 4(a)–4(c). It is observed that without control, both pest and predator survive in the environment. However, the optimal application of pesticide control can eradicate both susceptible and infected pest. Since the growth of predator also depends on the number of susceptible pest population, the population of predator is also reduced. In other words, the use of pesticide affects not only the pest population but also the population of predator. Furthermore, we also plot optimal pesticide control in Figure 5(d). From this figure we see that the use of pesticide is optimal when it is applied at maximum level for about 50 units of time and then reduce gradually the pesticide control until it stops at final time.
(a)
(b)
(c)
(d)
4. Conclusion
We have proposed an ecoepidemic pest-predator model describing the interaction of pest with its natural enemy where the pest may be infected by pathogens (bacteria, fungi, or viruses). We also consider the effect of pest harvesting. The presence of natural enemy and infectious pathogens is considered as biological control for the pest population, while pest harvesting can be considered as one of mechanical pest control. Our model is based on the Leslie-Gower equation in which the predator (natural enemy) grows logistically where its carrying capacity depends on the number of susceptible pest and other natural resources for the predator. Thus, although there is no pest, the predator still survives due to other natural resources. This explains why the extinction of both pest and predator equilibrium (), the survival of susceptible prey equilibrium (), and the predator free equilibrium () cannot be stable. The primary aim of controlling pest is to eliminate the pest from agriculture. Thus, the pest free equilibrium () is the most important equilibrium due to the fact that if this equilibrium is stable then the pest is totally eradicated. will be globally asymptotically stable if . It is observed that the stability of may be achieved if the rate of pest harvesting is suitably high.
In some cases, the combination of biological and mechanical control to eliminate the pest may be not very effective and therefore we need to include chemical control such as the use of pesticide. To consider the effect of such chemical control, we modify our model by including the effect of the use of pesticide. Then we formulate the optimal control problem to minimize the total number of pest, the cost of applying pesticide, and also the side effect of pesticide. Our numerical simulation shows that the optimal use of pesticide can eliminate totally the pest population.
This article presents numerical simulations with artificial parameters so that the results can only provide qualitative behaviour, rather than quantitative description. Nonetheless, the presented dynamical behaviour, the optimal control approach, and numerical simulations can provide description of the possible outcomes of the model. Appropriate mathematical model which is suitable for a specific case of natural phenomena can be obtained by estimating real parameters. Such estimation can be performed by fitting the real world data with the proposed model.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the Directorate of Research and Community Service, The Directorate General of Strengthening Research and Development, and the Ministry of Research, Technology and Higher Education (Brawijaya University), Indonesia, Contract no. 137/SP2H/LT/DRPM/III/2016 dated March 10, 2016, and Contract no. 460.18/UN10.C10/PN/2017 dated April 18, 2017.