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.

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, [58] 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, [912] 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 [1315]. Other researchers were interested in the study of disease within the predator population only [1619]. There are also some studies where both the prey and predator populations are infected with the disease [2023].

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, 2426]. 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 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 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 and carrying capacity ; 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 that leaves 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 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 and maximum attack rates and for susceptible prey and infected prey, respectively. Since there is a vulnerability of infected prey relative to susceptible prey the vulnerability constant rate is used in the functional response. Moreover the constants and 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 , while in the absence of prey, the predator decays exponentially with natural death rate .(7)Only the predator population is assumed to be harvested with the Michael-Mentence harvesting function [30], where represents hunting effort, is the catchability coefficient of the predator, and , , 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 , , and .

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 and hence they are Lipschitzian functions. Therefore system (3) has a unique solution. Furthermore the solution of system (3) that initials in 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 thatThen by solving this differential inequality, we obtain that , and then for , we have . Now assume that ; then where . Therefore by solving the last differential inequality we obtain thatHence for , we get that . 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 and the axial equilibrium point always exist. The first planar , where represents a unique positive root of the following third-order polynomial equation:where , , , , while from now onward . Straightforward computation shows that exists uniquely in the interior of positive quadrant of plane provided that the following sufficient conditions hold:The second planar equilibrium point , whereexists uniquely in the interior of positive quadrant of plane provided thatFinally the coexistence equilibrium point is denoted by wherewhile 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.Consequently, the coexistence equilibrium point exists uniquely in the interior of provided that in addition to any one of conditions (14a), (14b), (14c), and (14d) the following sufficient conditions hold:Now in order to study the local stability of these equilibrium points, the Jacobian matrix of system (3) at the point is computed as follows: where Here and . Consequently the Jacobian matrix at vanishing equilibrium point is written as Obviously the eigenvalues of are given by , , . Therefore is a saddle point.

The Jacobian matrix at axial equilibrium point is written asThen the eigenvalues of are given by , and . Accordingly the axial point is locally asymptotically stable if and only ifNow the Jacobian matrix at first planar equilibrium point can be writtenClearly The eigenvalues of can be determined as follows:whereTherefore the eigenvalues of are easily determined byStraightforward computation shows that these eigenvalues have negative real parts and hence the first planar equilibrium point is locally asymptotically stable if and only if the following sufficient conditions hold:Similarly the Jacobian matrix at second planar equilibrium point can be written Then the eigenvalues of are given by whereClearly the first two eigenvalues resulting from second term in (27) have negative real parts, while the third eigenvalues in the direction which are written ashave negative real part and then the second planar equilibrium point is locally asymptotically stable if and only ifFinally the Jacobian matrix at coexistence equilibrium point can be written in the formwhereHere and . The characteristic equation of can be written in the formwherewhileTherefore straightforward computation shows thatwhilewithKeeping the above in view, it is well known that from Routh-Hurwitz criterion (33) has three roots with negative real parts provided that , and . Consequently the following theorem can be proved easily.

Theorem 2. The coexistence equilibrium point 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 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 and , respectively. These two positive equilibrium points are given by and for subsystems (40) and (41), respectively. In fact they coincide with the predator free equilibrium point and disease free equilibrium point 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 function , called the Dulac function, such that the expression has the same sign and is not identically zero almost everywhere in the simply connected region of the plane, then system (40) has no nonconstant periodic dynamic lying entirely in the . Thus by choosing we obtain that . Therefore subsystem (40) has no periodic dynamic lying entirely in interior and hence is a globally asymptotically stable in the interior of. Similarly by choosing we can show that the positive equilibrium point of the second subsystem (41) given by is globally asymptotically stable in the interior of the provided that one of the following two conditions holds.

Theorem 3. Assume that condition (42) or (43) holds. Then system (3) is persistent provided that

Proof. Consider the following function , where , , and are positive constants; clearly for all and when . Furthermore, it is clear that Now the proof follows if for all the boundary equilibrium points, for suitable choice of constants , , and . Clearly for suitable choice of positive constant sufficiently large with respect to the constants and , while is positive under conditions (44a) and (44b). Moreover under condition (44c), while 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 point is locally asymptotically stable in . Then, it is globally asymptotically stable, provided that the following conditions are outstanding:

Proof. Recognize the following function:where , , are positive constants to be determined later on. Clearly is a continuously differentiable real valued function for all with and . Moreover we have thatHere , . Further simplification gives thatSo by choosing the positive constants , , as given below, we getClearly conditions (47a)-(47b) guarantee that ; hence is negative definite, and hence the axial equilibrium point is globally asymptotically stable and the proof is complete.

Theorem 5. Assume that the disease free equilibrium point is locally asymptotically stable in . Then, it is a globally asymptotically stable, provided that the following conditions hold:where .

Proof. Consider the functionwhere , , are positive constants to be determined later on. Clearly is a continuously differentiable real valued function for all with and . Moreover we have thatFurther simplification gives Here and .
Now rearranging the terms of the last equation further givesSo by choosing the positive constants aswe get that Accordingly, using the given conditions (52a)–(52d) we obtain Clearly is negative definite, and hence the disease free equilibrium point is globally asymptotically stable under the given conditions and hence the proof is complete.

Theorem 6. Assume that the free predator equilibrium point is locally asymptotically stable in . Then it is globally asymptotically stable, provided that the following condition holds:

Proof. Consider the next functionwhere , , are positive constants to be determined later on. Clearly is a continuously differentiable real valued function for all with and . Moreover we have thatHere , , and are given above. Now straightforward computations give By choosing the positive constants as then we obtain thatAccordingly, using the given conditions (60a) and (60b) we obtain Clearly , which means it is negative semi-definite, and hence the predator free equilibrium point is globally stable (but not asymptotically stable) under the given conditions. Moreover, since system (3) has the maximum invariant set for if and only if conditions (60a)-(60b) hold and , by Lyapunov-Lasalle’s theorem, all the solutions starting in approach the singleton set , which is the positively invariant subset of the set where . Hence becomes attracting too; hence it is globally asymptotically stable and that completes the proof.

Theorem 7. Assume that the coexistence equilibrium point of system (3) is locally asymptotically stable in . Then it is globally asymptotically stable, provided that the following conditions hold:Here , , , and are given in the proof.

Proof. Consider the real valued functionHere , , are positive constants to be determined. Clearly is a continuously differentiable real valued function for all with and . Moreover we have thatHere , , and are given above, while , , , and . Further simplification for the above equation leads toConsequently, by choosing the positive constants , , asthe following is obtained:where and . Further simplification for the last equation gives Clearly, is negative definite function under the given conditions and hence the coexistence equilibrium point 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.It is observed that, for the set of data (74), system (3) has a globally asymptotically stable positive equilibrium point 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 27.

It is observed that for the set of data (74) with the trajectory of system (3) approaches asymptotically predator free equilibrium point, while it approaches disease free point when ; see Figure 2. Further the trajectory of system (3) approaches periodic dynamics in the plane for data (74) with as shown in Figure 3; however it approaches asymptotically the positive equilibrium point otherwise.

On the other hand, for the data (74) with and , the solution of system (3) approaches asymptotically as shown in Figures 4(a)-4(b), while it approaches asymptotically when 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 too with decreasing the infection rate , that depends on our hypothetical set of data in which we assumed that the conversion rate of predator from susceptible () is less than that from infected species () and once () enter the range the system approaches disease free point as shown in the typical figure given by Figure 4(d).

Moreover it is observed that varying the parameter 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 . Now for the parameters values given by (74) with and the trajectory of system (3) approaches asymptotically and , 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 to that of system (3) using data (74) with , while the system still is stable at positive equilibrium point otherwise.

For the parameters (74) with (similarly when ) the trajectory of system (3) approaches asymptotically , while it approaches for data (74) with (similarly when or ) as shown in Figure 6. Otherwise the system has a globally asymptotically stable positive point.

Finally for the parameters (74) with (similarly when or ) the trajectory of system (3) approaches asymptotically , while it approaches for data (74) with (similarly when or ) 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 .(2)The system has no periodic dynamics lying in the interior of ; 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 by 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 makes the solution approaches asymptotically the disease free equilibrium point. Finally further decreasing of 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 plane.(4)Increasing the infection rate parameter () 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 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 and , respectively. Finally decreasing the infection rate further () leads to approaches to axial equilibrium point .(5)It is observed that varying the parameter , 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 .(6)Increasing the death rate of the infected prey () so that and 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 or decreasing it below 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 , which stands for the half saturation constant of harvesting in Michael-Mentence harvesting function, as that obtained in case of increasing or decreasing . 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 too.(8)Finally increasing the parameter , 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 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 as that effect obtained with varying .

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

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.