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

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 equilibriumpoints are calculatedwith 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.


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][6][7][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][10][11][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][14][15].Other researchers were interested in the study of disease within the predator population only [16][17][18][19].There are also some studies where both the prey and predator populations are infected with the disease [20][21][22][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][25][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 Section 2, the mathematical model is formulated

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 time  given by () and infected with density at time  represented by (), while the predator consists of only one compartment with density at time  denoted by ().
(2) The prey population grows logistically with intrinsic growth rate  > 0 and carrying capacity  > 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 rate  ∈ (0, 1) that leaves (1 − ) 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  > 0 and maximum attack rates  1 > 0 and  2 > 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  1 ∈ (0, 1) and  2 ∈ (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 rate  1 > 0, while in the absence of prey, the predator decays exponentially with natural death rate  2 > 0.
(7) Only the predator population is assumed to be harvested with the Michael-Mentence harvesting function [30], where  > 0 represents hunting effort,  > 0 is the catchability coefficient of the predator, and   ,  = 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.
where (0) ≥ 0, (0) ≥ 0, and (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: System (1) reduces to the following dimensionless system: According 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 + = {(, , ) ∈ R 3 : (0) ≥ 0, (0) ≥ 0; (0) ≥ 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.

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  0 = (0, 0, 0) and the axial equilibrium point  1 = (1, 0, 0) always exist.The first planar  2 = (, 0, ), where  represents a unique positive root of the following third-order polynomial equation: where from now onward  = (1 − ).Straightforward computation shows that  2 exists uniquely in the interior of positive quadrant of  −plane provided that the following sufficient conditions hold: The second planar equilibrium point  3 = (ŝ, î, 0), where exists uniquely in the interior of positive quadrant of  −plane provided that Finally the coexistence equilibrium point is denoted by  4 = ( * ,  * ,  * ) where while  * is a unique positive root of the following third-order polynomial equation: Here with Clearly (11c) has a unique positive root provided that one set of the following sets of conditions holds.
The Jacobian matrix at axial equilibrium point  1 = (1, 0, 0) is written as Then the eigenvalues of ( 1 ) are given by  11 = −1 < 0,  12 =  2 −  4 and  13 =  5 /( 1 + ) −  7 / 8 −  9 .Accordingly the axial point  1 is locally asymptotically stable if and only if Now the Jacobian matrix at first planar equilibrium point  2 = (, 0, ) can be written Clearly The eigenvalues of ( 2 ) can be determined as follows: where Therefore the eigenvalues of ( 2 ) are easily determined by Straightforward computation shows that these eigenvalues have negative real parts and hence the first planar equilibrium point  2 is locally asymptotically stable if and only if the following sufficient conditions hold: Similarly the Jacobian matrix at second planar equilibrium point  3 = (ŝ, î, 0) can be written Then the eigenvalues of ( 3 ) are given by where Clearly the first two eigenvalues resulting from second term in ( 27) have negative real parts, while the third eigenvalues in the  −direction which are written as have negative real part and then the second planar equilibrium point is locally asymptotically stable if and only if Finally the Jacobian matrix at coexistence equilibrium point  4 = ( * ,  * ,  * ) can be written in the form where Here  * 1 =  1 +  * +  * and  * 2 =  8 +  * .The characteristic equation of ( 4 ) can be written in the form Therefore straightforward computation shows that while with Keeping the above in view, it is well known that from Routh-Hurwitz criterion (33) has three roots with negative real parts provided that  1 > 0,  3 > 0 and û > 0. Consequently the following theorem can be proved easily.
Theorem 2. The coexistence equilibrium point  4 of system ( 3) is locally asymptotically stable in the interior of positive cone provided that the following sufficient conditions hold: 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  > 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: Note that it is easy to verify that these two subsystems have unique positive equilibrium points in the interior of their positive domains Σ 1 = {(, ) ∈ R 2 : (0) ≥ 0, (0) ≥ 0} and Σ 2 = {(, ) ∈ R 2 : (0) ≥ 0, (0) ≥ 0}, respectively.These two positive equilibrium points are given by  1 = (ŝ, î) and  2 = (, ) for subsystems ( 40) and ( 41), respectively.In fact they coincide with the predator free equilibrium point  3 and disease free equilibrium point  2 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  1 function (, ), called the Dulac function, such that the expression (/)( 11 ) + (/)( 12 ) 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 (, ) = 1/ we obtain that (/)( 11 ) + (/)( 12 ) = −1/ < 0. Therefore subsystem (40) has no periodic dynamic lying entirely in interior Σ 1 and hence  1 ≡  3 is a globally asymptotically stable in the interior of Σ 1 .Similarly by choosing (, ) = 1/ we can show that the positive equilibrium point of the second subsystem (41) given by  2 ≡  2 is globally asymptotically stable in the interior of the Σ 2 provided that one of the following two conditions holds.
Journal of Applied Mathematics Now the proof follows if   / > 0 for all the boundary equilibrium points, for suitable choice of constants  > 0,  > 0, and  > 0.

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 5. Assume that the disease free equilibrium point  2 = (, 0, ) is locally asymptotically stable in R 3 + .Then, it is a globally asymptotically stable, provided that the following conditions hold: where  1 =  1 + .
Proof.Consider the function where   ,  = 1, 2, 3, are positive constants to be determined later on.Clearly  2 (, , ) > 0 is a continuously differentiable real valued function for all (, , ) ∈ R 3 + with (, , ) ̸ = (, 0, ) and  2 (, 0, ) = 0. Moreover we have that Further simplification gives Here  2 =  8 +  and  2 =  8 + .Now rearranging the terms of the last equation further gives So by choosing the positive constants as Accordingly, using the given conditions (52a)-(52d) we obtain Clearly  2 / < 0 is negative definite, and hence the disease free equilibrium point E 2 is globally asymptotically stable under the given conditions and hence the proof is complete.

Journal of Applied Mathematics
Here ,  1 , and  2 are given above.Now straightforward computations give By choosing the positive constants as then we obtain that Accordingly, using the given conditions (60a) and (60b) we obtain Clearly  3 / ≤ 0, which means it is negative semidefinite, and hence the predator free equilibrium point E 3 is globally stable (but not asymptotically stable) under the given conditions.Moreover, since system (3) has the maximum invariant set for  3 / = 0 if and only if conditions (60a)-(60b) hold and (, , ) = (ŝ, î, 0), by Lyapunov-Lasalle's theorem, all the solutions starting in R 3 + approach the singleton set { 3 }, which is the positively invariant subset of the set where  3 / = 0. Hence E 3 becomes attracting too; hence it is globally asymptotically stable and that completes the proof.Theorem 7. Assume that the coexistence equilibrium point  4 = ( * ,  * ,  * ) of system ( 3) is locally asymptotically stable in R 3 + .Then it is globally asymptotically stable, provided that the following conditions hold: Here Λ  ,  = 1, 2, 3,  * 1 , and  * 2 are given in the proof.

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.It is observed that, for the set of data (74), system (3) has a globally asymptotically stable positive equilibrium point  4 = (0.45, 0.28, 0.17) as shown in Figure 1.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.
It is observed that for the set of data (74) with  1 ≥ 0.47 the trajectory of system (3) approaches asymptotically predator free equilibrium point, while it approaches disease free point when 0.15 <  1 ≤ 0.25; see Figure 2. Further the trajectory of system (3) approaches periodic dynamics in the  −plane for data (74) with  1 ≤ 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  2 ≥ 0.88 and 0.1 <  2 ≤ 0.17, the solution of system (3) approaches asymptotically  3 = (ŝ, î, 0) as shown in Figures 4(a)-4(b), while it approaches asymptotically  1 = (1, 0, 0) when  2 ≤ 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  3 too with decreasing the infection rate  2 , that depends on our hypothetical set of data in which we assumed that the conversion rate of predator from susceptible ( 5 ) is less than that from infected species ( 6 ) and once ( 5 ) enter the range  5 ≥ 0.47 the system approaches disease free point  2 = (, 0, ) as shown in the typical figure given by Figure 4(d).
Moreover it is observed that varying the parameter  3 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  3 .Now for the parameters values given by (74) with 0.2 ≤  4 < 0.5 and 0.5 ≤  4 the trajectory of system (3) approaches asymptotically  3 = (ŝ, î, 0) and  1 = (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 ≤  4 < 0.5, while the system still is stable at positive equilibrium point otherwise.

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 of R 3 + ; rather it either persists at the positive  infected prey represented by  5 and  6 , respectively.Finally decreasing the infection rate further ( 2 ≤ 0.1) leads to approaches to axial equilibrium point  1 .

Susceptible prey Infected prey Predator
(5) It is observed that varying the parameter  3 , 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  3 .
(6) Increasing the death rate of the infected prey ( 4 ) so that 0.2 ≤  4 < 0.5 and 0.5 ≤  4 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 above  5 ≥ 0.48 or decreasing it below  5 ≤ 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  8 , which stands for the half saturation constant of harvesting in Michael-Mentence harvesting function, as that obtained in case of increasing or decreasing  5 .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  5 too.
(8) Finally increasing the parameter  7 , 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  7 below 0.15.Similar effect on the dynamical behavior of system (3) is obtained in case of increasing or decreasing the refuge rate  and predator death rate  9 as that effect obtained with varying  7 .
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 parameter  and harvesting represented by parameters  7 and  8 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.

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.

Figure 3 :
Figure 3: System (3) approaches asymptotically to periodic dynamics using data (74) with  1 = 0.1.(a) Periodic attractor in the  −plane.(b) Trajectories of the three species approach to periodic in the  −plane.