# A Study of a Diseased Prey-Predator Model with Refuge in Prey and Harvesting from Predator

Authors: Ahmed Sami Abdulghafour; Raid Kamel Naji
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

In this paper, a mathematical model of a prey-predator system with infectious disease in the prey population is proposed and studied. It is assumed that there is a constant refuge in prey as a defensive property against predation and harvesting from the predator. The proposed mathematical model is consisting of three first-order nonlinear ordinary differential equations, which describe the interaction among the healthy prey, infected prey, and predator. The existence, uniqueness, and boundedness of the system’ solution are investigated. The system's equilibrium points are calculated with studying their local and global stability. The persistence conditions of the proposed system are established. Finally the obtained analytical results are justified by a numerical simulation.


---

## Body

## 1. Introduction

The interaction between the prey and predator species has a long history since Lotka-Volterra model; see [1]. Similarly, the interaction of the susceptible–infected–recovered population is an interesting subject of research work since the pioneering work of Kermack and McKendrick [2]. The dynamics of disease within the ecological systems now becomes an important subject of research. In fact Anderson and May [3] were the first who combined these two systems, while Chattopadhyay and Arino [4] were the first who used the term “eco-epidemiology” for such type of models. On one hand several studies of prey-predator dynamics have been done within the last decades taking into account the effects of variety of biological factors; see, for example, [5–8] and the references therein. On the other hand variety of mathematical models have been proposed and studied in the field of epidemiology taking into account different types of incidence rates and disease; see, for example, [9–12] and the references therein.The existence of disease in the prey population, predator population, or both is real situation in the ecological species. It has been observed that this type of incidence occurs through infection by some viral disease, bacterial disease, or parasite disease. Many of these studies were focused on the study of disease in prey population only [13–15]. Other researchers were interested in the study of disease within the predator population only [16–19]. There are also some studies where both the prey and predator populations are infected with the disease [20–23].It is well known that the harvesting of the species is necessary for the coexistence of the species and hence it took a lot of interest from the researchers in their proposed ecological models. Different types of harvesting have been proposed and studied including constant harvesting, density dependent proportional harvesting, and nonlinear harvesting [8, 24–26]. The refuge of prey species is also a biological factor necessary for the coexistence of the species and hence it is another factor of great interest as defensive properties of the prey against the predation; see [8, 27, 28].Keeping the above in view, in this paper we propose and study an eco-epidemiological prey-predator model involving a prey refuge and harvesting from the predator. It is assumed that the disease exists only in prey and it will not transfer to predator through the feeding process. The paper is organized as follows. In Section2, the mathematical model is formulated and its dimensionless variables and parameters are determined; moreover the existence, uniqueness, and uniform boundedness of its solution are discussed. In Section 3, the local stability of all possible equilibrium points and persistence of the system are studied. Section 4 deals with the global stability analysis of the system using suitable Lyapunov functions. However Section 5 provides a numerical simulation of the proposed system for suitable chose of parameters values. Finally Section 6 gives some conclusions and discussed the obtained results.


## 2. Mathematical Model

In this section, a prey-predator system involving infected disease in prey is proposed for study. It is assumed that there is a harvesting effort applied on the predator individuals only. Accordingly, the following hypotheses are adopted to formulate the mathematical model.(1)
The prey population is divided into two compartments, susceptible with density at timeT given by S(T) and infected with density at time T represented by I(T), while the predator consists of only one compartment with density at time T denoted by Y(T).(2)
The prey population grows logistically with intrinsic growth rater>0 and carrying capacity K>0; it is assumed that the infected cannot reproduce; rather it competes with the susceptible individuals for food and space.(3)
There is a type of protection of prey population from the predation by predator, represented by a constant prey’s refuge ratem∈(0,1) that leaves 1-mS of prey available to be hunted by the predator.(4)
The susceptible prey population becomes infected by contact with the infected prey according to the simple mass action kinetics withβ>0 as the rate of infection.(5)
The predator population consumes both the prey populations according to modified Holling type II functional response for the predation [29] with half saturation constant b>0 and maximum attack rates a1>0 and a2>0 for susceptible prey and infected prey, respectively. Since there is a vulnerability of infected prey relative to susceptible prey the vulnerability constant rate θ>0 is used in the functional response. Moreover the constants e1∈(0,1) and e2∈(0,1) are the conversion rates from susceptible and infected preys to predator, respectively.(6)
The disease causes a death in the infected population represented by diseased death rated1>0, while in the absence of prey, the predator decays exponentially with natural death rate d2>0.(7)
Only the predator population is assumed to be harvested with the Michael-Mentence harvesting function [30], where E>0 represents hunting effort, c>0 is the catchability coefficient of the predator, and li, i=1,2, are suitable positive constants.According to the above set of hypotheses the dynamics of an eco-epidemic prey-predator model with refuge in prey and harvesting in predator can be described in the following set of first-order nonlinear differential equations.(1)dSdT=rS1-S+IK-a11-mSYb+1-mS+θ1-mI-βSIdIdT=βSI-a21-mIYb+1-mS+θ1-mI-d1IdYdT=e1a11-mS+e2a21-mIb+1-mS+θ1-mIY-cEYl1E+l2Y-d2Ywhere S(0)≥0, I(0)≥0, and Y(0)≥0.Now in order to reduce the number of parameters and specify the control set of parameters the following dimensionless variables and parameters are used:(2)t=rT,x=XK,i=IK,y=a1YrK,w1=bK,w2=βKr,w3=a2a1,w4=d1r,w5=e1a1r,w6=e2a2r,w7=a1l1Erl2K,w8=ca1Er2l2K,w9=d2rSystem (1) reduces to the following dimensionless system:(3)dsdt=s1-s+i-1-myw1+1-ms+θ1-mi-w2i=sf1s,i,ydidt=iw2s-w31-myw1+1-ms+θ1-mi-w4=if2s,i,ydydt=yw51-ms+w61-miw1+1-ms+θ1-mi-w7w8+y-w9=yf3s,i,yAccording to the above dimensionless form, it is easy to verify that the differential equations are continuous and have continuous partial derivatives on the domain R+3=s,i,y∈R3:s0≥0,i0≥0;y0≥0 and hence they are Lipschitzian functions. Therefore system (3) has a unique solution. Furthermore the solution of system (3) that initials in R+3 is uniformly bounded as shown in the following theorem.Theorem 1.
All solutions of system (3) are uniformly bounded.Proof.
According to the prey population which is given in the first equation of system (3), it is observed that(4)dsdt≤s1-sThen by solving this differential inequality, we obtain that s(t)≤s(0)/e-t1-s0+s(0), and then for t→∞, we have st≤1. Now assume that Ω=s+i+y; then (5)dΩdt=dsdt+didt+dydt→yieldsdΩdt≤2-μΩwhere μ=min⁡1,w4,w9. Therefore by solving the last differential inequality we obtain that(6)Ωt≤Ω0e-μt-2μe-μt-1Hence for t→∞, we get that Ωt≤2/μ. Thus all the solutions of system (3) are uniformly bounded.


## 3. The Local Stability and Persistence

In the following the existence and local stability of all possible equilibrium points are investigated and then the persistence conditions of the system are established. Obviously system (3) has at most five nonnegative equilibrium points. The vanishing equilibrium point E0=(0,0,0) and the axial equilibrium point E1=(1,0,0) always exist. The first planar E2=(s-,0,y-), where s- represents a unique positive root of the following third-order polynomial equation:(7a)A1s3+A2s2+A3s+A4=0where A1=w9-w5R2, A2=(w9-w5)(w1-R)+w9w1R, A3=(w5w8-w7)R2+w1w5-2w1w9-w8w9R+w12w9, A4=-w1w1w9+w7+w8w9R<0, while (7b)y-=1-s-w1+Rs-Rfrom now onward R=(1-m). Straightforward computation shows that E2 exists uniquely in the interior of positive quadrant of sy-plane provided that the following sufficient conditions hold:(8a)s-<1(8b)w9>w5;A2>0orw9>w5;A3<0The second planar equilibrium point E3=(s^,i^,0), where(9)s^=w4w2;i^=w2-w4w21+w2exists uniquely in the interior of positive quadrant of si-plane provided that(10)w2>w4Finally the coexistence equilibrium point is denoted by E4=(s∗,i∗,y∗) where(11a)i∗=w3+w4-w2+w3s∗w31+w2(11b)y∗=w2s∗-w4w1w31+w2+w3+w4θR+w31+w2R-w2+w3θRsw321+w2Rwhile s∗ is a unique positive root of the following third-order polynomial equation:(11c)D1s3+D2s2+D3s+D4=0Here (12)D1=σ3-w2+w3σ6w2σ5D2=w2σ4-w4σ5σ3-w2w2+w3σ4σ6+w2w3+w4+w4w2+w3σ5σ6-w2σ5σ8D3=σ1-σ2w2+w3-w4σ3σ4+w2w3+w4+w4w2+w3σ4σ6-w4w3+w4σ5σ6-w2σ4σ8+w4σ5σ8D4=σ2w3+w4-σ7-w4w3+w4σ4σ6+w4σ4σ8with (13)σ1=w5w8-w7+w8w9w331+w22R2σ2=w6w8-θw7+w8w9w321+w2R2σ3=w5-w9w31+w2R2σ4=w1w31+w2+w3+w4θRσ5=w31+w2R-w2+w3θRσ6=w6-w9θRσ7=w1w7+w8w9w331+w22Rσ8=w1w9w31+w2RClearly (11c) has a unique positive root provided that one set of the following sets of conditions holds.(14a)D1>0;D2>0;D4<0(14b)D1>0;D3<0;D4<0(14c)D1<0;D2<0;D4>0(14d)D1<0;D3>0;D4>0Consequently, the coexistence equilibrium point E4=(s∗,i∗,y∗) exists uniquely in the interior of R+3 provided that in addition to any one of conditions (14a), (14b), (14c), and (14d) the following sufficient conditions hold:(15a)s∗<w3+w4w2+w3(15b)w2s∗-w4σ4+σ5s>0Now in order to study the local stability of these equilibrium points, the Jacobian matrix of system (3) at the point E=(s,i,y) is computed as follows: (16)JE=Cij3×3where (17)C11=f1s,i,y+sR2yN12-1;C12=sθR2yN12-1+w2;C13=s-RN1;C21=iw2+w3R2yN12;C22=f2s,i,y+iθw3R2yN12;C23=i-w3RN1;C31=yw1w5R+w5θ-w6R2iN12;C32=yw1w6R+w6-θw5R2sN12;C33=yw7N22+f3s,i,y.Here N1=w1+Rs+θRi and N2=w8+y. Consequently the Jacobian matrix at vanishing equilibrium point E0=(0,0,0) is written as (18)JE0=1000-w4000-w7w8+w9Obviously the eigenvalues of JE0 are given by λ01=1>0, λ02=-w4<0, λ03=-w7/w8+w9<0. Therefore E0 is a saddle point.The Jacobian matrix at axial equilibrium pointE1=(1,0,0) is written as(19)
J
E
1
=
-
1
-
1
+
w
2
-
R
w
1
+
R
0
w
2
-
w
4
0
0
0
w
5
R
w
1
+
R
-
w
7
w
8
-
w
9Then the eigenvalues of JE1 are given by λ11=-1<0, λ12=w2-w4 and λ13=w5R/w1+R-w7/w8-w9. Accordingly the axial point E1 is locally asymptotically stable if and only if(20a)w2<w4(20b)w5Rw1+R<w7w8+w9Now the Jacobian matrix at first planar equilibrium point E2=(s-,0,y-) can be written(21)
J
E
2
=
R
2
s
-
y
-
w
1
+
R
s
-
2
-
s
-
θ
R
2
s
-
y
-
w
1
+
R
s
-
2
-
1
+
w
2
s
-
-
R
s
-
w
1
+
R
s
-
0
w
2
s
-
-
w
3
R
y
-
w
1
+
R
s
-
-
w
4
0
w
1
w
5
R
y
-
w
1
+
R
s
-
2
w
1
w
6
R
y
-
+
w
6
-
θ
w
5
R
2
s
-
y
-
w
1
+
R
s
-
2
w
7
y
-
w
8
+
y
-
2Clearly The eigenvalues of JE2 can be determined as follows:(22)w2s--w3Ry-w1+Rs--w4-λλ2-T2λ+D2λ=0where(23)T2=R2s-y-w1+Rs-2-s-+w7y-w8+y-2D2=R2s-y-w1+Rs-2-s-w7y-w8+y-2+Rs-w1+Rs-w1w5Ry-w1+Rs-2Therefore the eigenvalues of J(E2) are easily determined by(24)λ21=T22-12T22-4D2λ22=w2s--w3Ry-w1+Rs--w4λ23=T22+12T22-4D2Straightforward computation shows that these eigenvalues have negative real parts and hence the first planar equilibrium point E2 is locally asymptotically stable if and only if the following sufficient conditions hold:(25a)w2s-<w3Ry-w1+Rs-+w4(25b)R2s-y-w1+Rs-2+w7y-w8+y-2<s-(25c)w7R2y-w1+Rs-+w1w5R2w8+y-2>w7w1+Rs-3Similarly the Jacobian matrix at second planar equilibrium point E3=(s^,i^,0) can be written (26)
J
E
3
=
-
s
^
-
s
^
1
+
w
2
-
R
s
^
w
1
+
R
s
^
+
θ
R
i
^
w
2
i
^
0
-
w
3
R
i
^
w
1
+
R
s
^
+
θ
R
i
^
0
0
w
5
R
s
^
+
w
6
R
i
^
w
1
+
R
s
^
+
θ
R
i
^
-
w
7
w
8
-
w
9Then the eigenvalues of JE3 are given by (27)w5Rs^+w6Ri^w1+Rs^+θRi^-w7w8-w9-λλ2-T3λ+D3=0where(28)T3=-s^<0D3=w21+w2s^i^>0Clearly the first two eigenvalues resulting from second term in (27) have negative real parts, while the third eigenvalues in the y-direction which are written as(29)λ33=w5Rs^+w6Ri^w1+Rs^+θRi^-w7w8-w9have negative real part and then the second planar equilibrium point is locally asymptotically stable if and only if(30)w5Rs^+w6Ri^w1+Rs^+θRi^<w7w8+w9Finally the Jacobian matrix at coexistence equilibrium point E4=(s∗,i∗,y∗) can be written in the form(31)JE4=aiji,j=1,2,3where(32)a11=R2s∗y∗-N1∗2s∗N1∗2;a12=θR2s∗y∗-1+w2N1∗2s∗N1∗2;a13=-Rs∗N1∗<0;a21=w2i∗+w3R2i∗y∗N1∗2>0;a22=θw3R2i∗y∗N1∗2>0;a23=-w3Ri∗N1∗<0;a31=w1w5Ry∗+w5θ-w6R2i∗y∗N1∗2;a32=w1w6Ry∗+w6-θw5R2s∗y∗N1∗2;a33=w7y∗N2∗2>0.Here N1∗=w1+Rs∗+θRi∗ and N2∗=w8+y∗. The characteristic equation of JE4 can be written in the form(33)λ3+A1λ2+A2λ+A3=0where(34)A1=-a11+a22+a33A2=a11a22-a12a21+a11a33-a13a31+a22a33-a23a32A3=a33a12a21-a11a22-a23a12a31-a11a32-a13a21a32-a22a31while(35)∆=A1A2-A3=-a11+a22a11a22-a12a21-2a11a22a33-a11+a33a11a33-a13a31-a22+a33a23a32-a22a33Therefore straightforward computation shows that(36)A1=-1N1∗2N2∗2N2∗2s∗R2y∗-N1∗2+θw3R2N2∗2i∗y∗+w7N1∗2y∗A3=w7s∗i∗y∗N1∗2N2∗2w2+w3θR2y∗-1+w2w2N1∗2+w3R2y∗+w3Rs∗i∗y∗N1∗5w5θ-w6R2N1∗Ry∗-N1∗s∗+1+w2i∗+w1w6-w51+w2RN1∗2+R2s∗2y∗N1∗4w1w2w6N1∗+w6Rw2N1∗s∗+w3Ry∗-w5θRw2N1∗s∗+w3Ry∗while(37)∆=-M1∗+M2∗+M3∗-M4∗with(38)M1∗=1N1∗6N1∗2s∗i∗y∗R2s∗y∗-N1∗2s∗+w3θR2i∗y∗1+w2w2N1∗2+w3R2y∗-w2+w3θR2y∗M2∗=w3w7θR2s∗i∗y∗2N1∗4N2∗2N1∗2-R2y∗M3∗=s∗y∗N1∗6N2∗4-R2N2∗2s∗y∗+N1∗2N2∗2s∗-w7N1∗2y∗w7R2N1∗2y∗-w7N1∗4+w1w5R2N2∗2+w5θ-w6R3N2∗2i∗M4∗=w3R2i∗y∗2N1∗6N2∗2w3θR2N2∗2i∗+w7N2∗2w1w6N1∗+w6-θw5RN1∗s∗+θw7y∗Keeping the above in view, it is well known that from Routh-Hurwitz criterion (33) has three roots with negative real parts provided that A1>0, A3>0 and ∆>0. Consequently the following theorem can be proved easily.Theorem 2.
The coexistence equilibrium pointE4 of system (3) is locally asymptotically stable in the interior of positive cone provided that the following sufficient conditions hold:(39a)max⁡1+w2w5,θw5Rs∗w1+Rs∗<w6<w5θ(39b)max⁡1+w2w2N1∗2+w3R2y∗θw2+w3,RN1∗s∗+1+w2i∗<R2y∗<1+w2θN1∗2(39c)R2N2∗2s∗y∗+θw3R2N2∗2i∗y∗+w7N1∗2y∗<N1∗2N2∗2s∗(39d)w1w2w6N1∗+w6Rw2N1∗s∗+w3Ry∗>w5θRw2N1∗s∗+w3Ry∗(39e)w7R2N1∗2y∗+w1w5R2N2∗2+w5θ-w6R3N2∗2i∗>w7N1∗4(39f)M2∗+M3∗>M1∗+M4∗Now before we start studying the global stability of the system, the persistence of system (3), which represents biologically the coexistence of all the species for all the time while mathematically it indicates that the solution of system (3) for all t>0 does not have omega limit set on the boundary planes, is investigated in the following theorem when there are no periodic dynamics on the boundary planes.
Clearly system (3) has two subsystems. The first subsystem exists in case of absence of predator species, while the second subsystem exists in case of absence of disease. These subsystems can be written, respectively, as follows:(40)dsdt=s1-s+i-w2i=f11s,ididt=iw2s-w4=f12s,i(41)dsdt=s1-s-Ryw1+Rs=f21s,ydydt=yw5Rsw1+Rs-w7w8+y-w9=f23s,yNote that it is easy to verify that these two subsystems have unique positive equilibrium points in the interior of their positive domains Σ1=s,i∈R2:s0≥0,i0≥0 and Σ2=s,y∈R2:s0≥0,y0≥0, respectively. These two positive equilibrium points are given by e1=(s^,i^) and e2=(s-,y-) for subsystems (40) and (41), respectively. In fact they coincide with the predator free equilibrium point E3 and disease free equilibrium point E2 of system (3), respectively, and have the same existence conditions.
Therefore according to the Bendixson–Dulac theorem on dynamical system (40) “if there exists a C1 function φ(s,i), called the Dulac function, such that the expression ∂/∂s(φf11)+∂/∂i(φf12) has the same sign and is not identically zero (≠0) almost everywhere in the simply connected region of the plane, then system (40) has no nonconstant periodic dynamic lying entirely in the Σ1. Thus by choosing φs,i=1/si we obtain that ∂/∂s(φf11)+∂/∂iφf12=-1/i<0. Therefore subsystem (40) has no periodic dynamic lying entirely in interior Σ1 and hence e1≡E3 is a globally asymptotically stable in the interior ofΣ1. Similarly by choosing φs,i=1/sy we can show that the positive equilibrium point of the second subsystem (41) given by e2≡E2 is globally asymptotically stable in the interior of the Σ2 provided that one of the following two conditions holds.(42)R2w1+Rs2+w7sw8+y2>1y(43)R2w1+Rs2+w7sw8+y2<1yTheorem 3.
Assume that condition (42) or (43) holds. Then system (3) is persistent provided that(44a)w2>w4(44b)w5Rw1+R>w7w8+w9(44c)w2s->w3Ry-w1+Rs-+w4(44d)w5Rs^+w6Ri^w1+Rs^+θRi^>w7w8+w9Proof.
Consider the following functionφs,i,y=sαiβyγ, where α, β, and γ are positive constants; clearly φs,i,y>0 for all s,i,y∈IntR+3 and φs,i,y→0 when s,iory→0. Furthermore, it is clear that (45)φ′φ=α1-s+i-Ryw1+Rs+θRi-w2i+βw2s-w3Ryw1+Rs+θRi-w4+γw5Rs+w6Riw1+Rs+θRi-w7w8+y-w9Now the proof follows if φ′/φ>0 for all the boundary equilibrium points, for suitable choice of constants α>0, β>0, and γ>0. (46)φ′φE0=α-w4β-w7w8+w9γφ′φE1=w2-w4β+w5Rw1+R-w7w8-w9γφ′φE2=w2s--w3Ry-w1+Rs--w4βφ′φE3=w5Rs^+w6Ri^w1+Rs^+θRi^-w7w8-w9γClearly φ′/φE0>0 for suitable choice of positive constant α sufficiently large with respect to the constants β>0 and γ>0, while φ′/φE1 is positive under conditions (44a) and (44b). Moreover φ′/φE2>0 under condition (44c), while φ′/φE3 is positive too under condition (44d). Hence the proof is complete.


## 4. Global Stability

In this section the global stability of each equilibrium point of system (3) is studied using suitable Lyapunov function as given in the following theorems.Theorem 4.
Assume that the axial equilibrium pointE1=1,0,0=(s0,0,0) is locally asymptotically stable in R+3. Then, it is globally asymptotically stable, provided that the following conditions are outstanding:(47a)w6w3w2<w51+w2<w6w3w4(47b)w9>w5Rw1Proof.
Recognize the following function:(48)V1=c1∫s0sτ-s0τdτ+c2i+c3ywhere ci, i=1,2,3, are positive constants to be determined later on. Clearly V1s,i,y>0 is a continuously differentiable real valued function for all s,i,y∈R+3 with s,i,y≠s0,0,0 and V1s0,0,0=0. Moreover we have that(49)dV1dt=c1s-s0sdsdt+c2didt+c3dydt=-c1s-s02-c1RsyN1+c1Rs0yN1-c11+w2si+c11+w2s0i+c2w2si-c2w3RiyN1-c2w4i+c3w5Rs+w6RiN1y-c3w7yw8+y-c3w9yHere R=1-m, N1=w1+Rs+θRi. Further simplification gives that(50)dV1dt≤-c1s-s02-c1-c3w5RsyN1-c11+w2-c2w2si-c2w3-c3w6RiyN1-c2w4-c11+w2s0i-c3w9-c1Rs0w1ySo by choosing the positive constants ci, i=1,2,3, as given below, we get(51)c1=w5;c2=w6w3;c3=1dV1dt≤-w5s-12-w51+w2-w6w3w2si-w6w3w4-w51+w2i-w9-w5Rw1yClearly conditions (47a)-(47b) guarantee that dV1/dt<0; hence dV1/dt is negative definite, and hence the axial equilibrium point E1 is globally asymptotically stable and the proof is complete.Theorem 5.
Assume that the disease free equilibrium pointE2=s-,0,y- is locally asymptotically stable in R+3. Then, it is a globally asymptotically stable, provided that the following conditions hold:(52a)w4w11+w2N-1>w4R2θy-+w1w21+w2N-1s-(52b)w5θRs-N-12<w6<w1w3w51+w2s-+w4w5θRs-w4N-1(52c)w1N-1>R2y-(52d)w7N-1w1w5w8N-2y-y-2<1-R2y-w1N-1s-s-2where N-1=w1+Rs-.Proof.
Consider the function(53)V2=c-1∫s-sτ-s-τdτ+c-2i+c-3∫y-yτ-y-τdτwhere c-i, i=1,2,3, are positive constants to be determined later on. Clearly V2s,i,y>0 is a continuously differentiable real valued function for all s,i,y∈R+3 with s,i,y≠s-,0,y- and V2s-,0,y-=0. Moreover we have that(54)dV2dt=c-1s-s-1-s+i-RyN1-w2i+c-2iw2s-w3RyN1-w4+c-3y-y-w5Rs+w6RiN1-w7w8+y-w9Further simplification gives (55)dV2dt=c-1s-s--s-s--1+w2i-RN1y-y-+R2y-N1N-1s-s-+R2θy-N1N-1i+w2c-2is-w3Rc-2N1iy-w4c-2i+c-3y-y-w5RsN1-w5Rs-N-1+w6RiN1+w7N2N-2y-y-Here N2=w8+y and N-2=w8+y-.
Now rearranging the terms of the last equation further gives(56)dV2dt≤-c-11-R2y-w1N-1s-s-2-RN1c-1-c-3w1w5N-1s-s-y-y-+c-3w7w8N-2y-y-2-c-11+w2-R2θy-w1N-1-c-2w2si-RN1c-2w3+c-3w5θRs-N-1-w6iy-RN1c-1Rθs-y-N-1-c-3w5Rθs-y-N-1-w6y-i-c-2w4-c-11+w2s-iSo by choosing the positive constants as(57)c-1=1,c-2=1+w2s-w4,c-3=N-1w1w5we get that (58)dV2dt≤-1-R2y-w1N-1s-s-2+w7N-1w1w5w8N-2y-y-2-Ry-w1N1w6N-12-w5θRs-2w5N-1i-w4w11+w2N-1-R2θy--w1w21+w2N-1s-w1w4N-1si-RN1w1w3w51+w2s-+w4w5θRs--w4w6N-1w1w4w5iyAccordingly, using the given conditions (52a)–(52d) we obtain (59)dV2dt≤-1-R2y-w1N-1s-s-2+w7N-1w1w5w8N-2y-y-2-Ry-w1N1w6N-12-w5θRs-2w5N-1iClearly dV2/dt<0 is negative definite, and hence the disease free equilibrium point E2 is globally asymptotically stable under the given conditions and hence the proof is complete.Theorem 6.
Assume that the free predator equilibrium pointE3=(s^,i^,0) is locally asymptotically stable in R+3. Then it is globally asymptotically stable, provided that the following condition holds:(60a)w9w5>Rw2s^+w31+w2i^w1w2(60b)w3w51+w2>w2w6Proof.
Consider the next function(61)V3=c^1∫s^sτ-s^τdτ+c^2∫i^iτ-i^τdτ+c^3ywhere c^i, i=1,2,3, are positive constants to be determined later on. Clearly V3s,i,y>0 is a continuously differentiable real valued function for all s,i,y∈R+3 with s,i,y≠s^,i^,0 and V3s^,i^,0=0. Moreover we have that(62)dV3dt=c^1s-s^1-s+i-RyN1-w2i+c^2i-i^w2s-w3RyN1-w4+c^3yw5Rs+w6RiN1-w7N2-w9Here R, N1, and N2 are given above. Now straightforward computations give (63)dV3dt≤-c^1s-s^2-c^11+w2-c^2w2i-i^s-s^-c^3w7N2y-c^3w9-c^2w3Rw1i^-c^1Rs^w1y-RN1c^1-c^3w5sy-RN1c^2w3-c^3w6iyBy choosing the positive constants as (64)c^1=1,c^2=1+w2w2,c^3=1w5then we obtain that(65)dV3dt≤-s-s^2-w9w5-Rw2s^+w31+w2i^w1w2y-RN1w3w51+w2-w2w6w2w5iyAccordingly, using the given conditions (60a) and (60b) we obtain (66)dV3dt≤-s-s^2-w9w5-Rw2s^+w31+w2i^w1w2yClearly dV3/dt≤0, which means it is negative semi-definite, and hence the predator free equilibrium point E3 is globally stable (but not asymptotically stable) under the given conditions. Moreover, since system (3) has the maximum invariant set for dV3/dt=0 if and only if conditions (60a)-(60b) hold and s,i,y=s^,i^,0, by Lyapunov-Lasalle’s theorem, all the solutions starting in R+3 approach the singleton set E3, which is the positively invariant subset of the set where dV3/dt=0. Hence E3 becomes attracting too; hence it is globally asymptotically stable and that completes the proof.Theorem 7.
Assume that the coexistence equilibrium pointE4=(s∗,i∗,y∗) of system (3) is locally asymptotically stable in R+3. Then it is globally asymptotically stable, provided that the following conditions hold:(67a)R2y∗<min⁡w1N1∗,1+w2w1N1∗θ(67b)θw5Rs∗w1+Rs∗<w6<w5w1+θRi∗Ri∗(67c)Λ1=c1∗1+w2N1N1∗-θR2y∗-c2∗w2N1N1∗+w3R2y∗>0(67d)2Λ2-Λ1>0(67e)2Λ3+Λ1N1N1∗i-i∗2+w7N2N2∗y-y∗2<2Λ2-Λ1N1N1∗s-s∗2Here Λi, i=1,2,3, c1∗, and c2∗ are given in the proof.Proof.
Consider the real valued function(68)V4=c1∗∫s∗sτ-s∗τdτ+c2∗∫i∗iτ-i∗τdτ+c3∗∫y∗yτ-y∗τdτ.Here ci∗, i=1,2,3, are positive constants to be determined. Clearly V4s,i,y>0 is a continuously differentiable real valued function for all s,i,y∈R+3 with s,i,y≠s∗,i∗,y∗ and V3s∗,i∗,y∗=0. Moreover we have that(69)dV4dt=c1∗s-s∗-s-s∗-1+w2i-i∗-RyN1+Ry∗N1∗+c2∗i-i∗w2s-s∗-w3RyN1+w3Ry∗N1∗+c3∗y-y∗N3N1-N3∗N1∗-w7N2+w7N2∗Here R, N1, and N2 are given above, while N1∗=w1+Rs∗+θRi∗, N2∗=w8+y∗, N3=w5Rs+w6Ri, and N3∗=w5Rs∗+w6Ri∗. Further simplification for the above equation leads to(70)dV4dt≤-c1∗w1N1∗-R2y∗w1N1∗s-s∗2+c2∗θw3R2y∗w1N1∗i-i∗2+c3∗w7w8N2∗y-y∗2-s-s∗i-i∗c1∗w2+1N1N1∗-θR2y∗N1N1∗-c2∗w2N1N1∗+w3R2y∗N1N1∗-s-s∗y-y∗RN1c1∗-c3∗N1∗w1w5+θw5Ri∗-w6Ri∗-i-i∗y-y∗RN1c2∗w3-c3∗N1∗w1w6-θw5Rs∗+w6Rs∗Consequently, by choosing the positive constants ci∗, i=1,2,3, as(71)c1∗=w1w5+θw5Ri∗-w6Ri∗N1∗;c2∗=w1w6-θw5Rs∗+w6Rs∗w3N1∗andc3∗=1the following is obtained:(72)dV4dt=-Λ2N1N1∗s-s∗2+Λ3N1N1∗i-i∗2+w7N2N2∗y-y∗2-Λ1N1N1∗s-s∗i-i∗where Λ2=c1∗N1N1∗-R2y∗>0 and Λ3=c2∗θw3R2y∗>0. Further simplification for the last equation gives (73)dV4dt=-2Λ2-Λ1N1N1∗s-s∗2+2Λ3+Λ1N1N1∗i-i∗2+w7N2N2∗y-y∗2-Λ12N1N1∗s-s∗+Λ12N1N1∗i-i∗2Clearly, dV4/dt<0 is negative definite function under the given conditions and hence the coexistence equilibrium point E4=(s∗,i∗,y∗) is globally asymptotically stable and this completes the proof.


## 5. Numerical Simulation

In this section the global dynamics of system (3) is studied numerically to verify the obtained analytical results in addition to specifying the control set of parameters. For the following hypothetical set of parameters, system (3) is solved numerically and the obtained trajectories are drawn in the form of phase portrait and time series.(74)w1=0.4,w2=0.5,w3=1,w4=0.1,w5=0.4w6=0.7,w7=0.2,w8=1,w9=0.1,m=0.4,θ=1It is observed that, for the set of data (74), system (3) has a globally asymptotically stable positive equilibrium point E4=0.45,0.28,0.17 as shown in Figure 1.Figure 1
Time series for the trajectories of system (3) using data (74) with different sets of initial points. (a) Trajectories of Susceptible prey. (b) Trajectories of infected prey. (c) Trajectories of predator.
(a) (b) (c)Now in order to discover the impact of varying the parameters values on the dynamics of system (3), the system is solved numerically with varying one parameter each time and then the attractors of the obtained trajectories are present in the form of figures as shown in Figures 2–7.Figure 2
Time series for the trajectory of the system (3) using data (74) with typical values of w1. (a) Trajectories of three species approach asymptotically to E3=(0.2,0.53,0) when w1=0.55. (b) Trajectories of three species approach asymptotically to E2=(0.5,0,0.41) when w1=0.2.
(a) (b)Figure 3
System (3) approaches asymptotically to periodic dynamics using data (74) with w1=0.1. (a) Periodic attractor in the sy-plane. (b) Trajectories of the three species approach to periodic in the sy-plane.
(a) (b)Figure 4
Time series for the trajectory of the system (3) using data (74) with typical values of w2. (a) Trajectories of three species approach asymptotically to E3=(0.1,0.45,0) when w2=0.95. (b) Trajectories of three species approach asymptotically to E3=(0.66,0.29,0) when w2=0.15. (c) Trajectories of three species approach asymptotically to E1=(1,0,0) when w2=0.05. (d) Trajectories of three species approach asymptotically to E2=(0.65,0,0.46) when w2=0.15 and w5=0.48.
(a) (b) (c) (d)Figure 5
Time series for the trajectory of the system (3) using data (74) with typical values of w4. (a) Trajectories of three species approach asymptotically to E3=(0.5,0.33,0) when w4=0.25. (b) Trajectories of three species approach asymptotically to E2=(1,0,0) when w4=0.55.
(a) (b)Figure 6
Time series for the trajectory of the system (3) using data (74) with typical values of w5. (a) Trajectories of three species approach asymptotically to E2=(0.57,0,0.53) when w5=0.5. (b) Trajectories of three species approach asymptotically to E3=(0.2,0.53,0) when w5=0.15.
(a) (b)Figure 7
Time series for the trajectory of the system (3) using data (74) with typical values of w7. (a) Trajectories of three species approach asymptotically to E3=(0.2,0.53,0) when w7=0.3. (b) Trajectories of three species approach asymptotically to E2=(0.45,0,0.61) when w7=0.1.
(a) (b)It is observed that for the set of data (74) with w1≥0.47 the trajectory of system (3) approaches asymptotically predator free equilibrium point, while it approaches disease free point when 0.15<w1≤0.25; see Figure 2. Further the trajectory of system (3) approaches periodic dynamics in the sy-plane for data (74) with w1≤0.15 as shown in Figure 3; however it approaches asymptotically the positive equilibrium point otherwise.On the other hand, for the data (74) with w2≥0.88 and 0.1<w2≤0.17, the solution of system (3) approaches asymptotically E3=(s^,i^,0) as shown in Figures 4(a)-4(b), while it approaches asymptotically E1=(1,0,0) when w2≤0.1 as shown in Figure 4(c). Otherwise the system still has a globally asymptotically stable positive equilibrium point. Note that although it looks confusing as the trajectory of system (3) approaches the predator free equilibrium point E3 too with decreasing the infection rate w2, that depends on our hypothetical set of data in which we assumed that the conversion rate of predator from susceptible (w5) is less than that from infected species (w6) and once (w5) enter the range w5≥0.47 the system approaches disease free point E2=(s-,0,y-) as shown in the typical figure given by Figure 4(d).Moreover it is observed that varying the parameterw3 with the rest of parameters as in (74) has a quantitative effect on the dynamics of system (3) and the solution still approaches a positive equilibrium point that depends on the value of w3. Now for the parameters values given by (74) with 0.2≤w4<0.5 and 0.5≤w4 the trajectory of system (3) approaches asymptotically E3=(s^,i^,0) and E1=(1,0,0), respectively, as given in Figure 5. Otherwise the system still is stable at positive equilibrium point.‎ There is similar behavior of the parameter θ in the range θ≥1.3 to that of system (3) using data (74) with 0.2≤w4<0.5, while the system still is stable at positive equilibrium point otherwise.For the parameters (74) with w5≥0.48 (similarly when w8≥1.5) the trajectory of system (3) approaches asymptotically E2=(s-,0,y-), while it approaches E3=(s^,i^,0) for data (74) with w5≤0.25 (similarly when w6≤0.65 or w8≤0.8) as shown in Figure 6. Otherwise the system has a globally asymptotically stable positive point.Finally for the parameters (74) with w7≥0.23 (similarly when w9≥0.13 or m≥0.49) the trajectory of system (3) approaches asymptotically E3=(s^,i^,0), while it approaches E2=(s-,0,y-) for data (74) with w7≤0.15 (similarly when w9≤0.06 or m≤0.06) as shown in Figure 7. Otherwise the system has a globally asymptotically stable positive point.


## 6. Discussion and Conclusions

The dynamics of a refuged prey-predator system, involving infectious disease in prey species and harvesting from a predator species, is mathematically simulated through a mathematical model consisting of three nonlinear ordinary differential equations of the first order. The existence, uniqueness, and boundedness of the solution of the proposed model are discussed analytically. All feasible equilibrium points are determined and then the local stability analysis for them is carried out. The persistence conditions of the system are established. Suitable Lyapunov functions are used to show the global stability of the system’s equilibrium points. Finally the proposed dynamical system is solved numerically in order to confirm the obtained analytical results and specify the control set of parameters too. It is observed that for the hypothetical set of parameters given by (74) the following results are obtained; different sets of parameters values may be used too.(1)
For the data (74), system (3) has a globally asymptotically stable positive equilibrium point in the interior of R+3.(2)
The system has no periodic dynamics lying in the interior ofR+3; rather it either persists at the positive equilibrium point or else loses its persistence and the system approaches asymptotically a specific attractor in the boundary planes.(3)
Increasing the half saturation constant represented byw1 above 0.47 leads to losing the persistence of system (3) and the system approaches asymptotically the predator free equilibrium point, while decreasing half saturation constant in the range 0.15<w1≤0.25 makes the solution approaches asymptotically the disease free equilibrium point. Finally further decreasing of w1 below the value 0.15 leads to losing the stability of the disease free point and the solution approaches periodic dynamics in the interior of sy-plane.(4)
Increasing the infection rate parameter (w2) above the value 0.88 leads to losing the persistence of system (3) too and the solution of system (3) approaches asymptotically the predator free equilibrium point. However decreasing the value of this parameter to the range 0,1<w2≤0.17 makes the system approaches asymptotically either disease free equilibrium point or predator free equilibrium point depending on the values of conversion rates from susceptible prey and infected prey represented by w5 and w6, respectively. Finally decreasing the infection rate further (w2≤0.1) leads to approaches to axial equilibrium point E1.(5)
It is observed that varying the parameterw3, which represents the ratio of the predator’s attack rate of infected prey to predator’s attack rate of susceptible prey, has a quantitative effect on the dynamics of system (3) and the system still persists at the positive equilibrium point that depends on the value of w3.(6)
Increasing the death rate of the infected prey (w4) so that 0.2≤w4<0.5 and 0.5≤w4 leads to losing the persistence of system (3) and the trajectory of system approaches asymptotically the predator free equilibrium point and axial equilibrium point, respectively. However the system still persists at the positive equilibrium point otherwise. Similar behavior of the dynamics has been obtained when the vulnerability constant rate increases above 1.3 with the rest of parameters as in (74).(7)
Now increasing the conversion rate from the susceptible prey abovew5≥0.48 or decreasing it below w5≤0.25 causes losing of persistence of the system and the trajectory approaches asymptotically the disease free equilibrium point and predator free equilibrium point, respectively. Similar dynamical behavior has been obtained when increasing or decreasing the parameter w8, which stands for the half saturation constant of harvesting in Michael-Mentence harvesting function, as that obtained in case of increasing or decreasing w5. The decreasing of the conversion rate from the infected prey below 0.65 has the same dynamical effects on system (3) as that obtained with decreasing w5 too.(8)
Finally increasing the parameterw7, which stands for the maximum harvesting rate in Michael-Mentence harvesting function, above 0.23 leads to losing the persistence of system (3) and the solution approaches asymptotically the predator free equilibrium point, while it approaches asymptotically the disease free equilibrium point with decreasing the parameter w7 below 0.15. Similar effect on the dynamical behavior of system (3) is obtained in case of increasing or decreasing the refuge rate m and predator death rate w9 as that effect obtained with varying w7.Keeping the above in view, it is easy to verify that all the analytical stability conditions are satisfied for each case in the above-mentioned point. Furthermore, the refuge represented by parameterm and harvesting represented by parameters w7 and w8 have a vital effect on the dynamical behavior of system (3) including losing the persistence and moving between the disease free equilibrium point and predator free equilibrium point.



---

*Source: 2952791-2018-12-06.xml*