# Dynamics of Leslie-Gower Pest-Predator Model with Disease in Pest Including Pest-Harvesting and Optimal Implementation of Pesticide

Authors: Agus Suryanto; Isnani Darti
Journal: International Journal of Mathematics and Mathematical Sciences
 (2019)

---

## 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.


---

## Body

## 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 by(1a)dSdt=r1-S+IKS-βSI-c1SYk1+S-αS(1b)dIdt=βSI-cI(1c)dYdt=a2-c2Yk2+SYwith initial conditions (2)S0=S0>0,I0=I0>0,Y0=Y0>0.Here, S and I are the susceptible and infected pest, respectively, while Y 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 S-I epidemic model with bilinear incidence rate. Furthermore, the susceptible pest is continuously harvested with rate α.Table 1
Biological interpretation of parameters.
Parameter Biological interpretation r Intrinsic growth rate ofS K Carrying capacity ofS β Infection rate c Death rate ofI induced by the disease c1 Predation rate c2 Competition rate ofY a2 Intrinsic growth rate ofY k1 Constant of environment protection toS k2 Constant of environment protection toY α Harvesting rate ofS


## 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 ∀t>0.Proof.
Based on system (1a), (1b), and (1c), we directly get(3)St=S0exp∫0tr1-Sz+IzK-βIz-c1Yzk1+Sz-αdz>0,It=I0exp∫0tβSz-cdz>0,Yt=Y0exp∫0ta2-c2Yzk2+Szdz>0.
It is clear that every solution of initial value problem (1a), (1b), and (1c)-(2) will always be positiveLemma 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) ForS0=0, we directly get the result. For S0>0, the positivity solution gives St>0. Furthermore, from (1a) and (1b) we get(4)dSdt≤r1-SKS,dS+Idt≤r1-S+IKS,
which imply thatlimt→∞⁡supS≤K and limt→∞supS+I≤K.
(ii) The result is directly obtained forY0=0. Similarly, for Y0>0 we have Yt>0. From (1c) and previous result, we obtain(5)dYdt≤a2-c2Yk2+K+εY.
Therefore, we have thatlimt→∞supY≤k2+K+εa2/c2. Since ε>0 is arbitrary, we can conclude that limt→∞supY≤k2+Ka2/c2.
(iii) We setp(t)=S(t)+I(t)+Y(t) and calculate the derivative of pt with respect to t to get (6)dptdt=r1-S+IKS-c1SYk1+S-αS-cI+a2-c2Yk2+SY.
Now, for anyq<c we have that(7)dptdt+qpt≤r+qS+a2+qY.
For a givenε>0, we have t0>0 such that(8)dptdt+qpt≤σ+ε, for any t>t0 and σ=r+qK+a2+qk2+Ka2/c2. Applying Grönwall’s differential inequality, we get(9)pt≤pt0exp-qt-t0+σ+εq1-exp-qt-t0,
and fort→∞ and taking ε→0, we get(10)limt→∞suppt≤σq.
Hence, 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)
E1=(0,0,0), where all populations are extinct. The extinction of all populations equilibrium (E1) always exists.(ii)
E2=Kr-α/r,0,0, where only susceptible prey survives. The survival of susceptible prey equilibrium (E2) is feasible if r>α, i.e., when the intrinsic growth rate of prey is larger than the harvesting rate of susceptible prey.(iii)
E3=0,0,a2k2/c2, where prey is extinct and predator survives. The prey-free equilibrium (E3) is always feasible.(iv)
E4=S4∗,I4∗,0, where the predator is extinct with S4∗=c/β and I4∗=βKr-α-cr/ββK+r. The predator-free equilibrium (E4) is feasible if βKr-α>cr.(v)
E5=S5∗,0,Y5∗, where there is not infected prey. Here S5∗=-a2c1K+c2Kα-r+c2k1r+D/2c2r, D=a2c1K+c2Kα-r+c2k1r2-4c2rKa2k2c1+c2k1α-r, and Y5∗=a2S5∗+k2/c2. The free of infected prey equilibrium (E5) exists if r-α>c1a2k2/c2k1.(vi)
E∗=S∗,I∗,Y∗, where prey and predator coexist. Here, S∗=c/β,Y∗=a2βk2+c/βc2 and I∗=k1+S∗βKr-α-cr-c1βKY∗/c+βk1βK+r. The coexistence equilibrium is feasible if I∗>0.


### 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 E^=S^,I^,Y^ is(11)JE^=r1-2S^+I^K-βI^-c1k1Y^S^+k12-α-rK+βS^-c1S^S^+k12βI^βS^-c0c2Y^2S^+k220-2c2Y^S^+k2+a2. The direct evaluation of the Jacobian matrix at E1,E2,E3, and E4 is, respectively,(12)JE1=r-α000-c000a2,JE2=α-r-βK+rr-αrc1Kα-rk1r+Kr-α0βKr-α-crr000a2,JE3=r-α-c1a2k2c2k1000-c0a22c20-a2 and (13)JE4=r1-2c+βI4∗βK-α-βI4∗-cβK+rβK-cc1βk1+cβI4∗0000a2. All eigenvalues of JE1,JE2,JE3, and JE4 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 pointsE1,E2, and E4 are always unstable. If r<α+c1a2k2/c2k1, then E3 is locally asymptotically stable.The evaluation of Jacobian matrix (11) at E5 gives(14)JE5=r1-2S5∗K-α-c1k1Y5∗k1+S5∗2-rS5∗K-βS5∗-c1S5∗k1+S5∗0βS5∗-c0a22/c20-a2. The characteristics equation of JE5 is given by the following cubic equation: (15)λ3+A1λ2+A2λ+A3=0, where A1=-c1Y5∗/k1+S5∗2-r/KS5∗-βS5∗-c+a2,A2=S5∗c1Y5∗/k1+S5∗2-r/KβS5∗-c-a2S5∗c1Y5∗/k1+S5∗2-r/K+c1S5∗a22/c2k1+S5∗-a2βS5∗-c andA 3 = a 2 S 5 ∗ β S 5 ∗ - c c 1 Y 5 ∗ / k 1 + S 5 ∗ 2 - r / K - c 1 a 2 / c 2 k 1 + S 5 ∗. It is obvious that if c1Y5∗/k1+S5∗2<r/K and βS5∗<c, then A1>0,A2>0, A3>0 and also (16)A4=A1A2-A3=S5∗βS5∗-c-c1Y5∗k1+S5∗2-rK-βS5∗-c+a2c1Y5∗k1+S5∗2-rK-c1S5∗a22c2k1+S5∗c1Y5∗k1+S5∗2-rK-a2+a2βS5∗-c2-a22βS5∗-c>0.From the Routh-Hurwitz criterion, we can conclude the following result.Theorem 4.
Equilibrium pointE5 is asymptotically stable (locally) if c1Y5∗/k1+S5∗2<r/K and βS5∗<c.Using similar arguments, the stability of coexistence equilibrium can be obtained as follows.Theorem 5.
Ifc1Y∗/k1+S∗2<r/K, then the coexistence point (E∗) 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.
Ifr<α+c1a2k2/c2k1+K, then the pest-free equilibrium (E3) is globally asymptotically stable.Proof.
From Theorem3, we know that if r<α+c1a2k2/c2k1+K<α+c1a2k2/c2k1 then the equilibrium E3 is locally asymptotically stable. To show the global stability of E3, we have to prove that limt→∞St=0,limt→∞It=0, and limt→∞Yt=a2k2/c2. From (1c), we notice that (17)dYdt=a2-c2Yk2+SY≥a2-c2Yk2Y, which gives that limt→∞infY≥a2k2/c2. Based on this result, Lemmas 1 and 2(i), and (1a), we have that dS/dt≤rS-c1SY/k1+S-αS≤r-c1a2k2/c2k1+K-αS. Then, limt→∞St=0 if r<α+c1a2k2/c2k1+K. For any ε>0, there exists some sufficiently large t1>0 such that -ε<S(t)<ε for t≥t1. Using (1b), we obtain(18)-βε-cI≤dIdt≤βε-cI. By taking ε→0, we get the following comparison equation:(19)dMtdt=-cMt. Obviously M=0 is the only equilibrium of (19). Furthermore, M=0 is globally asymptotically stable equilibrium. Based on the comparison theorem, we get limt→∞It=0. Using similar argument, we can also show that limt→∞Yt=a2k2/c2.


### 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 r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5. Using these parameters, we show in Figure 1 the stability areas of E3,E5, and E∗ with respect to a2 and α. It is shown that, for a relatively small intrinsic predator growth rate (a2) and a relatively low harvesting rate of susceptible prey (α), the system will be convergent to the coexistence equilibrium (E∗). Increasing the value of a2 or α may cause the instability of E∗ and simultaneously may induce the stability of (E5). If the value of a2 or α is further increased, then equilibrium (E3) may be stable. Thus, there exists a transcritical bifurcation phenomenon driven by a2 or α. The dependence of stability properties of equilibrium can also be seen in Figure 2. It is seen that, for α=0, the coexistence point (E∗) is stable if a2≤0.825. Furthermore, the free of infected prey equilibrium (E5) is stable for 0.825<a2<1, while the prey free equilibrium (E3) is stable if a2≥1. If we take α=0.1, then the coexistence equilibrium, the free of infected prey, and the prey free equilibrium will be stable if a2≤0.7333,0.7333<a2<0.9 and a2≥0.9, respectively.Figure 1
The stability regions ofE3,E5 and E∗ with respect to a2 and α. The parameter values are: r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5.Figure 2
The stable equilibrium point as function ofa2 for ((a)-(c)) α=0 and ((d)-(f)) α=0.1. The parameter values are r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5.
(a) (b) (c) (d) (e) (f)Next we plot in Figure3 the stability regions of E3,E5, and E∗ using the same parameter values as before, except a2=0.25 and varying the death rate of prey induced by the disease (c). It is observed that, for relatively small rate of harvesting and c<0.2, 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 α=0 and α=0.8; see Figure 4. It can be seen from Figures 4(a)–4(c), for α=0, the coexistence point and the free of infected prey point will be stable if c<0.0651 and c≥0.0651, respectively. On the other hand, if we take α=0.8, then the prey free point will be stable for any c; see Figures 4(d)–4(f). Hence system (1a), (1b), and (1c) also exhibits transcritical bifurcation driven by c.Figure 3
The stability regions ofE3,E5 and E∗ with respect to c and α. The parameter values are r=1,K=1,β=0.1,a2=0.25,c1=0.5,c2=0.25,k1=1, and k2=0.5.Figure 4
The stable equilibrium point as function ofc for ((a)-(c)) α=0 and ((d)-(f)) α=0.8. The parameter values are r=1,K=1,β=0.1,a2=0.25,c1=0.5,c2=0.25,k1=1, and k2=0.5.
(a) (b) (c) (d) (e) (f)


## 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 ∀t>0.Proof.
Based on system (1a), (1b), and (1c), we directly get(3)St=S0exp∫0tr1-Sz+IzK-βIz-c1Yzk1+Sz-αdz>0,It=I0exp∫0tβSz-cdz>0,Yt=Y0exp∫0ta2-c2Yzk2+Szdz>0.
It is clear that every solution of initial value problem (1a), (1b), and (1c)-(2) will always be positiveLemma 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) ForS0=0, we directly get the result. For S0>0, the positivity solution gives St>0. Furthermore, from (1a) and (1b) we get(4)dSdt≤r1-SKS,dS+Idt≤r1-S+IKS,
which imply thatlimt→∞⁡supS≤K and limt→∞supS+I≤K.
(ii) The result is directly obtained forY0=0. Similarly, for Y0>0 we have Yt>0. From (1c) and previous result, we obtain(5)dYdt≤a2-c2Yk2+K+εY.
Therefore, we have thatlimt→∞supY≤k2+K+εa2/c2. Since ε>0 is arbitrary, we can conclude that limt→∞supY≤k2+Ka2/c2.
(iii) We setp(t)=S(t)+I(t)+Y(t) and calculate the derivative of pt with respect to t to get (6)dptdt=r1-S+IKS-c1SYk1+S-αS-cI+a2-c2Yk2+SY.
Now, for anyq<c we have that(7)dptdt+qpt≤r+qS+a2+qY.
For a givenε>0, we have t0>0 such that(8)dptdt+qpt≤σ+ε, for any t>t0 and σ=r+qK+a2+qk2+Ka2/c2. Applying Grönwall’s differential inequality, we get(9)pt≤pt0exp-qt-t0+σ+εq1-exp-qt-t0,
and fort→∞ and taking ε→0, we get(10)limt→∞suppt≤σq.
Hence, 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)
E1=(0,0,0), where all populations are extinct. The extinction of all populations equilibrium (E1) always exists.(ii)
E2=Kr-α/r,0,0, where only susceptible prey survives. The survival of susceptible prey equilibrium (E2) is feasible if r>α, i.e., when the intrinsic growth rate of prey is larger than the harvesting rate of susceptible prey.(iii)
E3=0,0,a2k2/c2, where prey is extinct and predator survives. The prey-free equilibrium (E3) is always feasible.(iv)
E4=S4∗,I4∗,0, where the predator is extinct with S4∗=c/β and I4∗=βKr-α-cr/ββK+r. The predator-free equilibrium (E4) is feasible if βKr-α>cr.(v)
E5=S5∗,0,Y5∗, where there is not infected prey. Here S5∗=-a2c1K+c2Kα-r+c2k1r+D/2c2r, D=a2c1K+c2Kα-r+c2k1r2-4c2rKa2k2c1+c2k1α-r, and Y5∗=a2S5∗+k2/c2. The free of infected prey equilibrium (E5) exists if r-α>c1a2k2/c2k1.(vi)
E∗=S∗,I∗,Y∗, where prey and predator coexist. Here, S∗=c/β,Y∗=a2βk2+c/βc2 and I∗=k1+S∗βKr-α-cr-c1βKY∗/c+βk1βK+r. The coexistence equilibrium is feasible if I∗>0.


## 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 E^=S^,I^,Y^ is(11)JE^=r1-2S^+I^K-βI^-c1k1Y^S^+k12-α-rK+βS^-c1S^S^+k12βI^βS^-c0c2Y^2S^+k220-2c2Y^S^+k2+a2. The direct evaluation of the Jacobian matrix at E1,E2,E3, and E4 is, respectively,(12)JE1=r-α000-c000a2,JE2=α-r-βK+rr-αrc1Kα-rk1r+Kr-α0βKr-α-crr000a2,JE3=r-α-c1a2k2c2k1000-c0a22c20-a2 and (13)JE4=r1-2c+βI4∗βK-α-βI4∗-cβK+rβK-cc1βk1+cβI4∗0000a2. All eigenvalues of JE1,JE2,JE3, and JE4 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 pointsE1,E2, and E4 are always unstable. If r<α+c1a2k2/c2k1, then E3 is locally asymptotically stable.The evaluation of Jacobian matrix (11) at E5 gives(14)JE5=r1-2S5∗K-α-c1k1Y5∗k1+S5∗2-rS5∗K-βS5∗-c1S5∗k1+S5∗0βS5∗-c0a22/c20-a2. The characteristics equation of JE5 is given by the following cubic equation: (15)λ3+A1λ2+A2λ+A3=0, where A1=-c1Y5∗/k1+S5∗2-r/KS5∗-βS5∗-c+a2,A2=S5∗c1Y5∗/k1+S5∗2-r/KβS5∗-c-a2S5∗c1Y5∗/k1+S5∗2-r/K+c1S5∗a22/c2k1+S5∗-a2βS5∗-c andA 3 = a 2 S 5 ∗ β S 5 ∗ - c c 1 Y 5 ∗ / k 1 + S 5 ∗ 2 - r / K - c 1 a 2 / c 2 k 1 + S 5 ∗. It is obvious that if c1Y5∗/k1+S5∗2<r/K and βS5∗<c, then A1>0,A2>0, A3>0 and also (16)A4=A1A2-A3=S5∗βS5∗-c-c1Y5∗k1+S5∗2-rK-βS5∗-c+a2c1Y5∗k1+S5∗2-rK-c1S5∗a22c2k1+S5∗c1Y5∗k1+S5∗2-rK-a2+a2βS5∗-c2-a22βS5∗-c>0.From the Routh-Hurwitz criterion, we can conclude the following result.Theorem 4.
Equilibrium pointE5 is asymptotically stable (locally) if c1Y5∗/k1+S5∗2<r/K and βS5∗<c.Using similar arguments, the stability of coexistence equilibrium can be obtained as follows.Theorem 5.
Ifc1Y∗/k1+S∗2<r/K, then the coexistence point (E∗) 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.
Ifr<α+c1a2k2/c2k1+K, then the pest-free equilibrium (E3) is globally asymptotically stable.Proof.
From Theorem3, we know that if r<α+c1a2k2/c2k1+K<α+c1a2k2/c2k1 then the equilibrium E3 is locally asymptotically stable. To show the global stability of E3, we have to prove that limt→∞St=0,limt→∞It=0, and limt→∞Yt=a2k2/c2. From (1c), we notice that (17)dYdt=a2-c2Yk2+SY≥a2-c2Yk2Y, which gives that limt→∞infY≥a2k2/c2. Based on this result, Lemmas 1 and 2(i), and (1a), we have that dS/dt≤rS-c1SY/k1+S-αS≤r-c1a2k2/c2k1+K-αS. Then, limt→∞St=0 if r<α+c1a2k2/c2k1+K. For any ε>0, there exists some sufficiently large t1>0 such that -ε<S(t)<ε for t≥t1. Using (1b), we obtain(18)-βε-cI≤dIdt≤βε-cI. By taking ε→0, we get the following comparison equation:(19)dMtdt=-cMt. Obviously M=0 is the only equilibrium of (19). Furthermore, M=0 is globally asymptotically stable equilibrium. Based on the comparison theorem, we get limt→∞It=0. Using similar argument, we can also show that limt→∞Yt=a2k2/c2.


## 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 r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5. Using these parameters, we show in Figure 1 the stability areas of E3,E5, and E∗ with respect to a2 and α. It is shown that, for a relatively small intrinsic predator growth rate (a2) and a relatively low harvesting rate of susceptible prey (α), the system will be convergent to the coexistence equilibrium (E∗). Increasing the value of a2 or α may cause the instability of E∗ and simultaneously may induce the stability of (E5). If the value of a2 or α is further increased, then equilibrium (E3) may be stable. Thus, there exists a transcritical bifurcation phenomenon driven by a2 or α. The dependence of stability properties of equilibrium can also be seen in Figure 2. It is seen that, for α=0, the coexistence point (E∗) is stable if a2≤0.825. Furthermore, the free of infected prey equilibrium (E5) is stable for 0.825<a2<1, while the prey free equilibrium (E3) is stable if a2≥1. If we take α=0.1, then the coexistence equilibrium, the free of infected prey, and the prey free equilibrium will be stable if a2≤0.7333,0.7333<a2<0.9 and a2≥0.9, respectively.Figure 1
The stability regions ofE3,E5 and E∗ with respect to a2 and α. The parameter values are: r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5.Figure 2
The stable equilibrium point as function ofa2 for ((a)-(c)) α=0 and ((d)-(f)) α=0.1. The parameter values are r=1,K=1,β=0.1,c=0.01,c1=0.5,c2=0.25,k1=1, and k2=0.5.
(a) (b) (c) (d) (e) (f)Next we plot in Figure3 the stability regions of E3,E5, and E∗ using the same parameter values as before, except a2=0.25 and varying the death rate of prey induced by the disease (c). It is observed that, for relatively small rate of harvesting and c<0.2, 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 α=0 and α=0.8; see Figure 4. It can be seen from Figures 4(a)–4(c), for α=0, the coexistence point and the free of infected prey point will be stable if c<0.0651 and c≥0.0651, respectively. On the other hand, if we take α=0.8, then the prey free point will be stable for any c; see Figures 4(d)–4(f). Hence system (1a), (1b), and (1c) also exhibits transcritical bifurcation driven by c.Figure 3
The stability regions ofE3,E5 and E∗ with respect to c and α. The parameter values are r=1,K=1,β=0.1,a2=0.25,c1=0.5,c2=0.25,k1=1, and k2=0.5.Figure 4
The stable equilibrium point as function ofc for ((a)-(c)) α=0 and ((d)-(f)) α=0.8. The parameter values are r=1,K=1,β=0.1,a2=0.25,c1=0.5,c2=0.25,k1=1, and k2=0.5.
(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 u reduces the growth rate of susceptible and infected prey, respectively, by the amount of ϵ1S and ϵ2I, where ϵ1<ϵ2. Therefore, system (1a), (1b), and (1c) is now modified as(20a)dSdt=r1-S+IKS-βSI-c1SYk1+S-αS-ϵ1S(20b)dIdt=βSI-cI-ϵ1I(20c)dYdt=a2-c2Yk2+SY 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: (21)J=∫0TAS+BI+Cu2dt subject to system of differential equations (20a), (20b), and (20c) and initial condition (2). Here, A and B are, respectively, the weight factors which correspond to susceptible and infected pest, while C 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(22)Ω=ut:0≤ut≤1,t∈0,T. We solve our optimal control problem using Pontryagin’s maximum principle. For that aim, we first define the Hamiltonian(23)H=AS+BI+Cu2+λ1dSdt+λ2dIdt+λ3dYdt where λit,i=1,2,3 are the adjoint (or costate) variables. The adjoint equations are(24a)dλ1dt=-∂H∂S=-A-λ1r1-2S+IK-βI-c1k1Yk1+S2-α-ϵ1u-λ2βI-λ3c2Yk2+S2(24b)dλ2dt=-∂H∂I=-B+λ1rSK+βS-λ2βS-c-ϵ2u(24c)dλ3dt=-∂H∂Y=λ1c1Sk1+S+λ32c2Yk2+S-a2, with transversality conditions (25)λ1T=λ2T=λ3T=0. Using the optimality condition, (26)∂H∂u=0, the optimal control is characterized by (27)u∗=λ1ϵ1S∗+λ2ϵ2I∗2C. Since the control variable has to be in Ω, the characterization is then given by(28)u∗=0ifλ1ϵ1S∗+λ2ϵ2I∗<0λ1ϵ1S∗+λ2ϵ2I∗2Cif0≤λ1ϵ1S∗+λ2ϵ2I∗≤2C1ifλ1ϵ1S∗+λ2ϵ2I∗>2C.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 ut. 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 asr=1,K=1,β=0.02,c=0.01,c1=0.3,c2=0.2,k1=1,k2=0.75,a2=0.5, and α=0.2. In this case, system (1a), (1b), and (1c) is convergent to the free of infected pest equilibrium E5=0.205,0.0,2.388 because c1Y5∗/k1+S5∗2=0.4932<r/K=1, 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 A=B=1. By keeping in mind that the cost to kill a single pest is quite low, we take C=0.025. Using those values of parameter we perform a simulation with initial values S0=0.3,I0=0.1 and Y0=0.2, where the pesticide control is assumed to be applied in 100 units of time; i.e., T=100. 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.Figure 5
The curves of state variables for system (20a), (20b), and (20c) without and with control: (a) susceptible pest, (b) infected pest, (c) predator, and (d) optimal control variable.
(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 (E1), the survival of susceptible prey equilibrium (E2), and the predator free equilibrium (E4) cannot be stable. The primary aim of controlling pest is to eliminate the pest from agriculture. Thus, the pest free equilibrium (E3) is the most important equilibrium due to the fact that if this equilibrium is stable then the pest is totally eradicated. E3 will be globally asymptotically stable if r<α+c1a2k2/c2k1+K. It is observed that the stability of E3 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.



---

*Source: 5079171-2019-06-18.xml*