On the Dynamics of an Eco-Epidemiological System Incorporating a Vertically Transmitted Infectious Disease

An eco-epidemiological system incorporating a vertically transmitted infectious disease is proposed and investigated. Micheal-Mentence type of harvesting is utilized to study the harvesting effort imposed on the predator. All the properties of the solution of the system are discussed. The dynamical behaviour of the system, involving local stability, global stability, and local bifurcation, is investigated. The work is finalized with the numerical simulation to observe the global behaviour of the solution.


1.
Introduction Most of the real world problems, including biological and epidemiological problems, can be formulated mathematically using the differential equations or difference equations. The application of mathematical modeling in biology provides models known as mathematical biology models or ecological models. However, the application of mathematical models in epidemiology provides models known as epidemiological models. Since the environment contains millions of species that interact with each other and may have different types of diseases, the mathematical models that combine both ecology and epidemiology are known as eco-epidemiology models. It is well known that there are two different modes of pathogen transmission, namely the horizontal and vertical transmission modes. Horizontal transmission means the transmission of disease between the individuals of the same generation, while the transmission of disease from parent to offspring is known as vertical transmission [1]. Although most of the epidemic models are interested in the horizontal transmission type of disease, there are few studies that are interested to study epidemic models with diseases transmitted vertically; see for example [2][3][4] and the references therein. Later on, ISSN: 0067-2904 Satar et al. Iraqi Journal of Science, 2021, Vol. 62, No. 5, pp: 1642-1658 3461 Naji and Hussien [5] studied the dynamics of the spread of infectious diseases within an epidemic system. They considered both horizontal and vertical transmission in the host population. Keeping the above in view, many studies have been performed, in which the researchers presented eco-epidemiological models where the diseases exert horizontal transmission [6][7][8][9]. Later on, Sieber et al. [10] considered an eco-epidemiological model incorporating differential competition. They reported that the existence of differential competition can tremendously change the stability and persistence of predator-prey systems. Kant and Kumar [11] suggested and investigated an ecoepidemiological model involving stages that are structured with linear functional response. They assumed that the stages are existing in prey and predator, while the infection occurs in the prey population only. Das [12] proposed and studied an eco-epidemiological system in which the disease exists in the predator population. He studied the effect of alternative food on the system dynamics. Saifuddin et al. [13] considered the existence of strong-Allee effect on a simple eco-epidemiological model. They studied the dynamics of the suggested system under the combined influence of strong-Allee parameter and competition coefficients. Abdulghafour and Naji [14] proposed and studied an eco-epidemiological model incorporating a prey refuge and nonlinear harvesting from the predator. They assumed that the feeding process do not transfer the disease from prey to predator. Shaikh et al. [15] proposed and studied an eco-epidemiological model involving a virus disease. In all these proposed eco-epidemiological models and many others, in addition to the consideration of horizontal transmission type of diseases, different factors are included, such as harvesting, vaccines, toxicants, etc. Recently, Abdul Star and Naji [16] suggested and studied a prey-predator system incorporating a vertically transmitted infectious disease in predator population only. Recalling the above studies, in this paper, however, an eco-epidemiological system is suggested so that it involves a disease transmitted vertically as well as horizontally within predator species. It is assumed that the predator is falling under the effects of harvesting of nonlinear type. In fact, we used the harvesting function proposed by [17]. The paper is organized as follows; section (2) includes the formulation of the model and its dimensionless, in addition to the properties of the solution. In Section (3), the stability analysis of the system is carried out. The local bifurcation analysis is investigated in section (4). Section (5) provides a numerical simulation. Finally, Section (6) gives some conclusions and discussion on the obtained results.

The formulation and dimensionless of the model
An eco-epidemiological model incorporating a vertically transmitted infectious disease and harvesting in a predator population is formulated and studied. Consider the following hypotheses, which are adopted in the formulation of the model: 1.
Let ( ), ( ) and ( ) represent the densities at time of the prey, susceptible predators, and infected predators, respectively. 2.
Let ( ) grows logistically, in the absence of predation, with a growth rate and carrying capacity .

3.
Assume that ( ) and ( ) consume ( ) according to Holling-type II functional response, with maximum attack rates half saturation level and conversion rates . They decay exponentially in the absence of prey species, according to natural death rates .

4.
It is assumed that the disease is transmitted vertically in the predator species, in addition to transmission by contact, with infection rate .

5.
Finally, ( ) and ( ) are assumed to be harvested according to Micheal-Mentence type of harvesting function, where represents hunting efforts, while are the catchability coefficients of the predator and are positive constants. Accordingly, the dynamics of the above described eco-epidemiological model can be described using the following set of differential equations: . / (1) . Now, the number of parameters can be reduced using the following dimensionless variables and parameters Therefore, system (1) can be written in the following dimensionless form: Hence, system (2) has the following domain: *( ) +. Clearly, all the right-hand side functions are continuous and have continuous partial derivatives. Hence, the solution of system (2) exists and is unique. Now, the solution of system (2) is proved to be uniformly bounded, as shown in the following theorem. Theorem 1. The solution ( ( ) ( ) ( )) of system (2), starting at any initial condition that belong to is uniformly bounded in the region where is given in the proof. Analysis of the stability In this section, the stability analysis of system (2) is carried out through computing all the possible steady-state points, and then their type of stability is discussed. Direct computation shows that system (2) has the following steady-state points. The first steady-state point, denoted by ( ) and the second steady-state point, denoted by ( ) always exist. The third steady-state point is denoted by while ̅ is a positive root of the following three degree polynomial equation: , . Clearly, exists uniquely in the first quadrant of -plane, provided that any set of the following sets of conditions holds: . (5a) .
(5b) The fourth steady-state point is denoted by while ̂ is a positive root of the following equation: , -. Clearly, exists uniquely in the first quadrant of -plane, provided that any set of the following sets of conditions holds: . (8a) .
(8b) Finally, the fifth steady-state point is denoted by ( ) and located at the interior of positive octant , where ( )( ) , (9a) while ( ) represents the intersection point of the following two isoclines in the interior of first quadrant of plane.
Obviously, as , the above two isoclines become , . Direct computation shows that the above isoclines, which are represented by (9d) and (9e), intersect the axis at the positive points and respectively, provided that the following sets of sufficient conditions are satisfied: , Accordingly, the two isoclines (9b) and (9c) intersect each other at ( ), which belong to the interior of the first quadrant of plane, provided that the following sufficient conditions hold: Consequently, the fifth steady-state point exists in the interior of uniquely, provided that, in addition to conditions (10a)-(10e), the following condition should be hold: (10f) Now, the local stability around the above steady-state points is studied using the linearization technique. The variational matrix (VM) of system (2) around the first steady-state point ( ) has the following eigenvalues: .
Accordingly, the first steady-state point is a saddle point.
The VM evaluated at the second steady-state point ( ) is written as Therefore, the eigenvalues of ( ) are given by , Therefore, the second steady-state point ( ) is locally asymptotically stable (LAS) if the following necessary and sufficient conditions hold: . /, . /.
(14b) The VM evaluated at the third steady-state point Hence, one of the eigenvalues of ( ) is where Therefore, from the above eigenvalues and , the third steady-state point is LAS if the following conditions hold: Furthermore, the VM evaluated at the fourth steady-state point ( ̂ ) is written as Clearly, ( ) has three eigenvalues, one of them is written as , while the other two eigenvalues are the roots of the equation , where ̂ ) . Again, this equation has the roots √ and √ .
Hence, the fourth steady-state point is LAS if the following conditions hold: Finally, the stability analysis around the fifth steady-state point ( ) is studied in the next theorem. Theorem 2: The fifth steady-state point where all the symbols are clearly described in the proof. Proof: The VM around the fifth steady-state point where , , , Moreover, it is easy to verify that . Direct computation yields that the sufficient conditions (21a)-(21d) with the sign of VM elements guarantee that , , , and . Therefore, using the conditions (21e)-(21f) guarantees the positivity of , , and . Hence, due to the Routh-Hurwitz criterion, the fifth steady-state point ( ) is LAS.
In the following, Lyapunov method is used to determine the basin of attraction for each steady-state point. Then we will say that the steady-state is globally asymptotically stable (GAS) if its basin of attraction can be extended to the whole domain of system (2).

Theorem 3:
Let the second steady-state point is LAS. Then it is GAS if, in addition to condition (21d), the following condition holds:  (21d) and (24) guarantee that is negative definite. Also, since approaches to infinity if and only if any one of their variables or then it is radially unbounded Laypunov function. Therefore, the second steady-state point is GAS. Theorem 4: Assume that the third steady-state point is LAS. Then it is GAS if the following sufficient conditions hold: Note that the conditions (25a)-(25d) guarantee that is negative definite. Hence, is radially unbounded Laypunov function. Therefore, the third steady-state point is GAS. Theorem 5. Assume that the fourth steady-state point is LAS. Then it is GAS if the following sufficient conditions hold: It is easy to verify that conditions (26a)-(26d) guarantee that is negative definite. Then is radially unbounded Laypunov function. Accordingly, the fourth steady-state point is GAS. Theorem 6. Assume that the fifth steady-state point is LAS. Then it is GAS if the following sufficient conditions hold: , and . Note that conditions (27a)-(27c) ensure that is negative definite. Hence, is radially unbounded Laypunov function. Consequently, the fifth steady-state point is GAS.

4.
The occurrence of bifurcation Inthissection,theSotomoyor'stheorem [18]forlocalbifurcationisperformedtostudythepossibility of the occurrence of local bifurcation near the steady-state points of system (2). Now, for simplifying the notations of system (2), we rewrite it in the vector form as follows ( ), with ( ) and ( ) .
So, according to VM of system (2) at the point , it is easy to verify that, for any vector ( ) , we have that where .
On the other hand, we have that where ) . Then, in the following theorems, the occurrence of local bifurcation with specifying their type at each steady-state point is discussed.
where and are given in the proof. Otherwise, system (2) has a pitchfork bifurcation. Proof: By using similar arguments used in the proof of theorem (7), we obtain the following results. The VM of system (2) at with is ( ) The eigenvalues of are where and are given in the proof. Otherwise, it has a pitchfork bifurcation. Proof: It is easy to verify that the VM of system (2) at with is The eigenvalues of are , -, -, which are negative due to condition (20b) and .
(32b) Then, when the parameter passes through the value , system (2) around the fifth steady-state point  undergoes a saddle-node bifurcation if   , -, -, - Proof: From Equation (22), it is easy to verify that the VM of system (2) at with is where ,for all with , -( ). Also, it is clear that, at , we obtain that , where is given in Equation (23). Therefore, the VM of system (2) that is given by Equation (34)  Clearly, , ( )( )under condition (32c). Thus, the proof follows.

Simulation of the system
In this section, system (2) is simulated numerically using the following set of hypothetical data. The objective is to specify the types of attractors in system (2) and detect the control parameters.
For the data given by equation (35), it is observed that the solution of system (2) approaches asymptotically (APAS) the periodic attractor in the interior of as shown in Figure (1).  (2) coexists at a periodic attractor using equation (35). Now, for the same data with decreasing the parameter in the range , system (2) loses its persistence and the trajectory APAS to the periodic attractor in the interior of the first quadrant of plane. While upon increasing it in the range , system (2) loses its persistence, too, and the trajectory APAS the periodic attractor in the interior of the first quadrant of plane. However, for the range , the trajectory of system (2) APAS the fourth steady-state point in the interior of the first quadrant of plane. Otherwise, the system (2) still has periodic dynamics in the interior of see Figure (2) for the selected values of the parameter.   Also, it is noted that, for the range with the other of parameters being as in equation (35), system (2) loses the persistence and the solution APAS the periodic dynamics in the interior of the first quadrant of the plane. However, for system (2) faces extinction in species and it APAS the periodic dynamics in the interior of the first quadrant of the plane. While, otherwise, the solution of system (2) still approaches to the periodic attractor in the interior of , see Figure   Furthermore, the effect of harvesting on the species is shown in Figure (6), so that for with other parameters being as in equation (35), the solution of system (2) still persists at a periodic dynamics in the interior of . However, increasing the parameter further leads to extinction in the species and the solution approaches the periodic dynamics in the interior of the first quadrant of the plane. On the other hand, varying the parameters and , with the other parameters being as in equation (35), has quantitative effects on the level of population density of and respectively, in the 3D periodic dynamics of system (2). Also, it is noted that varying the parameters and has similar effects on the solution of system (2) as those shown in the case of varying and respectively. Now, for with the other parameters being as in equation (35), system (2) faces an extinction in the species and the solution of system (2) APAS the periodic dynamics in the interior of the first quadrant of the plane. Otherwise, system (2)

Discussion and conclusions
In this paper, an eco-epidemiological model is suggested and studied. The model is combining a prey-predator system with an infectious disease in the predator. It is assumed that the disease is transmitted vertically, in addition to the natural horizontal transmission within predator individuals. Moreover, there is a harvesting event on the predator using Micheal-Mentence type of harvesting function. The stability analysis of all possible steady-states is investigated using the linearization technique and Lyapunov functions. The possibility of the occurrence of local bifurcation around the steady-states of the system is investigated. Finally, the paper is ended with the numerical simulation of the model using a hypothetical set of data, as given in equation (35). Regarding the numerical simulation results that depend on the data given by equation (35), different sets of data can be used too. The following conclusions are obtained.

1.
System (2) APAS the periodic dynamics in the interior of . In fact, due to the complexity of the stability conditions given by (21a)-(21f), it is difficult to find data that satisfy all these conditions and then obtain an asymptotically stable coexistence equilibrium point, but it still exists analytically.

2.
Decreasing the half saturation constant ( ) below a specific value or increasing it above another specific value causes an extinction in one compartment of the predator species (susceptible or infected) and, then, system (2) loses its persistence and APAS the periodic dynamics or steady-state point in the boundary planes. However, it still persists at a periodic dynamics otherwise.

3.
Similar behavior as that observed in the half saturation constant is also obtained regarding the conversion rates ( and ) and the death rates of predator compartments ( and ).

4.
Increasingthecontact'sinfectionrate( ) above a specific value leads to an extinction in the susceptible predator and the solution APAS the periodic dynamics in the boundary plane. Otherwise, the system (2) persists at periodic dynamics in the interior of .

5.
Similar behavior as that observed in the contact'sinfectionrateisalsoobtainedregardingeach of the maximum harvesting rates ( and ).

6.
The parameters which stand for the prevention of the harvesting process, represented by ( and ) have quantitative effects on the levels of predator curves in the periodic dynamics that fall in the interior of . 7.
Decreasingthecontact'stransmission rate ( ) below a specific value leads to an extinction in the infected predator and the solution APAS the periodic dynamics in the boundary plane. Otherwise, the system (2) still persists at periodic dynamics in the interior of .