# Ecoepidemiological Model and Analysis of Prey-Predator System

Authors: Abayneh Fentie Bezabih; Geremew Kenassa Edessa; Koya Purnachandra Rao
Journal: Journal of Applied Mathematics
 (2021)

---

## Abstract

In this paper, the prey-predator model of five compartments is constructed with treatment given to infected prey and infected predator. We took predation incidence rates as functional response type II, and disease transmission incidence rates follow simple kinetic mass action function. The positivity, boundedness, and existence of the solution of the model are established and checked. Equilibrium points of the models are identified, and local stability analyses of trivial equilibrium, axial equilibrium, and disease-free equilibrium points are performed with the method of variation matrix and the Routh-Hurwitz criterion. It is found that the trivial equilibrium pointEo is always unstable, and axial equilibrium point EA is locally asymptotically stable if βk−t1+d2<0,qp1k−d3s+k<0 and qp3k−t2+d4s+k<0 conditions hold true. Global stability analysis of an endemic equilibrium point of the model has been proven by considering the appropriate Lyapunov function. The basic reproduction number of infected prey and infected predators are obtained as R01=qp1−d32kβd3s2/qp1−d3qp1−d32kst1+d2+rsqp2kqp1−kd3−d3s and R02=qp1−d3qp3d3k+αrsqkqp1−kd3−d3s/qp1−d32t2+d4k, respectively. If the basic reproduction number is greater than one, then the disease will persist in the prey-predator system. If the basic reproduction number is one, then the disease is stable, and if the basic reproduction number is less than one, then the disease dies out from the prey-predator system. Finally, simulations are done with the help of DEDiscover software to clarify results.


---

## Body

## 1. Introduction

Mathematical modeling of prey-predator systems of interaction of species have a long history since the original remarkable work done by the Lotka-Volterra Model in the 1920s [1–4], and the SIR model compartment of systems of population is another vital area of research after the pioneering work of Kermack and Mckendrick [1–9]. Anderson and May were the first who combined these two modeling systems, while Chattopadhyay and Arino were the first who used the term “ecoepidemiology” for such models [2, 3, 6]. The dynamics of disease in prey-predator systems now become an interesting area of research due to the fact that prey-predator interaction is rich and complex in nature [4, 6, 10–13]. Several mathematical models have been proposed and studied on prey-predator systems [1–6, 8–12]. Many studies focused on the study of disease in prey only [1–3, 5, 6, 10, 12, 14–20], other researchers were interested in the study of disease within the predator population only [18, 21], and there are also some studies on diseases in both prey and predators [4, 8, 11]. In this paper, we proposed and studied infectious disease in both prey and predator interaction of species with treatment given to infected prey and infected predator.


## 2. Mathematical Model Formulation and Assumptions

In this paper, the prey-predator population is divided into five compartments. Let us denoteXt as the susceptible prey, Wt as the infected prey, Yt as the susceptible predator, Zt as the infected predator, and Ht as both infected prey and infected predator population under treatment. In the absence of infectious disease, the susceptible prey population grows logistically with intrinsic growth rate r and environmental carrying capacity k, and only susceptible prey can reproduce. In the presence of infectious disease, susceptible predators become infected predators when they come into contact with infected predators, susceptible prey become infected prey when they come into contact with infected prey, and the contact process is assumed to follow bilinear functional with convolution rate α,β, respectively. The predation functional response of predator towards the prey is assumed to follow a different Holling type II functional response form with p1,p2 respective predation coefficients of Xt,Wt due to susceptible predator, and p3,p4 respective predation coefficient of Xt,Yt due to infected predator. Suppose consumed prey converted into predator with efficiency q and also half-saturated constants. It is also assumed that infected prey Wt and infected predator Zt can only recover through treatment and are treated at a treatment rate of t1,t2, respectively. The prey-predator population Ht,Wt,Yt, and Zt suffers from infectious disease with death rate d1,d2,d3, and d4, respectively. Moreover, assume that all variables and parameters used in the model are nonnegative.According to the above assumption, we have the following model flow diagram.From the model flow diagram in Figure1, we have the following set of differential equations:
(1)dXdt=gX,W+r1H−βXW−p1f1X,Y−p3f3X,Z,(2)dWdt=βXW−t1W−d2W−p2f2W,Y−p4f4W,Z,(3)dYdt=qp1f1X,Y+qp2f2W,Y+r2H−αYZ−d3Y,(4)dZdt=qp3f3X,Z+qp4f4W,Z+αYZ−t2Z−d4Z,(5)dHdt=t1W+t2Z−d1H−r1H−r2H,with initial conditions X0≥0;W0≥0;Y0≥0;Z0≥0;H0≥0;andpi>0,i=1,2,3,4,and0<q≤1.Figure 1
Model diagram.Depending on the assumptions of per capita growth of functiongX,W for susceptible prey, and different type II functional responses fi,i=1,2,3,4, we have a more feasible model (8) emanating from model (3) as
(6)dXdt=rX1−X+Wk+r1H−βXW−p1XYs+X−p3XZs+X=fX,W,Y,Z,H,(7)dWdt=βXW−t1W−d2W−p2WYs+W−p4WZs+W=gX,W,Y,Z,H,(8)dYdt=qp1XYs+X+qp2WYs+W+r2H−αYZ−d3Y=hX,W,Y,Z,H,(9)dZdt=qp3XZs+X+qp4WZs+W+αYZ−t2Z−d4Z=iX,W,Y,Z,H,(10)dHdt=t1W+t2Z−d1H−r1H−r2H=jX,W,Y,Z,H,with initial conditions X0≥0;W0≥0;Y0≥0;Z0≥0;H0≥0;p1,p2,p3,p4>0,and0<q≤1, it is possible to arrange model (6)–(10) and express it in more compact forms with excellent patterns in
(11)dXdt=Xr1−X+Wk−βW+p1Y+p3Zs+X+r1H,dWdt=WβX−t1+d2+p2Y+p4Zs+W,dYdt=Yqp1Xs+X+p2Ws+W−αZ+d3+r2H,dZdt=Zqp3Xs+X+p4Ws+W+αY−t2+d4,dHdt=t1W+t2Z−d1+r1+r2H,with initial conditions X0≥0;W0≥0;Y0≥0;Z0≥0;H0≥0;p1,p2,p3,p4>0;and0<q≤1.


## 3. Mathematical Analysis of the Model

In this section, positivity, boundedness, and existence of the solution of the model are checked. This mathematical analysis of the model could be considered as the primary results.Theorem 1 (boundedness).
All solutions of model (8) are bounded in feasible region ℝ+5.Proof.
Each solutionXt, Wt, Yt, Zt, and Ht of the model is bounded if and only if total population N is bounded. Let total population of prey-predator N=X+W+Y+Z+H.
ForΛ>0 to be constant,
(12)dNdt+ΛN=dXdt+dWdt+dYdt+dzdt+dHdt+ΛN.By substituting all model Equations (8) into (12) and removing all negative terms, we have the following results: dN/dt+ΛN≤rX+qp1XY/s+X+qp2WY/s+W+qp4WZ/s+W+qp3XZ/s+X+ΛN=μ. Then, solving the differential inequality dN/dt+ΛN≤μ yields Nt≤μ/Λ1−e−Λt+N0e−Λt for t→∞,N→μ/Λ. We know that total prey-predator population is nonnegative, and hence, 0≤Nt≤μ/Λ. So we have invariant feasible region Ω=X,W,Y,Z,H∈ℝ+5:0≤Nt≤μ/Λ. This proves the theorem and the model is mathematically well posedTheorem 2 (positivity).
All solutions of model (6)–(10) are positive.Proof.
To prove Theorem2, we have to show that variables Xt,Wt,Yt,Zt,Ht of the models (8) are nonnegative ∀t≥0.
(i)
Positivity ofXt: from the susceptible prey model in (8), dX/dt=rX1−X+W/k+r1H−βXW−p1XY/s+X−p3XZ/s+X without loss of generality. After removing all the positive terms from the right-hand side of the differential equation, we have the following differential inequality: dX/dt≥−rX2+rXW/k+βXW+P1XY+P3XZ/S+X divide both sides by negative yields −dX/dt≤RX2+RXW/k+βXW+p1XY+p3XZ/S+X. But it is also clear that the following inequality holds: rX2+rXW/k+βXW+P1XY+P3XZ/S+X≤rX2+rXW+βXW+p1XY+p3XZ=XrX+rW+βW+p1Y+p3Z. Assume that rW+βW+p1Y+p3Z=C, then the differential inequality is reduced to −dX/dt≤XrX+C. This inequality can be arranged for integration by partial fraction as ∫1/XrX+CdX≥∫−dt, integrating the integral inequality ∫1/C/X+−r/C/rX+CdX≥−∫dt will give us 1/ClnX−1/ClnrX+C≥−t+Q, where Q is integration constant. Using rules of logarithm, the inequality can be written as lnX/rX+C≥−Ct+CQ. Finally, solving for X will give us Xt≥ACe−Ct/1−rAe−Ct, for A=eCQ. Therefore, Xt>0 for 1−rAe−Ct>0. That is, Xt is nonnegative for t>1/ClnrA.(ii)
Positivity ofWt: from the infected prey model in (7), dW/dt=βXW−t1W−d2W−p2WY/s+W−p4WZ/s+W, without loss of original generality, after removing the positive term (βXW). We obtain the following differential inequality: dW/dt≥−t1W+d2W+p2WY/s+W+p4WZ/s+W if and only if −dW/dt≤t1W+d2W+p2WY/s+W+p4WZ/s+W. But it is clear that t1W+d2W+p2WY/s+W+p4WZ/s+W≤t1W+d2W+p2WY+p4WZ=t1+d2+p2Y+p4ZW holds true. Now assume that t1+d2+p2Y+p4Z=C. Then, we have −dW/dt≤CWt. Now applying integration yield lnW≥−Ct+Q, where Q is the integration constant, then solving for the variable Wt gives the equation Wt≥e−Ct+Q which is the exponential function and positive at all times. Hence, Wt is positive.(iii)
Positivity ofYt: from the susceptible predator model in (8), dY/dt=qp1XY/s+X+qp2WY/s+W+r2H−αYZ−d3Y, without loss of original generality, after removing all positive terms qp1XY/s+X+qp2WY/s+W+r2H, we obtain differential equation dy/dt≥−αz+d3y. Then, applying integration by separation of variable method results, lny≥−αz+d3t+Q, where Q is integration constant, and solving for variable Yt,we obtain the solution y≥e−αz+d3t+Q. Therefore, yt≥e−αz+d3t+Q is a positive exponential function; hence, yt is positive.(iv)
Positivity ofZt: from the infected predator model in (9), dZ/dt=qp2WZ/s+W+qp3XZ/s+W+αYZ−t2Z−d4Z after removing all positive terms qp2WZ/s+W+qp3XZ/s+W, we obtain the differential inequality dz/dt≥−t2+d4z. Applying integration by separation of variable method yields lnz≥−t2+d4t+Q where Q is integration constant by separation of variable method. Then, solving for Z will result in zt≥e−t2+d4t+Q which is the exponential function that is positive at all time. Hence, Zt is positive.
Alternative verification: in mode (6)–(10), dZ/dt=qp2WZ/s+W+qp3XZ/s+W+αYZ−t2Z−d4Z Can be written as dZ/dt=qp2W/s+W+qp3X/s+W+αY−t2−d4Z if only if 1/ZdZ=qp2W/s+W+qp3X/s+W+αY−t2−d4dt. This equation can be arranged as ∫1/ZdZ=∫qp2W/s+W+qp3X/s+W+αY−t2−d4dt, which after computing the integration yields, lnZ=qp2W/s+W+qp3X/s+W+αY−t2−d4t+Q, where Q is the integration constant. Thus, Zt=eqp2W/s+W+qp3X/s+W+αY−t2−d4t+Q is the exponential function which is positive, and hence, Zt is positive.
(v)
Positivity ofHt: from infected prey and infected predator population under treatment model in (8) dH/dt=t1W+t2Z−d1H−r1H−r2H, without loss of generality, after removing all positive terms, we have dH/dt≥−d1+r1+r2H if and only if dH/H≥−d1+r1+r2dt integrating results, lnH≥−d1+r1+r2t+Q. Then, solving the variable H yields H≥e−d1+r1+r2t+Q which is an exponential function which is positive at all times. Therefore, Ht>0, and hence, Ht is positive. Thus, variables Xt,Wt,Yt,Zt, and Ht are all positive quantities and remain in ℝ+5 for all t.Theorem 3 (existence).
All solutions of model (8) together with the initial conditions X0>0,W0≥0,Y0≥0,Z0≥0,andH0≥0 exist in ℝ+5, i.e., the model variables Xt,Wt,Yt,Zt, and Ht exist for all t and remain in ℝ+5.Proof.
Let the system of differential equation (8) be given as
(13)f1=rX1−X+Wk+r1H−βXW−p1XYs+X−p3XZs+X,f2=βXW−t1W−d2W−p2WYs+W−p4WZs+W,f3=qp1XYs+X+qp2WYs+W+r2H−αYZ−d3Y,f4=qp3XZs+X+qp4WZs+W+αYZ−t2Z−d4Z,f5=t1W+t2Z−d1H−r1H−r2H.
According to the Derrick and Groosman theorem, letΩ denote the region Ω=X,W,Y,Z,H∈ℝ+5;N≤μ/λ. Then, model (8) have a unique solution if ∂fi/∂xj,i,j=1,2,3,4,5 are continuous and bounded in Ω. Here, x1=X,x2=W,x3=Y,x4=Z,andx5=H. The continuity and the boundedness can be shown as follows.
Thus, all the partial derivatives in Table1, ∂fi/∂xj,i,j=1,2,3,4,5 exist, continuous, and bounded in a region Ω for all positive values of model variables in Table 2 and model parameters in Table 3. Hence, by the Derrick and Groosman theorem, a solution for the model (6)–(8) exists and is unique.Table 1
Partial derivatives.
Forf1:∂f1∂X=r−2rxk−rWk−βW−sp1Y+p3Zs+X2<∞,∂f1∂W=−rXk−βX<∞,∂f1∂Y=−p1Xs+X<∞,∂f1∂Z=−p3Xs+X<∞,∂f1∂H=r1<∞,Forf3:∂f3∂X=sqp1Ys+X2<∞,∂f3∂W=sqp2Ys+X2<∞,∂f3∂Y=qp1Xs+X+qp2Ws+W−αZ−d3<∞,∂f3∂Z=−αy−d3<∞,∂f3∂H=r2<∞,Forf2:∂f2∂X=βW<∞,∂f2∂W=βX−t1−d2−sp2Y+p4Zs+W2<∞,∂f2∂Y=−p2Ws+W<∞,∂f2∂Z=−p4Ws+w<∞,∂f2∂H=0<∞,Forf4:∂f4∂X=sqp3Zs+X2<∞,∂f4∂W=sqp4Zs+W2<∞,∂f4∂Y=αZ<∞,∂f4∂Z=qp3Xs+X+qp4Ws+W−αY−t2−d4<∞,∂f4∂H=0<∞.Forf5:∂f5∂X=0<∞,∂f5∂W=t1<∞,∂f5∂Y=0<∞,∂f5∂Z=t2<∞,∂f5∂H=−d1−r1−r2<∞.Table 2
Notations and description of variables.
VariablesDescriptionsXtPopulation size of susceptible preyWtPopulation size of infected preyYtPopulation size of susceptible predatorZtPopulation size of infected predatorHtPopulation size of infected population under treatmentTable 3
Notations and description of parameters.
ParametersDescription of parametersr,kIntrinsic growth rate and carrying capacity of susceptible preyα,βDisease transmission rates of prey and predatort1,t2Treatment rate of infected prey and infected predatorr1,r2Recovery rate of infected prey and infected predatorpi,i=1,2,3,4Predation coefficientsfi,i=1,2,3,4Functional responsedi,i=1,2,3,4Death ratesq,sEfficiency of predation and half-saturation constant


## 4. Stability Analysis

Stability analysis in the absence of predators in the model, that is, whenyt and Zt are zero, model (8) can be written as
(14)dXdt=rX1−X+Wk+r1H−βXW=fX,W,H,dWdt=βXW−t1W−d2W=gX,W,H,dHdt=t1W−d1H−r1H=hX,W,H.This system has trivialEo0,0,0, axial EAk,0,0, and positive EoX,W,H equilibrium points where
(15)X=k−kβr−βd2+t1+kβr1t1rd1+r1d2+t1,W=βd2+t1,H=βd1+r1d2+t1,with Jacobian matrix(16)JX,W,H=r−2rXk−rWk−βW−rXk−βXr1βWβX−t1−d200t1−d1−r1.Theorem 4.
The trivial equilibriumEo is a saddle point with unstable manifold in X-direction and stable manifold in the WY-plane.Proof.
The Jacobian matrix atEo is given by
(17)JEo=r0r10−t1−d200t1−d1−r1,to compute eigenvalues compute the detJEo−λI3=0,
(18)r−λ0r10−t1−d2−λ00t1−d1−r1−λ=0.
Then,r−λ−t1−d2−λ−d1−r1−λ=0 is the characteristic polynomial.
Thus, eigenvalues areλ1=r>0,λ2=−t1−d2<0,λ3=−d1−r1<0 which is a saddle point with unstable manifold in the X-direction and stable manifold in the WY-plane.Theorem 5.
The axial equilibriumEA is a saddle point if βk−t1−d2>0 and unstable manifold in X-direction if βk−t1−d2<0, then EA is stable.Proof.
The Jacobian matrix atEA is given by
(19)JEA=−r−r−βkr10βk−t1−d200t1−d1−r1,
to compute eigenvalues compute thedetJEo−λI3=0,
(20)−r−λ−r−βkr10βk−t1−d2−λ00t1−d1−r1−λ=0.
Then,−r−λβk−t1−d2−λ−d1−r1−λ=0 is characteristic polynomial.Thus,λ1=−r<0,λ2=βk−t1−d2,λ3=−d1−r1<0; hence, the axial equilibrium point is the saddle point if βk−t1−d2>0 and stable if βk−t1−d2<0.Stability analysis in the absence of infectious disease in system (8), that is, when there is no disease, Wt,Zt, and Ht are all zero and model (8) becomes
(21)dXdt=rX1−X+Wk−p1XYs+X=fX,Y,dYdt=qp1XYs+X−d3Y=gX,Y.This system contains trivialEo0,0, axial EAk,0, and positive EoX,Y equilibrium points, where
(22)X=rk−s+rk2r+2krs+rXs2−4kp12r,Y=qp1d3s+rk−s+rk2r+2krs+rs2−4kp1/2r,and Jacobian matrix is given by
(23)JX,Y=r−2rXk−p1YSs+X2qp1Xs+kqp1YSs+X2qp1Xs+X−d3.Theorem 6.
The trivial equilibriumEo is a saddle point with unstable manifold in the X-direction and stable manifold in the Y-direction.Proof.
The Jacobian matrix atEo is given by
(24)JEo=r00−d3,
hence, eigenvalues areλ1=r>0,λ2=−d3<0 which is a saddle point.Theorem 7.
The axial equilibriumEA is stable if qp1k/s+k−d3<0, otherwise unstable.Proof.
The Jacobian matrix atEA is given by
(25)JEA=−r−p1ks+k0−p1ks+k−d3.
To find eigenvalues, computedetJEA−λI3=0,
(26)−r−λ−p1ks+k0qp1ks+k−d3−λ=0.
Eigen values areλ1=−r<0,λ2=qp1k/s+k−d3.
Thus,EA is stable if qp1k/s+k−d3<0 and otherwise unstable.Theorem 8.
The positive equilibriumE is stable if r−2rX/k−p1YS/s+X2+qp1k/s+k−d3>0 and r−2rX/k−p1YS/s+X2qp1k/s+k−d3+qp12XYS/s+X3>0.Proof.
The Jacobian matrix atE is given by
(27)JX,Y=r−2rXk−p1YSs+X2qp1Xs+kqp1YSs+X2qp1Xs+X−d3.
Then, computedetJE−λI3=0,
(28)r−2rXk−p1YSs+X2−λqp1Xs+kqp1YSs+X2qp1Xs+X−d3−λ=0.
Then,(29)r−2rXk−p1YSs+X2⏟a−λqp1Xs+X−d3⏟b−λ+sqp12XYs+X3⏟c=0,
is a characteristic polynomial.
Using the Routh-Hurwitz criterion, the quadratic polynomial is stable ifa+b>0,ab+c>0, otherwise unstable.Equilibrium points in model (8) are steady state points of the form X,W,Y,Z,H of model ((8)) that satisfies dX/dt=dW/dt=dY/dt=dZ/dt=dH/dt=0, provided that each variable is nonnegative. In model ((8)), five steady state points are identified and listed here: trivial steady state Eo0,0,0,0,0, axial steady state EAk,0,0,0,0, disease-free steady state E¯X¯,0,Y¯,0,0, and endemic steady state E∗X∗,W∗,Y∗,Z∗,H∗. Computation of disease free and endemic equilibrium points is presented as follows.Disease-free equilibrium points (DFEP) of model (8) are steady state solutions when there is no infectious disease in the population. In the absence of infectious disease in the prey-predator system, the variables Wt=Zt=Ht=0 and dX/dt=dW/dt=dY/dt=dZ/dt=dH/dt=0. Then, model (8) become
(30)rX¯1−X¯k−p1X¯Y¯s+X=0,qp1X¯Y¯s+X¯−d3Y¯=0.Thus, solvingX¯ and Y¯, it is found that X¯=d3s/qp1−d3andY¯=rsqkqp1−kd3−d3s/qp1−d32k, and hence, disease-free equilibrium point (DFEP) of model (8) is given by
(31)E¯=X¯,0,Y¯,0,0=d3sqp1−d3,0,rsqkqp1−kd3−d3sqp1−d32k,0,0.The endemic equilibrium point (EEP) is positive equilibrium pointE∗X∗,W∗,Y∗,Z∗,H∗ obtained by solving model equation (8) as dX/dt=dW/dt=dY/dt=dZ/dt=dH/dt=0 for which all variables are nonzero
(32)rX∗1−X∗+W∗k+r1H∗−βX∗W∗−p1X∗Y∗+p3X∗Z∗s+X∗=0,βX∗W∗−t1W∗−d2W∗−p2W∗Y∗+p4W∗Z∗s+W∗=0,qp1X∗Y∗s+X∗+qp2W∗Y∗s+W∗+r2H∗−αY∗Z∗−d3Y∗=0,qp3X∗Z∗s+X∗+qp4W∗Z∗s+W∗+αY∗Z∗−t2+d4Z∗=0,t1W∗+t2Z∗−d1+r1+r2H∗=0.Then, solving for the variablesX∗,W∗,Y∗,Z∗,andH∗, the endemic equilibrium points of the model exists and a simplified result is obtained
(33)X∗=−rβd1+2kr−rs−2kβd1d2−rβr1+2krd2r1,W∗=β−sd2+t1−4sβ−p2−p4t1+d2+st1+d2−β22t1+d2,Y∗=1d3−α−rr2r1+rβr2r1+rβr22kr1t1+d2−rsd2r22kr1t1+d2−rsr2t12kr1t1+d2,Z∗=1t2+d4α+qp4s+β−t1+d2+4sβ−p2−p4t1+d2+st1+d2−β2/2t1+d2,H∗=1r1β−r+rβ−rsd2−rst12kt1+d2.To study the stability analysis of equilibrium points of model (8), it is better to linearize model (8) using variation matrix. Let us represent the model (8) as system of parametrical equations as follows:
(34)dXdt=fX,W,Y,Z,H=rX1−X+Wk+r1H−βXW−p1XYs+X−p3XZs+X,(35)dWdt=gX,W,Y,Z,H=βXW−t1W−d2W−p2WYs+W−p4WZs+W,(36)dYdt=hX,W,Y,Z,H=qp1XYs+X+qp2WYs+W+r2H−αYZ−d3Y,(37)dZdt=iX,W,Y,Z,H=qp3XZs+X+qp4WZs+W+αYZ−t2Z−d4Z,(38)dHdt=jX,W,Y,Z,H=t1W+t2Z−d1H−r1H−r2H.Then, the variation matrix of these functions (36) is given by
(39)VX,W,Y,Z,H=fXfWfYfZfHgXgWgYgZgHhXhWhYhZhHiXiWiYiZiHjXjWjYjZjH,where each element of the matrix represents partial derivatives of functions (36) with respect to model variables and computations of each element of the variation matrix are given as
(40)VX,W,Y,Z,H=fX−rXk−βX−p1Xs+X−p3Xs+Xr1βWgW−p2Ws+W−p4Ws+w0sqp1Ys+X2sqp2Ys+X2hY−αy−d3r2sqp3Zs+X2sqp4Zs+W2αZiz00t10t2−d1−r1−r2,where, fX=r−2rx/k−rW/k−βW−sp1Y+p3Z/s+X2,gW=βX−t1−d2−sp2Y+p4Z/s+W2,hY=qp1X/s+X+qp2W/s+W−αZ−d3,iz=qp3X/s+X+qp4W/s+W−αY−t2−d4.Theorem 9 (TEP).
Trivial equilibrium pointEo0,0,0,0,0 exists and always locally asymptotically unstable.Proof.
Consider the variation matrix (40) at Eo.
(41)VEo=r000r10−t1−d200000−d30r2000−t2−d400t10t2−d1−r1−r2.
Eigenvalue of variation matrix can be computed from the characteristic polynomialdetVEo−λI5=0(42)r−λ000r10−t1−d2−λ00000−d3−λ0r2000−t2−d4−λ00a10a2−d1−r1−r2−λ=0,
⟹r−λ−t1−d2−λ−d3−λ−t2−d4−λ−d1−r1−r2−λ=0 is the characteristic polynomial. The eigen values are λ1=r,λ2=−t1−d2,λ3=−d3,λ4=−t2−d4,λ5=−d1−r1−r2. Thus, the trivial equilibrium point is locally asymptotically unstable due to the parameter r is nonnegative.Theorem 10 (AEP).
Axial equilibrium pointEAk,0,0,0,0 exists and always locally asymptotically stable in model (8) if and only if model parameters satisfy the conditions: βk−t1−d2<0,qp1k/s+k−d3<0,andqp3k/s+k−t2−d4<0. Otherwise, EA is locally asymptotically unstable.Proof.
Consider the variation matrix (40) at JEA(43)VEA=−r−r−βk−p1ks+k−p3ks+kr10βk−t1−d200000qp1ks+k−d3d3r2000qp3ks+k−t2−d400t10t2−d1−r1−r2.
Then, it is possible to find eigenvalue of this variation matrix asdetVEA−λI5=0,
(44)−r−λ−r−βk−p1ks+k−p3ks+kr10βk−t1−d2−λ00000qp1ks+k−d3−λd3r2000qp3ks+k−t2−d4−λ00t10t2−d1−r1−r2−λ=0.
Use first column to find the determinant,(45)⟹−r−λβk−t1−d2−λ0000qp1ks+k−d3−λd3r200qp3ks+k−t2−d4−λ0t10t2−d1−r1−r2−λ=0.
Use first row to compute the remaining determinant,(46)⟹−r−λβk−t1−d2−λqp1ks+k−d3−λd3r20qp3ks+k−t2−d4−λ00t2−d1−r1−r2−λ=0.
Use first column to get the characteristic polynomial,(47)⟹−r−λβk−t1−d2−λqp1ks+k−d3−λqp3ks+k−t2−d4−λ0t2−d1−r1−r2−λ=0,⟹−r−λβk−t1−d2−λqp1ks+k−d3−λqp3ks+k−t2−d4−λ−d1−r1−r2−λ=0,is characteristic polynomial. Then, eigenvalues are
(48)λ1=−r,λ2=−d1−r1−r2,λ3=βk−t1−d2,λ4=qp1ks+k−d3,λ5=qp3ks+k−t2−d4.
The axial equilibrium pointEA is locally asymptotically stable, if
(49)βk−t1−d2<0,qp1ks+k−d3<0,qp3ks+k−t2−d4<0.
Otherwise,EA is locally asymptotically unstable.Theorem 11 (DFEP).
The disease-free equilibrium pointE¯X¯,0,Y¯,0,0=d3s/qp1−d3,0,rsqkqp1−kd3−d3s/qp1−d32k,0,0 exists and is always locally asymptotically stable if and only if the model parameter satisfy conditions (i) βX¯−t1−d2−sp2Y¯/s2≤0, (ii) qp3X¯/s+X−αY¯−t2−d4≤0, (iii) r−2rX¯/k−sp1Y¯/s+X¯2+qp1X¯/s+X¯−d3>0, and (iv) r−2rX¯/k−sp1Y¯/s+X¯2qp1X¯/s+X¯−d3+sqp12X¯2/s+X¯3>0.Proof.
Consider the variation matrix (40) at disease-free equilibrium point E¯X¯,0,Y,¯0,0 such that
(50)VE¯=r−2rX¯k−sp1Y¯s+X¯2−rX¯k−βX¯−p1X¯s+X¯−p3X¯s+X¯r10βX¯−t1−d2−sp2Y¯s2000sqp1X¯s+X¯2sqp2Y¯s+X¯2qp1X¯s+X¯−d3−αY¯−d3r2000qp3X¯s+X¯−αY¯−t2−d400t10t2−d1−r1−r2.
Then, it is possible to find the determinant of the variation matrix asdetVE¯−λI5(51)r−2rX¯k−sp1Y¯s+X¯2−λ−rX¯k−βX¯−p1X¯s+X¯−p3X¯s+X¯r10βX¯−t1−d2−sp2Y¯s2−λ000sqp1X¯s+X¯2sqp2Y¯s+X¯2qp1X¯s+X¯−d3−λ−αY¯−d3r2000qp3X¯s+X¯−αY¯−t2−d4−λ00t10t2−d1−r1−r2−λ.
To find eigenvalues, compute the determinant using second row(52)βX¯−t1−d2−sp2Y¯s2−λ∗r−2rX¯k−sp1Y¯s+X¯2−λ−p1X¯s+X¯−p3X¯s+X¯r1sqp1X¯s+X¯2qp1X¯s+X¯−d3−λ−αY−d3r200qp3X¯s+X¯−αY−t2−d4−λ0000−d1−r1−r2−λ.
Now again use third row to find determinant:(53)βX¯−t1−d2−sp2Y¯s2−λqp3X¯s+X¯−αY¯−t2−d4−λ∗r−2rX¯k−sp1Y¯s+X¯2−λ−p1X¯s+X¯r1sqp1X¯s+X¯2qp1X¯s+X¯−d3−λr200−d1−r1−r2−λ.
Then, use third row to find characteristic polynomial(54)βX¯−t1−d2−sp2Y¯s2−λqp3X¯s+X¯−αY¯−t2−d4−λ−d1−r1−r2−λ∗r−2rX¯k−sp1Y¯s+X¯2−λ−p1X¯s+Xsqp1X¯s+X¯2qp1X¯s+X¯−d3−λβX¯−t1−d2−sp2Y¯s2−λqp3X¯s+X−αY¯−t2−d4−λ−d1−r1−r2−λ∗r−2rX¯k−sp1Y¯s+X¯2−λqp1X¯s+X¯−d3−λ+sqp12X¯2s+X¯3=0.
Eigenvalues are(55)λ1=βX¯−t1−d2−sp2Y¯s2,λ2=qp3X¯s+X¯−αY¯−t2−d4,λ3=−d1−r1−r2,
and the remaining eigenvalues can be obtained from the roots of quadratic equation:(56)r−2rX¯k−sp1Y¯s+X¯2⏟a−λqp1X¯s+X¯−d3⏟b−λ+sqp12X¯2s+X¯3⏟c=0.
It is known that a quadratic equationa−λb−λ+c=0 is locally asymptotically stable if and only if a+b>0andab+c>0, using such Routh-Hurwitz criterion, the disease-free equilibrium point E¯X¯,0,Y,¯0,0 is locally asymptotically stable if
(57)r−2rX¯k−sp1Y¯s+X¯2+qp1X¯s+X¯−d3>0,r−2rX¯k−sp1Y¯s+X¯2qp1X¯s+X¯−d3+sqp12X¯2s+X¯3>0,λ1=βX¯−t1−d2−sp2Y¯s2≤0,λ2=qp3X¯s+X−αY¯−t2−d4≤0,where X¯=d3s/qp1−d3 and Y¯=rsqkqp1−kd3−d3s/qp1−d32k. Otherwise, the disease-free equilibrium point is asymptotically unstable.Now let us see again the global stability analysis of model (8) around the endemic equilibrium point or positive equilibrium point E∗X∗,W∗,Y∗,Z∗,H∗ which shows coexistence. For that, let us state the following theorem and prove by taking appropriate Lyapunov function L.Theorem 12 (global stability).
Endemic equilibrium pointE∗X∗,W∗,Y∗,Z∗,H∗ exists and globally asymptotically stable.Proof.
Define appropriate Lyapunov functionLX,W,Y,Z,H=1/2X−X∗2+α1/2W−W∗2+α2/2Y−Y∗2+α3/2Z−Z∗2+α4/2H−H∗2, where α1,α2,α3,α4>0 are chosen properly such that dL/dt=0∀X∗,W∗,Y∗,Z∗,H∗∈ℝ+5 and dL/dt≤0∀X,W,Y,Z,H∈ℝ+5.This implies E∗ of the system is Lyapunov stable and dL/dt<0, ∀X,W,Y,Z,H∈ℝ+5 near E∗. This implies E∗ is globally asymptotically stable point. Now differentiate the Lyapunov function L with respect to t(58)dLdt=X−X∗dXdt+α1W−W∗dIdt+α2Y−Y∗dYdt+α3Z−Z∗dZdt+α4H−H∗dHdt.
Now substitute the model equation (8) into (58), we have the following equation:
(59)dLdt=X−X∗rX1−X+Wk+r1H−βXW−p1XYs+X−p3XZs+X+α1W−W∗βXW−t1W−d2W−p2WYs+W−p4WZs+W+α2Y−Y∗qp1XYs+X+qp2WYs+W+r2H−αYZ−d3Y+α3Z−Z∗qp3XZs+X+qp4WZs+W+αYZ−t2Z−d4Z+α4H−H∗t1W+t2Z−d1H−r1H−r2H.
Take outX,W,Y,Z,H from each bracket and write a change as follows
(60)dLdt=X−X∗X−X∗r1−X+Wk+r1HX−βW−p1Ys+X−p3Zs+X+α1W−W∗W−W∗βX−t1−d2−p2Ys+W−p4Zs+W+α2Y−Y∗Y−Y∗qp1Xs+X+qp2Ws+W+r2HY−αZ−d3+α3Z−Z∗Z−Z∗qp3Xs+X+qp4Ws+W+αY−t2−d4+α4H−H∗H−H∗t1WH+t2ZH−d1−r1−r2.
By rearranging, it is obtained that(61)dLdt=−X−X∗2−r1−X+Wk−r1HX+βW+p1Ys+X+p3Zs+X−α1I−I∗2−βX+t1+d2+p2Ys+W+p4Zs+W−α2Y−Y∗2−qp1Ss+X−qp2Ws+W−r2HY+αZ+d3−α3Z−Z∗2−qp3Xs+X−qp4Ws+W−αY+t2+d4−α4H−H∗2−t1WH−t2ZH+d1+r1+r2.
Thus, it is possible to setα1,α2,α3,α4 such that dL/dt≤0 and endemic equilibrium point E∗ is globally stable point.


## 5. Basic Reproduction Number

The basic reproduction number denoted byR0 and defined as the expected number of people getting secondary infection among the whole susceptible population. This number shows a potential for spread of disease within a given population. When R0<1, each infected individual produces on average less than one new infected individual so that the disease is expected to die out. On the other hand, if R0>1, then each individual produces more than one new infected individual so that the disease is expected to continue spreading in the population.Theorem 13.
The basic reproduction number for infected prey at disease-free equilibrium point (DFEP)E¯X¯,0,Y¯,0,0=d3s/qp1−d3,0,rsqkqp1−kd3−d3s/qp1−d32k,0,0 is given by
(62)R01=qp1−d32kβd3s2qp1−d3qp1−d32kst1+d2+rsqp2kqp1−kd3−d3s.Proof.
Consider infected prey equation in (8), we have the following model equation
(63)dWdt=βXW−t1W−d2W−p2WYs+W−p2WZs+W=βX−t1+d2+p2Ys+W+p2Zs+WW.
Now let us define functionsF and V, F=βX,V=t1+d2+p2Y/s+W+p2Z/s+W.
EvaluateF and V at DFEP E¯X¯,0,Y¯,0,0=d3s/qp1−d3,0,rsqkqp1−kd3−d3s/qp1−d32k,0,0,
(64)FE¯=βX¯=βd3sqp1−d3,VE¯=t1+d2+p2Y¯s=t1+d2+rsqp2kqp1−kd3−d3sqp1−d32ks=qp1−d32t1+d2ks+rsqp2kqp1−kd3−d3sqp1−d32ks.
Then, the basic reproduction number of infected prey is(65)R01=FV−1=βd3sqp1−d3∗qp1−d32ksqp1−d32t1+d2ks+rsqp2kqp1−kd3−d3s,R01=qp1−d32kβd3s2qp1−d3qp1−d32t1+d2ks+rsqp2kqp1−kd3−d3s.Theorem 14.
The basic reproduction number for infected predators at disease-free equilibrium point (DFEP)E¯X¯,0,Y¯,0,0=d3s/qp1−d3,0,rsqkqp1−kd3−d3s/qp1−d32k,0,0 is given by
(66)R02=qp1−d3qp3d3k+αrsqkqp1−kd3−d3sqp1−d32t2+d4k.Proof.
Consider the infected predator model equation in (8), we have the following equation
(67)dZdt=qp2WZs+W+qp3XZs+W+αYZ−t2Z−d4Z=qp2Ws+W+qp3Xs+W+αY−t2+d4Z.
Now let us define functionsF and V as F=qp2W/s+W+qp3X/s+W+αY,V=t2+d4. Then, evaluate F and V at (DFEP) E¯X¯,0,Y¯,0,0=d3s/qp1−d3,0,rsqkqp1−kd3−d3s/qp1−d32k,0,0,
(68)FE¯=qp3d3ssqp1−d3+αrsqkqp1−kd3−d3sqp1−d32k=qp1−d3qp3d3k+αrsqkqp1−kd3−d3sqp1−d32k,VE¯=t2+d4.
Therefore, the basic reproduction number of infected predator isR02=FV−1, and hence,
(69)R02=qp1−d3qp3d3k+αrsqkqp1−kd3−d3sqp1−d32t2+d4k.


## 6. Simulation

In this section, simulation of model (8) is carried out using DEDiscover version: 2.6.4. software. For simulation, a set of meaningful values are assigned to model parameters in Table 4, and initial values for model variables are given in Table 5. The model is arranged in such a way for simulation purposes.
(70)dXdt=r∗X∗1−X+Wk+r¯1∗H−beta∗X∗W−P¯1∗X∗Ys+X−P¯3∗X∗Zs+X//susceptibleprey,dWdt=beta∗X∗W−t¯1∗W−d¯2∗W−P¯2∗W∗Ys+W−P¯4∗W∗Zs+W//infectedprey,dYdt=q∗P¯1∗X∗Ys+X+q∗P¯2∗W∗Ys+W+r¯2∗H−alpha∗Y∗Z−d¯3∗Y//susceptiblepredator,dZdt=q∗P¯3∗X∗Zs+X+q∗P¯4∗W∗Zs+W−alpha∗Y∗Z−t¯2∗Z−d¯4∗Z//infectedpredator,dHdt=t¯1∗W+t¯2∗Z−d¯1∗H−r¯1∗H−r¯2∗H//bothinfectedpopulationsundertreatment.Table 4
Parameter value used in simulation.
NameValueDescriptionr22.4000Growth rate of susceptible preyk1.0000E03Carrying capacity of susceptible preyr¯11.0000Recovery rate ofBeta2.4000Disease transmission rate in preyP¯11.0000Predation coefficient of susceptible prey due to susceptible predators1.0000Half-saturated rateP¯31.0000Predation coefficient of susceptible prey due to infected predatort¯11.0000Treatment rate of infected preyd¯21.0000Death rate of infected preyP¯21.0000Predation coefficient of infected prey due to predatorsP¯41.0000Predation rate of infected prey due to infected predatorq1.0000Efficiency of predationr¯21.0000Recovery rate of susceptible predatorAlpha2.6000Disease transmission rate in predatord¯31.0000Death rate of susceptible predatort¯21.0000Treatment rate of infected predatord¯41.0000Death rate of infected predatord¯11.0000Death rate of both infected and infected predator under treatmentTable 5
Initial conditions used for model variables.
NameValueDescriptionXt01.2000E04Initial # susceptible preyWt0200.0000Initial # infected preyHt01.0000Initial # under treated prey predatorYt0160.0000Initial # susceptible predatorZt0180.0000Initial # of infected predatorFrom simulation in Figures2 and 3, it can be concluded that treatment is a helpful tool to minimize or eradicate infection. It is shown that as the treatment rate increases on infected prey-predator, then the infected prey-predator population decreases rapidly. This shows due to the fact that the infected prey-predator population is recovering and moves to susceptible classes and that contributes to the susceptible prey-predator population to rise in number.Figure 2
Infected prey-predator with treatment.Figure 3
Plots forW, Z, and H.In the sample simulation in Figures4 and 5, it is shown that high infection and predation results in the whole prey-predator population declining to a certain level. Therefore, it is better to implement treatment mechanisms to sustain stability of the prey-predator system.Figure 4
High infection and predation.Figure 5
Individual plot forW, Y, and Z.


## 7. Conclusions and Recommendation

In this paper, It can be concluded that the formulated model is mathematically meaningful, valid, and biologically well posed by proving the boundedness, positivity, and existence of the solutions of the model. Trivial, axial, disease-free, and endemic equilibrium points are investigated. Moreover, It is observed that in our model, trivial equilibrium point is always locally asymptotically unstable. Axial equilibrium point is locally asymptotically stable if and only if the variables satisfy the following three conditions: (i)βk−t1+d2<0, (ii) qp1k−d3s+k<0, and (iii) qp3k−t2+d4s+k<0.Treatment is a helpful tool to minimize or eradicate infection in the prey-predator system. Therefore, providing treatment in an infected prey-predator system creates opportunity to recover from illness and the prey-predator population can be saved and exists in stable situation. Thus, it is recommended to apply treatment on infected prey-predator to make the whole prey-predator population safe and abundant in nature. One can extend this paper by assuming the predator grows logistically or by adding parameter like death rate on the prey or by including other variables like vaccination, immigration, migration on prey-predator system, and these things can be considered as limitation of this paper.



---

*Source: 6679686-2021-02-22.xml*