Study of LG-Holling type III predator-prey model with disease in predator

TL;DR

This paper develops and analyzes a Leslie-Gower Holling type III predator-prey model incorporating disease in predators, examining stability, persistence, and boundedness through mathematical analysis and numerical verification.

Contribution

It introduces a novel predator-prey model with disease dynamics and provides comprehensive stability and persistence analysis with numerical validation.

Findings

Local and global stability conditions established

Boundedness and existence of solutions proven

Numerical simulations confirm analytical results

Abstract

In this article, a Leslie-Gower Holling type III predator-prey model with disease in predator has been developed from both biological and mathematical point of view. The total population is divided into three classes, namely, prey, susceptible predator and infected predator. The local stability, global stability together with sufficient conditions for persistence of the ecosystem near biologically feasible equilibria is thoroughly investigated. Boundedness and existence of the system are established. All the important analytical findings are numerically verified using program software MATLAB and Maple.

Tables

Table 1. Table 1: The set of model parameters and variables, dimension and their biological description.

Variable or Parameter	Unit or Dimension	Description
$x$	$V$	Prey density
$y$	$V$	Density of susceptible predator
$z$	$V$	Density of infected predator
$a_1$	$T^{-1}$	Intrinsic growth rate of prey
$a_2$	$T^{-1}$	Intrinsic growth rate of susceptible predator
$a_3$	$T^{-1}$	Intrinsic growth rate of infected predator
$b_1$	$V^{-1}{T}^{-1}$	Intra-specific competition rate of prey
$c_1$	$T^{-1}$	Predation rate of susceptible predator
$c_2$	$T^{-1}$	Death rate due to intra-specific competition of susceptible predator
$c_3$	$T^{-1}$	Death rate due to intra-specific competition of infected predator
$\theta$	$V^{-1}{T}^{-1}$	Disease incidence rate
$k_1$	$V^2$	Half saturation constant of the prey
$k_2$	$V$	Measure of alternative food
$p$	$Dimensionless$	Constant lies between 0 to 1

Table 2. Table 2: Set of parameter values for numerical simulations; S ≡ 𝑆 absent S\equiv Parameter sets

$S$	$a_1$	$a_2$	$a_3$	$b_1$	$k_1$	$k_2$	$c_1$	$c_2$	$c_3$	$\theta$	$p$
$S_1$	$4.5$	$3.8$	$0.005$	$0.075$	$100$	$160$	$2.8$	$1.97$	$1.95$	$0.0937$	$0.047$
$S_2$	$4.5$	$3.8$	$0.005$	$0.075$	$100$	$20$	$2.8$	$1.97$	$0.005$	$0.0937$	$0.047$
$S_3$	$5.0$	$7.8$	$1.5$	$0.0005$	$50$	$55$	$1.7$	$1.05$	$1.0$	$0.0217$	$0.73$
$S_4$	$4.0$	$6.0$	$0.05$	$0.005$	$100$	$200$	$0.08$	$0.7$	$0.50$	$0.002537$	$0.93$

Table 3. Table 3: Schematic representation of our analytical results for set of parameter values S 1 subscript 𝑆 1 S_{1}

Eq. pts.	Existence	Feasibility	Stability	Figure
$E_0$	$(0,0,0)$	Feasible	Unstable	———
$E_1$	$(0,0,0.4102564103)$	Feasible	Unstable	——–
$E_2$	$(0,308.6294416,0)$	Feasible	Unstable	———-
$E_3$	$(0,5.207637519,35.24000437)$	Feasible	Unstable	——-
$E_4$	$(60,0,0)$	Feasible	Unstable	———
$E_5^I$	$(0.5159678886,309.6247096,0)$	Feasible	———	———
$E_5^{II}$	Does not exist	———–	———	———
$E_5^{III}$	Does not exist	———-	———-	———
$E_6^I$	$(59.98394610,0,0.5640614003)$	Feasible	Unstable	–
$E_6^{II}$	Does not exist	——	——	—–
$E_6^{III}$	Does not exist	——–	——-	—–
$E_{\ast }^I$	$(56.43479200,3.837197569,36.62450011)$	Feasible	Stability conditions are mentioned in section 5	Figure 1
$E_{\ast }^{II}$	Does not exist	——–	——–	——-
$E_{\ast }^{III}$	Does not exist	——–	——-	——–
$E_{\ast }^{IV}$	Does not exist	——–	——	——–

Table 4. Table 4: Schematic representation of our analytical results for set of parameter values S 2 subscript 𝑆 2 S_{2}

Eq. pts.	Existence	Feasibility	Stability	Figure
$E_0$	$(0,0,0)$	Feasible	Unstable	—
$E_1$	$(0,0,20)$	Feasible	Unstable	—
$E_2$	$(0,38.57868020,0)$	Feasible	Unstable	-
$E_3$	Does not exist	——–	——-	——–
$E_4$	$(60,0,0)$	Feasible	Unstable	——-
$E_5^I$	$(3.751381115,45.81484682,0)$	Feasible	Stability conditions are mentioned in section 4.6	Figure 8
$E_5^{II}$	Does not exist	—-	—-	—–
$E_5^{III}$	Does not exist	—–	——–	—-
$E_6^I$	$(57.70608071,0,77.70608071)$	Feasible	Stability conditions are mentioned in section 4.7	Figure 7
$E_6^{II}$	Does not exist	——	——	——
$E_6^{III}$	Does not exist	——–	——	—–
$E_{\ast }^I$	Does not exist	——–	——–	———
$E_{\ast }^{II}$	Does not exist	——–	——–	———
$E_{\ast }^{III}$	Does not exist	——	—–	——
$E_{\ast }^{IV}$	Does not exist	——	—-	——

Table 5. Table 5: Schematic representation of our analytical results for set of parameter values S 3 subscript 𝑆 3 S_{3}

Eq. pts.	Existence	Feasibility	Stability	Figure
$E_0$	$(0,0,0)$	Feasible	Unstable	—–
$E_1$	$(0,0,82.50000000)$	Feasible	Unstable	—–
$E_2$	$(0,408.5714286,0)$	Feasible	Unstable	—-
$E_3$	$(0,164.3476956,114.3012791)$	Feasible	Unstable	—–
$E_4$	$(10000,0,0)$	Feasible	Unstable	—-
$E_5^I$	$(0.3585095936,411.2346427,0)$	Feasible	Stability conditions are mentioned in section 4.6	Figure 3
$E_5^{II}$	Does not exist	——	——	—–
$E_5^{III}$	Does not exist	—–	—–	—-
$E_6^I$	$(6244.212048,0,9448.818072)$	Feasible	Stability conditions are mentioned in section 4.7	Figure 4
$E_6^{II}$	$(2.657590193,0,86.48638529)$	Feasible	Unstable	—-
$E_6^{III}$	$(30.13036196,0,127.6955429)$	Feasible	Unstable	——
$E_{\ast }^I$	$(0.5990751338,161.9321592,116.8375923)$	Feasible	Stability conditions are mentioned in section 5	Figure 2
$E_{\ast }^{II}$	$(73.69033011,33.00904773,252.2068593)$	Feasible	Unstable	—–
$E_{\ast }^{III}$	Does not exist	——-	——	—-
$E_{\ast }^{IV}$	Does not exist	—–	—–	—-

Table 6. Table 6: Schematic representation of our analytical results for set of parameter values S 4 subscript 𝑆 4 S_{4}

Eq. pts.	Existence	Feasibility	Stability	Figure
$E_0$	$(0,0,0)$	Feasible	Unstable	—–
$E_1$	$(0,0,20)$	Feasible	Unstable	——–
$E_2$	$(0,1714.285714,0)$	Feasible	Unstable	——-
$E_3$	$(0,1637.974544,44.24202326)$	Feasible	Unstable	—–
$E_4$	$(800,0,0)$	Feasible	Unstable	—–
$E_5^I$	$(3.142771393,1741.223755,0)$	Feasible	Unstable	—–
$E_5^{II}$	$(41.15227929,2067.019537,0)$	Feasible	Unstable	—-
$E_5^{III}$	$(618.5620922,7016.246504,0)$	Feasible	Unstable	——
$E_6^I$	$(798.1394249,0,99.81394249)$	Feasible	Unstable	——-
$E_6^{II}$	Does not exist	—-	——	——-
$E_6^{III}$	Does not exist	——-	—–	———
$E_{\ast }^I$	$(3.271695068,1618.749669,71.15684924)$	Feasible	Stability conditions are mentioned in section 5	Figure 5
$E_{\ast }^{II}$	$(756.7723630,426.6171136,1740.142426)$	Feasible	Stability conditions are mentioned in section 5	Figure 6
$E_{\ast }^{III}$	$(33.12282349,1461.965361,290.6548798)$	Feasible	Unstable	——-
$E_{\ast }^{IV}$	Does not exist	—-	——	——-

Equations

$$\frac{dx}{dt}=a_{1}x-b_1{x}^2-\frac{c_1{x}^{2}y}{k_1+x^2}-\frac{p{c}_1{x}^{2}z}{k_1+x^2}=f_1(x,y,z),$$

$$\frac{dy}{dt}=a_{2}y-\frac{c_{2}y(y+z)}{k_2+x}-\theta yz=f_2(x,y,z),$$

$$\frac{dz}{dt}=\theta yz+a_{3}z-\frac{c_{3}z(y+z)}{k_2+x}=f_3(x,y,z),$$

$$x(0)\ge 0,y(0)\ge 0,z(0)\ge 0,$$

$$\frac{d\omega }{dt}+\mu \omega$$

$$J_{5}=\left(\begin{array}[]{ccc}m_{11}&m_{12}&m_{13}\\ m_{21}&m_{22}&m_{23}\\ m_{31}&m_{32}&m_{33}\\ \end{array}\right),$$

$$J_{6}=\left(\begin{array}[]{ccc}m_{11}&m_{12}&m_{13}\\ m_{21}&m_{22}&m_{23}\\ m_{31}&m_{32}&m_{33}\\ \end{array}\right),$$

$$\frac{d{L}_1}{dt}$$

$$\frac{d{L}_2}{dt}$$

$$\frac{d{L}_3}{dt}$$

$$\frac{dL}{dt}$$

$$\Pi (x,y,z)$$

$$\displaystyle=\gamma_{1}\Big{(}a_{1}-b_{1}x-\frac{c_{1}xy}{x^{2}+k_{1}}-\frac{pc_{1}xz}{x^{2}+k_{1}}\Big{)}$$

$$\displaystyle+\gamma_{3}\Big{(}a_{3}-\frac{c_{3}(y+z)}{x+k_{2}}+\theta y\Big{)}.$$

$$\displaystyle\Pi\Big{(}0,\frac{a_{2}c_{3}-a_{3}c_{2}-a_{3}k_{2}\theta}{\theta(k_{2}\theta+c_{2}-c_{3})},-\frac{a_{2}c_{3}-a_{3}c_{2}-a_{2}k_{2}\theta}{\theta(k_{2}\theta+c_{2}-c_{3})}\Big{)}$$

$$=a_1{\gamma }_1>0.$$

$$\displaystyle\Pi\Big{(}x_{5},\frac{a_{2}(x_{5}+k_{2})}{c_{2}},0\Big{)}$$

$$\displaystyle+\gamma_{3}\Big{(}\frac{a_{3}c_{2}-a_{2}c_{3}+a_{2}k_{2}\theta+a_{2}x_{5}\theta}{c_{2}}\Big{)}>0\ \text{by the condition }3,$$

$$\displaystyle\Pi\Big{(}x_{6},0,\frac{a_{3}(x_{6}+k_{2})}{c_{3}}\Big{)}$$

$$\displaystyle+\gamma_{2}\Big{(}\frac{a_{2}c_{3}-a_{3}(c_{2}+(k_{2}+x_{6})\theta)}{c_{3}}\Big{)}>0\ \text{ by the condition }4.$$

Full text

Keywords: Eco-epidemic model, Intra-specific competition, Local and global stability, Lyapunov function, Persistence.

1 Introduction

The predators and the preys carry a dynamic relationship among themselves. And for its universal existence and importance, this relationship is one of the dominant themes in theoretical ecology. Mathematical modelling is considered to be very useful tool to understand and analyze the dynamic behavior of predator-prey systems. Predator functional response on prey population is the major element in predator-prey interaction. It describes the number of prey consumed per predator per unit time for given quantities of prey and predator. The most important and useful functional responses are Lotka-Volterra functional responses such as Holling type I functional response, Holling type II functional response and three species population models with such functional responses are widely researched in ecological literature [24], [23], [21], [19]. There are also many research works on three species systems like two preys one predator [1], [20], [12], [10], tritrophic food chain [5], [16], [2] etc.

The Mathematical modelling of epidemics has become a very important subject of research after the seminal model of Kermack-MacKendrick (1927) on SIRS (susceptible-infected-removed-susceptible) systems. It describes the evolution of a disease which gets transmitted upon contact. Important studies have been carried out with the aim of controlling the effects of diseases and of developing suitable vaccination strategies [4], [22], [3]. Eco-epidemic research describes disease that spread among interacting populations, where the epidemic and demographic aspects are merged within one model. During the last decade, this branch of science is developing and studied by the authors in [4], [15], [27]. In the natural world, species do not exist alone. It is of more biological significance to study persistence-extinction threshold of each population in systems of two or more interacting species subjected to parasitism. In mathematical biology the predator prey systems and models for transmissible disease are major field of study in their own right. In the growing ecoepidemic literature and from early papers [11], disease mainly spreading in the prey are examined in [29], [8], [9], but in [28], [17], [18], the epidemics are assumed to affect the predators. The predator-prey model with modified Lesli-Gower Holling type II Scheme was introduced in [6], [14], [26]. The LG model with Holling type II response function with disease in predator is discussed in [25]. But no one pay the attention for the modified LG model with Holling type III response function for predation with disease in predator.

Here we make an attempt to study the above said model with Holling type III response for predation and intra-specific competition among predators. The rest of the article is as follows. In Section $2$, we explain the formulation of the model under consideration and its assumptions. Section $3$ contains some preliminary results. In Section $4$, we analyze the system behavior of the trivial equilibria. Also the model with intra specific competition is analyzed for the system behavior around axial and boundary equilibria in Section $4$ . In Section $5$, local and global stability of the interior equilibria is analyzed. Section $6$ contains persistence of the system. Numerical simulation has been carried out in Section $7$ to support our analytical findings. The article comes to an end with a discussion of the results obtained in Section $8$.

2 Mathematical model formulation

We make the following assumptions:

• •

The disease spreads only among the predators. Let $y$ denotes the susceptible predators and $z$ the infected ones. The total predator population is $n(t)=y(t)+z(t)$.

• •

The disease spreads with a simple mass action law (with the disease incidence $\theta>0$). The prey population $x$ grows logistically with intrinsic growth rate $a_{1}>0$ and carrying capacity $a_{1}/b_{1}$ in the absence of predator population.

• •

We introduce intra-specific competition among the predator’s sound and infected sub-populations.

• •

Holling type-III response mechanism is considered for predation.

According to the above assumptions, we get the following model with non negative parameters

[TABLE]

where $a_{2},a_{3}(a_{2}\geq a_{3})$ are the per capita growth rates of each predator sub population. Thus from sick parents, the disease can be transmitted to their offspring. The parameter $k_{1}$ represents the half saturation constant of the prey and $k_{2}$ is the measure of alternative food. Hence the Jacobian matrix of the system (2.1) is $J=(m_{ij})\in\mathbb{R}^{3\times 3}$ with entries

$m_{11}=a_{1}-2b_{1}x-\frac{2c_{1}xy}{x^{2}+k1}+\frac{2c_{1}x^{3}y}{(x^{2}+k_{1})^{2}}-\frac{2pc_{1}xz}{x^{2}+k_{1}}+\frac{2pc_{1}x^{3}z}{(x^{2}+k_{1})^{2}},\qquad m_{12}=-\frac{c_{1}x^{2}}{x^{2}+k_{1}},\qquad m_{13}\\ =-\frac{pc_{1}x^{2}}{x^{2}+k_{1}},\qquad m_{21}=\frac{c_{2}y(y+z)}{(x+k_{2})^{2}},\qquad m_{22}=a_{2}-\frac{c_{2}(2y+z)}{x+k_{2}}-\theta z,\qquad m_{23}=-\frac{c_{2}y}{x+k_{2}}-\theta y,\qquad m_{31}=\frac{c_{3}z(y+z)}{(x+k_{2})^{2}},\qquad m_{32}=\theta z-\frac{c_{3}z}{x+k_{2}},\qquad m_{33}=\theta y+a_{3}-\frac{c_{3}(y+2z)}{x+k_{2}}$.

3 Preliminary results

3.1 Existence

Theorem 1**.**

Every solution of the system (2.1) with initial conditions exists in the interval $(0,+\infty)$ and $x(t)\geq 0$, $y(t)\geq 0$, $z(t)\geq 0$ for all $t\geq 0$.

Proof.

We have $\frac{dx}{dt}=f_{1}(x,y,z)$, $\frac{dy}{dt}=f_{2}(x,y,z)$, $\frac{dz}{dt}=f_{3}(x,y,z)$. Integrating we get $x(t)=x(0)e^{\int_{0}^{t}f_{1}(x,y,z)ds}$, $y(t)=y(0)e^{\int_{0}^{t}f_{2}(x,y,z)ds}$, $z(t)=z(0)e^{\int_{0}^{t}f_{3}(x,y,z)ds}$, where $x(0)=x_{0}>0$, $y(0)=y_{0}>0$, $z(0)=z_{0}>0$. Since $f_{1},\ f_{2},\ f_{3}$ are continuous function and hence locally Lipschitzian on $R^{3}_{+}$, the solution with positive initial condition exists and unique on $(0,\xi)$ where $0<\xi<\infty$. Hence the theorem. ∎

3.2 Boundedness

Theorem 2**.**

All the solutions of the system which initiate in $\mathbb{R}_{+}^{3}$ are uniformly bounded.

Proof.

Let us define a function $\omega=x+y+z$. Therefore, we have

[TABLE]

Hence we find $l>0$ such that $\frac{dw}{dt}+\mu\omega\leq l$, $\forall t\in(0,t_{b})$. Using the theory of differential inequality [7], we obtain 0<\omega(x,y,z)\leq\frac{l}{\mu}(1-e^{-ut})+\omega\Big{(}x(0),y(0),z(0)\Big{)}e^{-\mu t} and for $t_{b}\rightarrow\infty$ we have $0<\omega\leq\frac{l}{\mu}$.

Hence all the solutions of the system that initiate in $\mathbb{R}_{+}^{3}$ are confined in the region $\gamma=\{{(x,y,z)\in\mathbb{R}_{+}^{3}:\omega=\frac{l}{\mu}+\epsilon}\}$ for any $\epsilon>0$ and for $t$ large enough.

Hence the theorem. ∎

3.3 Equilibrium points

The system of equations (2.1) has the equilibrium points $E_{0}(0,0,0)$, $E_{1}(0,0,\frac{a_{3}k_{2}}{c_{3}})$, $E_{2}(0,\frac{a_{2}k_{2}}{c_{2}},0)$, $E_{3}(0,y_{3},z_{3})$, $E_{4}(\frac{a_{1}}{b_{1}},0,0)$, $E_{5}(x_{5},y_{5},0)$, $E_{6}(x_{6},0,z_{6})$ and $E_{*}(x_{*},y_{*},z_{*})$. The co-existence equilibrium is $E_{*}(x_{*},y_{*},z_{*})$ where $y_{*}=\frac{-a_{3}\theta x_{*}-a_{3}k_{2}\theta+a_{2}c_{3}-a_{3}c_{2}}{\theta(\theta x_{*}+k_{2}\theta+c_{2}-c_{3})}$ ,

$z_{*}=-\frac{-a_{2}\theta x_{*}-a_{2}k_{2}\theta+a_{2}c_{3}-a_{3}c_{2}}{\theta(\theta x_{*}+k_{2}\theta+c_{2}-c_{3})}$ and $x_{*}$ is root of the following equation

$P(x)=A_{4}x^{4}+A_{3}x^{3}+A_{2}x^{2}+A_{1}x+A_{0}=0$.

Here $A_{4}=b_{1}\theta^{2}$,

$A_{3}=b_{1}k_{2}\theta^{2}-a_{1}\theta^{2}+b_{1}c_{2}\theta-b_{1}c_{3}\theta$,

$A_{2}=-a_{1}k_{2}\theta^{2}+a_{2}c_{1}\theta p+b_{1}k_{1}\theta^{2}-a_{1}c_{2}\theta+a_{1}c_{3}\theta-a_{3}c_{1}\theta$,

$A_{1}=a_{2}c_{1}k_{2}\theta p+b_{1}k_{1}k_{2}\theta^{2}-a_{1}k_{1}\theta^{2}-a_{2}c_{1}c_{3}p+a_{3}c_{1}c_{2}p-a_{3}c_{1}k_{2}\theta+b_{1}c_{2}k_{1}\theta-b_{1}c_{3}k_{1}\theta+a_{2}c_{1}c_{3}-a_{3}c_{1}c_{2}$,

$A_{0}=-a_{1}k_{1}k_{2}\theta^{2}-a_{1}c_{2}k_{1}\theta+a_{1}c_{3}k_{1}\theta$.

We consider $P(x)=(\gamma_{1}x^{2}+\delta_{1}x+\phi_{1})(\gamma_{2}x^{2}+\delta_{2}x+\phi_{2})$.

Case-1: $P(x)$ has two real roots if either $\delta_{1}^{2}-4\gamma_{1}\phi_{1}>0$ or $\delta_{2}^{2}-4\gamma_{2}\phi_{2}>0$.

For the set of parameters $a_{1}=4.5,\ a_{2}=3.8,\ a_{3}=0.005,\ b_{1}=0.075,\ k_{1}=100,\ k_{2}=160,\ c_{1}=2.8,\ c_{2}=1.97,\ c_{3}=1.95,\ \theta=0.0937,\ p=0.047$, there are two real roots $[x=56.43479200,\ y=3.837197569,\ z=36.62450011]$ and $[x=-161.0505378,\ y=-1006.972848,\ z=1057.801828]$ in which first one is biologically feasible.

Case-2: $P(x)$ has four real roots if $\delta_{1}^{2}-4\gamma_{1}\phi_{1}>0$ and $\delta_{2}^{2}-4\gamma_{2}\phi_{2}>0$.

For the set of parameters $a_{1}=5.0,\ a_{2}=7.8,\ a_{3}=1.5,\ b_{1}=0.0005,\ k_{1}=50,\ k_{2}=55,\ c_{1}=1.7,\ c_{2}=1.05,\ c_{3}=1.0,\ \theta=0.0217,\ p=0.73$, there are four real roots $[x=0.5990751338,\ y=161.9321592,\ z=116.8375923]$, $[x=-65.33582457,\ y=-1734.893676,\ z=2108.504719]$, $[x=73.69033011,\ y=33.00904773,\ z=252.2068593]$, $[x=9933.742272,\ y=-67.78533292,\ z=358.0409590]$ in which first and third one are biologically feasible.

For the equilibrium point $E_{5}(x_{5},y_{5},0)$, $z_{5}=0$ gives $y_{5}=\frac{a_{2}(x_{5}+k_{2})}{c_{2}}$. Here $x_{5}$ is the root of

$Q(x)=A_{3}x^{3}+A_{2}x^{2}+A_{1}x+A_{0}=0$,

where $A_{3}=b_{1}c_{2}$, $A_{2}=-a_{1}c_{2}+a_{2}c_{1}$, $A_{1}=a_{2}c_{1}k_{2}+b_{1}c_{2}k_{1}$, $A_{0}=-a_{1}c_{2}k_{1}$.

We consider $Q(x)=(\mu_{1}x+\alpha_{1})(\nu_{1}x^{2}+\xi_{1}x+\eta_{1})$.

Case-1: $Q(x)$ has only one real root if $\xi_{1}^{2}-4\nu_{1}\eta_{1}<0$, which yields $x_{5}=-\frac{\alpha_{1}}{\mu_{1}}$.

For the set of parameters $a_{1}=4.5,\ a_{2}=3.8,\ a_{3}=0.005,\ b_{1}=0.075,\ k_{1}=100,\ k_{2}=160,\ c_{1}=2.8,\ c_{2}=1.97,\ c_{3}=1.95,\ \theta=0.0937,\ p=0.047$, there is only one real root $[x=0.5159678886,\ y=309.6247096,\ z=0]$ which is biologically feasible.

Case-2: $Q(x)$ has three real roots if $\xi_{1}^{2}-4\nu_{1}\eta_{1}>0$.

For the set of parameters $a_{1}=5.0,\ a_{2}=7.8,\ a_{3}=1.5,\ b_{1}=0.0005,\ k_{1}=50,\ k_{2}=55,\ c_{1}=1.7,\ c_{2}=1.05,\ c_{3}=1.0,\ \theta=0.0217,\ p=0.73$, there are three real roots $[x=0.3585095936,\ y=411.2346427,\ z=0]$, $[x=-91.96262935,\ y=-274.5795323,\ z=0]$, $[x=-15165.53874,\ y=-1.122497163*10^{5},\ z=0]$ in which first one is biologically feasible.

For the equilibrium point $E_{6}(x_{6},0,z_{6})$, $y_{6}=0$ gives $z_{6}=\frac{a_{3}(x_{6}+k_{2})}{c_{3}}$. Here $x_{6}$ is the root of

$R(x)=A_{3}x^{3}+A_{2}x^{2}+A_{1}x+A_{0}=0$,

where $A_{3}=b_{1}c_{3}$, $A_{2}=-a_{1}c_{3}+a_{3}c_{1}p$, $A_{1}=a_{3}c_{1}k_{2}p+b_{1}c_{3}k_{1}$, $A_{0}=-a_{1}c_{3}k_{1}$.

We consider $R(x)=(\mu_{2}x+\alpha_{2})(\nu_{2}x^{2}+\xi_{2}x+\eta_{2})$.

Case-1: $R(x)$ has only one real root if $\xi_{2}^{2}-4\nu_{2}\eta_{2}<0$, which yields $x_{6}=-\frac{\alpha_{2}}{\mu_{2}}$.

For the set of parameters $a_{1}=4.5,\ a_{2}=3.8,\ a_{3}=0.005,\ b_{1}=0.075,\ k_{1}=100,\ k_{2}=160,\ c_{1}=2.8,\ c_{2}=1.97,\ c_{3}=1.95,\ \theta=0.0937,\ p=0.047$, there is only one real root $[x=59.98394610,\ y=0,\ z=0.5640614003]$ which is biologically feasible.

Case-2: $R(x)$ has three real roots if $\xi_{2}^{2}-4\nu_{2}\eta_{2}>0$.

For the set of parameters $a_{1}=5.0,\ a_{2}=7.8,\ a_{3}=1.5,\ b_{1}=0.0005,\ k_{1}=50,\ k_{2}=55,\ c_{1}=1.7,\ c_{2}=1.05,\ c_{3}=1.0,\ \theta=0.0217,\ p=0.73$, there are three real roots $[x=2.657590193,y=0,z=86.48638529]$, $[x=30.13036196,\ y=0,\ z=127.6955429]$, $[x=6244.212048,\ y=0,\ z=9448.818072]$ which are biologically feasible.

4 System behaviour around boundary equilibria

4.1 Stability for $E_{0}$

The characteristic equation for $E_{0}$ is given by

\left|\begin{array}[]{ccc}a_{1}-\lambda&0&0\\ 0&a_{2}-\lambda&0\\ 0&0&a_{3}-\lambda\\ \end{array}\right|=0.

The equilibrium point $E_{0}$ has the eigenvalues $a_{1}$, $a_{2}$, $a_{3}$. All the eigenvalues are positive and it is unstable.

4.2 Stability for $E_{1}$

The characteristic equation for $E_{1}$ is given by

\left|\begin{array}[]{ccc}a_{1}-\lambda&0&0\\ 0&(\frac{a_{2}c_{3}-a_{3}(c_{2}+k_{2}\theta)}{c_{3}})-\lambda&0\\ \frac{a_{3}^{2}}{c_{3}}&a_{3}(-1+\frac{k_{2}\theta}{c_{3}})&-a_{3}-\lambda\\ \end{array}\right|=0.

Since one of the eigenvalues of $E_{1}$ is $a_{1}$, which is always positive and so, $E_{1}$ is unstable.

4.3 Stability for $E_{2}$

The characteristic equation for $E_{2}$ is given by

\left|\begin{array}[]{ccc}a_{1}-\lambda&0&0\\ \frac{a_{2}^{2}}{c_{2}}&-a_{2}-\lambda&-a_{2}-\frac{a_{2}k_{2}\theta}{c_{2}}\\ 0&0&(a_{3}-\frac{a_{2}c_{3}}{c_{2}}+\frac{a_{2}k_{2}\theta}{c_{2}})-\lambda\\ \end{array}\right|=0.

Since one of the eigenvalues of $E_{2}$ is $a_{1}$, which is always positive and therefore, $E_{2}$ is unstable.

4.4 Stability for $E_{3}$

The characteristic equation of the equilibrium $E_{3}(0,y_{3},z_{3})$ is given by

\left|\begin{array}[]{ccc}a_{1}-\lambda&0&0\\ \frac{c_{2}y_{3}(y_{3}+z_{3})}{k_{2}^{2}}&(a_{2}-\frac{c_{2}(2y_{3}+z_{3})}{k_{2}}-z_{3}\theta)-\lambda&-\frac{c_{2}y_{3}}{k_{2}}-y_{3}\theta\\ \frac{c_{3}z_{3}(y_{3}+z_{3})}{k_{2}^{2}}&-\frac{c_{3}z_{3}}{k_{2}}+z_{3}\theta&(a_{3}-\frac{c_{3}(y_{3}+2z_{3})}{k_{2}}+y_{3}\theta)-\lambda\\ \end{array}\right|=0,

where $y_{3}=\frac{-a_{3}k_{2}\theta+a_{2}c_{3}-a_{3}c_{2}}{\theta(k_{2}\theta+c_{2}-c_{3})}$ and $z_{3}=-\frac{-a_{2}k_{2}\theta+a_{2}c_{3}-a_{3}c_{2}}{\theta(k_{2}\theta+c_{2}-c_{3})}$. As eigenvalue $a_{1}$ for $E_{3}$ is always positive, the equilibrium point $E_{3}$ is unstable.

4.5 Stability for $E_{4}$

The characteristic equation for $E_{4}$ is given by

\left|\begin{array}[]{ccc}-a_{1}-\lambda&-\frac{a_{1}^{2}c_{1}}{b_{1}^{2}(\frac{a_{1}^{2}}{b_{1}^{2}}+k_{1})}&-\frac{a_{1}^{2}c_{1}p}{b_{1}^{2}(\frac{a_{1}^{2}}{b_{1}^{2}}+k_{1})}\\ 0&a_{2}-\lambda&0\\ 0&0&a_{3}-\lambda\\ \end{array}\right|=0.

Eigenvalues $a_{2}$ and $a_{3}$ are always positive of the equilibrium $E_{4}$. Hence $E_{4}$ is unstable.

4.6 Stability for $E_{5}(x_{5},y_{5},0)$

At $E_{5}(x_{5},y_{5},0)$, the Jacobian matrix for the system is given by

[TABLE]

where $m_{11}=-b_{1}x_{5}-\frac{c_{1}x_{5}y_{5}}{x_{5}^{2}+k_{1}}+\frac{2c_{1}x_{5}^{3}y_{5}}{(x_{5}^{2}+k_{1})^{2}}$,  $m_{12}=-\frac{c_{1}x_{5}^{2}}{x_{5}^{2}+k_{1}}$,  $m_{13}=-\frac{pc_{1}x_{5}^{2}}{x_{5}^{2}+k_{1}}$,  $m_{21}=\frac{c_{2}y_{5}^{2}}{(x_{5}+k_{2})^{2}}$,  $m_{22}=-\frac{c_{2}y_{5}}{x_{5}+k_{2}}$,  $m_{23}=-\frac{c_{2}y_{5}}{x_{5}+k_{2}}-\theta y_{5}$,  $m_{31}=0$,  $m_{32}=0$,  $m_{33}=\theta y_{5}+a_{3}-\frac{c_{3}y_{5}}{x_{5}+k_{2}}$.

The Characteristic equation for $J_{5}$ is given by

$(m_{33}-\lambda)\{(m_{11}-\lambda)(m_{22}-\lambda)-m_{21}m_{12}\}=0$

$\Rightarrow(m_{33}-\lambda)\{\lambda^{2}-(m_{11}+m_{22})\lambda+m_{11}m_{22}-m_{21}m_{12}\}=0$

$\Rightarrow\lambda_{1,2}=\frac{m_{11}+m_{22}\pm\sqrt{(m_{11}+m_{22})^{2}-4(m_{11}m_{22}-m_{21}m_{12})}}{2}$ and $\lambda_{3}=m_{33}$.

We choose $m_{11}<0$ and $m_{33}<0$. Then $E_{5}$ will be stable if

(i) $b_{1}x_{5}+\frac{c_{1}x_{5}y_{5}}{x_{5}^{2}+k_{1}}>\frac{2c_{1}x_{5}^{3}y_{5}}{(x_{5}^{2}+k_{1})^{2}}$,

(ii)$\frac{c_{3}y_{5}}{x_{5}+k_{2}}>\theta y_{5}+a_{3}$.

4.7 Stability for $E_{6}(x_{6},0,z_{6})$

At $E_{6}(x_{6},0,z_{6})$, the Jacobian matrix for the system is given by

[TABLE]

where $m_{11}=-b_{1}x_{6}-\frac{pc_{1}x_{6}z_{6}}{x_{6}^{2}+k_{1}}+\frac{2pc_{1}x_{6}^{3}z_{6}}{(x_{6}^{2}+k_{1})^{2}}$,  $m_{12}=-\frac{c_{1}x_{6}^{2}}{x_{6}^{2}+k_{1}}$,  $m_{13}=-\frac{pc_{1}x_{6}^{2}}{x_{6}^{2}+k_{1}}$,  $m_{21}=0$,  $m_{22}=a_{2}-\frac{c_{2}z_{6}}{x_{6}+k_{2}}-\theta z_{6}$,  $m_{23}=0$,  $m_{31}=\frac{c_{3}z_{6}^{2}}{(x_{6}+k_{2})^{2}}$,  $m_{32}=\theta z_{6}-\frac{c_{3}z_{6}}{x_{6}+k_{2}}$,  $m_{33}=-\frac{c_{3}z_{6}}{x_{6}+k_{2}}$.

The Characteristic equation for $J_{6}$ is given by

$(m_{22}-\lambda)\{(m_{11}-\lambda)(m_{33}-\lambda)-m_{31}m_{13}=0$

$\Rightarrow(m_{22}-\lambda)\{{\lambda^{2}-(m_{11}+m_{33})\lambda+m_{11}m_{33}-m_{31}m_{13}}\}=0$

$\Rightarrow\lambda_{1,2}=\frac{m_{11}+m_{33}\pm\sqrt{(m_{11}+m_{33})^{2}-4(m_{11}m_{33}-m_{31}m_{13})}}{2}$ and $\lambda_{3}=m_{22}$.

We choose $m_{22}<0$ and $m_{11}<0$. Then $E_{6}$ will be stable if

(i) $b_{1}x_{6}+\frac{pc_{1}x_{6}z_{6}}{x_{6}^{2}+k_{1}}>\frac{2pc_{1}x_{6}^{3}z_{6}}{(x_{6}^{2}+k_{1})^{2}}$,

(ii)$\frac{c_{2}z_{6}}{x_{6}+k_{2}}+\theta z_{6}>a_{2}$.

5 System behaviour near the coexistence equilibrium $E_{*}(x_{*},y_{*},z_{*})$

The entries for the Jacobian at $E_{*}(x_{*},y_{*},z_{*})$ are

$m_{11}=a_{1}-2b_{1}x_{*}+\frac{2c_{1}x_{*}^{3}y_{*}}{(x_{*}^{2}+k_{1})^{2}}+\frac{2pc_{1}x_{*}^{3}z_{*}}{(x_{*}^{2}+k_{1})^{2}}$, $m_{12}=-\frac{c_{1}x_{*}^{2}}{x_{*}^{2}+k_{1}}$, $m_{13}=-\frac{pc_{1}x_{*}^{2}}{x_{*}^{2}+k_{1}}$,

$m_{21}=\frac{c_{2}y_{*}(y_{*}+z_{*})}{(x_{*}+k_{2})^{2}}$, $m_{22}=-\frac{c_{2}y_{*}}{x_{*}+k_{2}}$, $m_{23}=-\frac{c_{2}y_{*}}{x_{*}+k_{2}}-\theta y_{*}$,

$m_{31}=\frac{c_{3}z_{*}(y_{*}+z_{*})}{(x_{*}+k_{2})^{2}}$, $m_{32}=\theta z_{*}-\frac{c_{3}z_{*}}{x_{*}+k_{2}}$, $m_{33}=-\frac{c_{3}z_{*}}{x_{*}+k_{2}}$.

5.1 Local Stability

The characteristic equation for $J_{*}$ is $\theta^{3}+A_{1}\theta^{2}+A_{2}\theta+A_{3}=0$, where

$A_{1}=-m_{11}-m_{22}-m_{33}$,

$A_{2}=m_{11}m_{22}+m_{11}m_{33}+m_{22}m_{33}-m_{13}m_{31}-m_{23}m_{32}-m_{21}m_{12}$,

$A_{3}=m_{11}m_{23}m_{32}+m_{12}m_{21}m_{33}+m_{13}m_{22}m_{31}-m_{11}m_{22}m_{33}-m_{12}m_{23}m_{31}-m_{13}m_{21}m_{32}$,

$A_{1}A_{2}-A_{3}=-m_{11}^{2}m_{22}-m_{11}^{2}m_{33}-m_{22}^{2}m_{33}-m_{11}m_{22}^{2}-m_{11}m_{33}^{2}-m_{22}m_{33}^{2}-2m_{11}m_{22}m_{33}+m_{11}m_{13}m_{31}+m_{11}m_{12}m_{21}+m_{22}m_{12}m_{21}+m_{22}m_{23}m_{32}\\ +m_{33}m_{23}m_{32}+m_{33}m_{13}m_{31}+m_{12}m_{23}m_{31}+m_{13}m_{21}m_{32}$.

Assume $(i)$ $m_{11}<0$, $(ii)$ $m_{32}>0$, $(iii)$ $m_{11}m_{31}+m_{33}m_{31}+m_{21}m_{32}<0$ or $m_{11}m_{12}+m_{22}m_{12}+m_{13}m_{32}>0$ or $m_{22}m_{23}+m_{33}m_{23}+m_{13}m_{21}>0$ and we have

$A_{1}>0$, $A_{3}>0$ and $A_{1}A_{2}-A_{3}>0$. By Routh-Hurwitz criterion, the interior or co-existence equilibrium $E_{*}(x_{*},y_{*},z_{*})$ is locally asymptotically stable.

5.2 Global Stability

Theorem 3**.**

The co-existence equilibrium point $E_{*}$ is globally asymptotically stable if $P_{1}>0$, $P_{2}>0$ and $P_{3}>0$ where $P_{1}$ , $P_{2}$, $P_{3}$ are defined latter.

Proof.

Let us define the function $L(x,y,z)=L_{1}(x,y,z)+L_{2}(x,y,z)+L_{3}(x,y,z)$,

where $L_{1}=x-x_{*}-x_{*}ln\frac{x}{x_{*}}$,  $L_{2}=y-y_{*}-y_{*}ln\frac{y}{y_{*}}$,  $L_{3}=z-z_{*}-z_{*}ln\frac{z}{z_{*}}$.

It is to be shown that $L$ is a Lyapunav function and $L$ vanishes at $E_{*}$ and it is positive for all $x,y,z>0$. Hence $E_{*}$ represents its global minimum. We have

[TABLE]

[TABLE]

[TABLE]

We consider

[TABLE]

where V=\Big{(}(x-x_{*}),(y-y_{*}),(z-z_{*})\Big{)}^{T} and $Q$ is symmetric quadratic form given by

Q=\left(\begin{array}[]{ccc}A&H&G\\ H&B&F\\ G&F&C\end{array}\right)

with the entries that are functions only of the variable $x$ and

$A=b_{1}-\frac{c_{1}(y_{*}+pz_{*})(xx_{*}-k_{1})}{(x_{*}^{2}+k_{1})(x^{2}+k_{1})}$, $B=\frac{c_{2}}{x+k_{2}}$, $C=\frac{c_{3}}{x+k_{2}}$, $F=\frac{c_{2}+c_{3}}{2(x+k_{2})}$,

$H=\frac{1}{2}[\frac{c_{1}x}{x^{2}+k_{1}}-\frac{c_{2}(y_{*}+z_{*})}{(x_{*}+k_{2})(x+k_{2})}]$, $G=\frac{1}{2}[\frac{pc_{1}x}{x^{2}+k_{1}}-\frac{c_{3}(y_{*}+z_{*})}{(x_{*}+k_{2})(x+k_{2})}]$.

If the matrix $Q$ is positive definite, then $\frac{dL}{dt}<0$. So, all the principal minors of $Q$, namely, $P_{1}\equiv A$,  $P_{2}\equiv AB-H^{2}$,  $P_{3}\equiv ABC+2FGH-AF^{2}-BG^{2}-CH^{2}=C(AB-H^{2})+G(FH-BG)+F(GH-AF)$ to be positive , i.e. , $P_{1}>0,\ P_{2}>0,\ P_{3}>0$.

∎

6 Persistence

If a compact set $D\subset\Omega=\{(x,y,z):x>0,y>0,z>0\}$ exists such that all solution of the system (2.1) eventually enter and remain in $D$, the system is called persistent.

Proposition 4**.**

*The system (2.1) is persistent if

$\frac{a_{2}c_{3}}{a_{3}}>(c_{2}+k_{2}\theta)$, 2. 2.

$a_{3}c_{2}>a_{2}c_{3}$, 3. 3.

$a_{3}c_{2}>a_{2}c_{3},\ a_{1}>b_{1}x_{5}+\frac{a_{2}c_{1}x_{5}(k_{2}+x_{5})}{c_{2}(k_{1}+x_{5}^{2})}$, 4. 4.

\frac{a_{2}c_{3}}{a_{3}}>\Big{(}c_{2}+(k_{2}+x_{6})\theta\Big{)},\ a_{1}c_{3}(k_{1}+x_{6}^{2})>x_{6}\Big{(}a_{3}c_{1}p(k_{2}+x_{6})+b_{1}c_{3}(k_{1}+x_{6}^{2})\Big{)}.

Proof.

Let us consider the method of average Lyapunav function, see [13] , considering a function of the form $V(x,y,z)=x^{\gamma_{1}}y^{\gamma_{2}}z^{\gamma_{3}}$, where $\gamma_{i}=1,2,3$ are positive constant to be determined. We define

[TABLE]

We are to prove that this function is positive at each boundary equilibrium. We have $\Pi(0,0,0)=\gamma_{1}a_{1}+\gamma_{2}a_{2}+\gamma_{3}a_{3}>0$ and $\Pi(\frac{a_{1}}{b_{1}},0,0)=\gamma_{2}a_{2}+\gamma_{3}a_{3}>0$.

Here $\Pi(0,0,\frac{a_{3}k_{2}}{c_{3}})=a_{1}\gamma_{1}+\frac{a_{2}c_{3}-a_{3}(c_{2}+k_{2}\theta)}{c_{3}}\gamma_{2}>0$ follows by condition 1.

With the condition 2, we have

$\Pi(0,\frac{a_{2}k_{2}}{c_{2}},0)=a_{1}\gamma_{1}+\frac{(a_{3}c_{2}-a_{2}c_{3}+a_{2}k_{2}\theta)}{c_{2}}\gamma_{3}>0$. Also

[TABLE]

We find

[TABLE]

[TABLE]

Hence a suitable choice of $\gamma_{1},\gamma_{2},\gamma_{3}$ is required to ensure $\Pi>0$ at the boundary equilibria. Hence $V$ is an average Lyapunav function and thus the system (2.1) is persistent. ∎

7 Numerical simulation

Analytical studies can never be completed without numerical verification of the derived results. In this section, we present computer simulations of some solutions of the system (2.1). Beside verification of our analytical findings, these numerical simulations are very important from practical point of view. We use four different set of numerical values for support of analytical results mentioned in Table 2.

8 Conclusions and comments

In this paper, we have proposed and analyzed an eco-epidemiological model where only the predator population is infected by an infectious disease. Here we have considered a modified Leslie-Gower and Holling type-III predator-prey model. We have divided the predator population into two sub classes: susceptible and infected. Then we study the behaviour of the system at various equilibrium points and their stability. The conditions for existence and stability of all the equilibria of the system have been given. The system (2.1) has eight equilibrium points: one trivial equilibrium $E_{0}$, three axial equilibrium points $E_{1}$, $E_{2}$, $E_{4}$, three planar equilibrium points $E_{3}$, $E_{5}$, $E_{6}$ and one coexistence equilibrium $E_{*}$. For our model: $E_{i}$, $i=0,1,2,4$ exist and are unstable for all times. $E_{3}$ exists if $a_{3}k_{2}\theta<a_{2}c_{3}-a_{3}c_{2}<a_{2}k_{2}\theta$ and $c_{2}>c_{3}$ but unstable. The equilibrium point $E_{5}$ is locally asymptotically stable if $b_{1}x_{5}+\frac{c_{1}x_{5}y_{5}}{x_{5}^{2}+k_{1}}>\frac{2c_{1}x_{5}^{3}y_{5}}{(x_{5}^{2}+k_{1})^{2}}$, $\frac{c_{3}y_{5}}{x_{5}+k_{2}}>\theta y_{5}+a_{3}$. Also $E_{6}$ is locally asymptotically stable if $b_{1}x_{6}+\frac{pc_{1}x_{6}z_{6}}{x_{6}^{2}+k_{1}}>\frac{2pc_{1}x_{6}^{3}z_{6}}{(x_{6}^{2}+k_{1})^{2}}$, $\frac{c_{2}z_{6}}{x_{6}+k_{2}}+\theta z_{6}>a_{2}$. The coexistence equilibrium point $E_{*}$ is locally as well as globally asymptotically stable under some conditions. Persistence of the system is also shown.

At last, we conclude that our eco-epidemic predator–prey model with infected predator exhibits very interesting dynamics. Here we have assumed Holling type III response mechanism for predation. So, we can refine the model considering other type of functional response. We can also consider the disease infection in the prey population, which can give us a very rich dynamics. There must be some time lag, called gestation delay. So, as part of future work to improve the model we can incorporate the gestation delay in our model to make it more realistic.