Dynamical of Ratio-Dependent Eco-epidemical Model with Prey Refuge CAUCHY –Jurnal Matematika Murni dan Aplikasi Volume 6(4) (2021), Pages 227-237 p-ISSN: 2086-0382; e-ISSN: 2477-3344 Submitted: November 17, 2020 Reviewed: February 11, 2021 Accepted: March 16, 2021 DOI: http://dx.doi.org/10.18860/ca.v6i4.10827 Dynamical of Ratio-Dependent Eco-epidemical Model with Prey Refuge Adin Lazuardy Firdiansyah Department of Tadris Matematika, STAI Muhammadiyah Probolinggo Jl. Soekarno Hatta 94B, Sukabumi, Mayangan, Probolinggo, Indonesia Email: adin.lazuardy@gmail.com ABSTRACT This paper discusses the dynamic analysis of three species in the eco-epidemiology model by considering the ratio-dependent function and prey refuge. The prey refuge is applied under the fact that infected prey has protection instincts that allow it to reduce predation risk. Here, we get the boundedness and three equilibrium points where are existence under certain conditions. In the model, three equilibrium points are locally asymptotically stable and one of the equilibrium points is globally asymptotically stable. We find that the system undergoes Hopf bifurcation around the interior equilibrium point by choosing prey refuge as a bifurcation parameter. We also find a condition for uniform persistence. Finally, several simulations of numerical are performed not only to illustrate the analytical results but also to illustrate the effect of the prey refuge. Keywords: eco-epidemiology model; global stable; Hopf bifurcation; local stable; persistence INTRODUCTION One of the natural phenomena that described the interaction between one species and another individual is the prey-predator interaction. This interaction depends on whether the effects are profitable or detrimental. In the real world, prey-predator interaction also can be influenced by infectious diseases. These diseases can affect population size in the predator-prey interaction. Since then, the combination of epidemiological and ecological becomes important issues that are often discussed by many researchers. Mathematical studies have considered this issue into an eco- epidemiology model that contains the class of susceptible and infective populations. Currently, several studies have focused on the spread of disease in prey only, e.g. [1]–[3]. It is well known that predators prefer to capture infected prey because they are easier to catch than susceptible prey. However, the predator can become infected after eating them. Therefore, several researchers are interested to investigate the spread of disease not only in prey but also in predators, e.g. [4][5]. Moreover, some studies have reviewed the spread of disease in both populations, e.g. [6]. Base on several experiments, the spread of infectious disease becomes an important factor to know the regulation of population density [7]. In this paper, we focus on the situation where predators can eat infected prey only. It appropriates to the fact that infected prey tends to change its behavior. Infected prey shall live in an area that is accessible to predators [8]. Moreover, infected prey is less agile than healthy prey and can be predated by predators easily [1]. Hudson et al. [9] have http://dx.doi.org/10.18860/ca.v6i4.10827 mailto:adin.lazuardy@gmail.com Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 228 observed that predators selectively capture heavily infected red grouse. Predators that consumed infected prey populations are described by the response functions. It is well known that response functions are an important component in the eco-epidemiology model. Generally, several researchers use Holling types as a response function in their model. According to [10], Holling types are divide into three types, namely Holling type I, Holling type II, and Holling type III. Holling type I means that the consumption rate of predator increases linearly with the density of prey but it achieves a constant value if predators are surfeit, e.g. [6][5]. Holling type II means that the consumption rate of predators increases if the density of prey is low, e.g. [1]–[4]. Meanwhile, Holling type III means that the consumption rate of predator increases when the density of prey is large but it decreases when the density of prey is low, e.g. [11]. In Holling type III, predators easily switch to eat one prey to another or they focus to eat prey in a location where it is most abundant [12]. The response function for the Holling type depends on the density of prey. This is unrealistic because it ignores the effects of predator interference. Base on the experiment, the density of predators can influence the consumption rate. In the modeling, the consumption rate that depended on the density of prey can’t describe the dynamics behavior when the density of predator influences the system [12]. Currently, several researchers have considered the density of both populations as a ratio-dependent function. This function depends on the ratio of prey to predator density [13]. According to [14], the ratio-dependent model is a more reasonable dynamic than the previous model. One of the phenomena that reduce predation risk is prey refuge. It can avoid the extinction of prey and influence the stability of the dynamic behavior [15]. According to [3], the prey refuge that is incorporated in the model is divided into two types, namely the refuge for a constant-number of prey and the refuge for a constant-proportion of prey. It is well known that the refuge for a constant-number of prey has a stronger stabilizing effect than the refuge for a constant-proportion of prey [14]. Therefore, the model is more realistic by incorporating prey refuge and it gives an accessible factor to the predator. In this study, we modify the eco-epidemiology model from [1] by changing Holling type II into a ratio-dependent function. Here, we also observe the effect of prey refuge in the system. Further, this article presents the results in the form of model analysis. Moreover, it is well observed that there are Hopf bifurcations around the positive equilibrium point and the condition for uniform persistence. Finally, several simulations are performed to illustrate the analytical results. METHODS We use several methods to modify the eco-epidemiology model from [1]. The method is presented as follows. 1. Reviewing and studying the eco-epidemiology model from previous literature. 2. Modifying the eco-epidemiology model by changing the Holling II type into the ratio- dependent function. 3. Investigating the boundedness, equilibrium points, and dynamical behavior in the modified model. 4. Investigating Hopf bifurcation and persistence in the modified model. 5. Performing numerical simulation by using the 5th-order predictor-corrector method as a numerical method to support the analytical results. Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 229 RESULTS AND DISCUSSION The Mathematical Model In this paper, the eco-epidemiology model consists of three populations, namely susceptible prey, infected prey, and predator. Let 𝑋𝑆(𝑑), 𝑋𝐼(𝑑), and π‘Œ(𝑑) is respectively defined as the density of susceptible prey, infected prey, and predator at the time 𝑑. According to [1], the general eco-epidemiology model is presented as follows. 𝑋�̇� = 𝑅 βˆ’ 𝛽𝑋𝑆𝑋𝐼 βˆ’ 𝛿𝑋𝑆, 𝑋𝐼̇ = 𝛽𝑋𝑆𝑋𝐼 βˆ’ 𝑓(𝑋𝐼, π‘Œ)π‘Œ βˆ’ πœ‚π‘‹πΌ, οΏ½Μ‡οΏ½ = 𝑒𝑓(𝑋𝐼, π‘Œ)π‘Œ βˆ’ π›Ύπ‘Œ, (1) with 𝑋𝑆 β‰₯ 0, 𝑋𝐼 β‰₯ 0, π‘Œ β‰₯ 0. The first equation of system (1) expresses that in the absence of disease, the prey population grows by following the equation as below. 𝑋�̇� = 𝑅 βˆ’ 𝛿𝑋𝑆, where 𝑅 is expressed as the level of recruitment in prey populations such as immigrants and new-born and 𝛿 is defined as the natural death rate of susceptible prey. Here, the growth of the prey population is affected by factors such as disease. The spread of disease is denoted by bilinear incidence rate 𝛽𝑋𝑆𝑋𝐼 with 𝛽 is the transmission rate. We assume that the prey is not infected because of inherited disease but other sources. Moreover, the disease spreads on susceptible prey only. We also assume that the infected prey populations do not recover nor reproduce. The second equation of system (1) describes the development of the infected prey population. They will be erased by the natural death rate πœ‚ and the predation of the predator. Here, predation is denoted by the response function 𝑓(𝑋𝐼, π‘Œ). We assume that infected prey can hide. Hiding behavior gives protection for infected prey which can protect from predation. The protection of infected prey is denoted by constant π‘š. Predators can capture infected prey by following a ratio-dependent function as in [14]. 𝑓(𝑋𝐼, π‘Œ) = { 0 , if 0 ≀ 𝑋𝐼 ≀ π‘š, π‘Ž(𝑋𝐼 βˆ’ π‘š) 𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ , if 𝑋𝐼 > π‘š, with π‘Ž is expressed as the predation rate and πœ‰ is defined as half capturing saturation constant. According to [1], if the density of infected prey populations is below the constant π‘š, then predators cannot eat them and will die exponentially. Meanwhile, if the density of infected prey populations is above the constant π‘š, then predators can eat them. Thus, when 𝑋𝐼 > π‘š, then system (1) shall become as follows. 𝑋�̇� = 𝑅 βˆ’ 𝛽𝑋𝑆𝑋𝐼 βˆ’ 𝛿𝑋𝑆, 𝑋𝐼̇ = 𝛽𝑋𝑆𝑋𝐼 βˆ’ π‘Ž(𝑋𝐼 βˆ’ π‘š)π‘Œ 𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ βˆ’ πœ‚π‘‹πΌ, οΏ½Μ‡οΏ½ = π‘Žπ‘’(𝑋𝐼 βˆ’ π‘š)π‘Œ 𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ βˆ’ π›Ύπ‘Œ, (2) with 𝑋𝑆(0) β‰₯ 0, 𝑋𝐼(0) β‰₯ 0, π‘Œ(0) β‰₯ 0. The last equation represents the behavior of the predator population. Predators only consume the infected prey and don’t consume the susceptible prey. When the infected prey population is absent, then the predator experiences a natural death rate 𝛾. Here, it is assumed that disease does not spread from infected prey to predator. In the model, all parameters are positive values. The meaning of parameter and their units are summarized in Table 1. Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 230 Table 1. Units and the meaning of parameters in system (2) Parameters Biological meaning Units 𝑋𝑆 The number of susceptible prey Number 𝑋𝐼 The number of infected prey Number π‘Œ The number of predators Number 𝑅 The level of recruitment in prey time-1 𝛽 The level of infection of disease mass-1 time-1 𝛿 The level of natural mortality in susceptible prey time-1 πœ‚ The level of natural mortality in infected prey time-1 𝛾 The level of natural mortality in predators time-1 πœ‰ Half saturation constant Number π‘Ž The level of predation in predators mass-1 time-1 𝑒 The level of alteration from prey into predators time-1 π‘š The measure of prey in the refuge Number The Boundedness To show the biological validity, we shall prove the boundedness of the system (2) as follows. Theorem 1. All solutions of the system (2) are uniformly bounded. Proof: We define π‘Š = 𝑋𝑆 + 𝑋𝐼 + 1 𝑒 π‘Œ. By differentiating π‘Š to 𝑑, we get π‘‘π‘Š 𝑑𝑑 = 𝑑𝑋𝑆 𝑑𝑑 + 𝑑𝑋𝐼 𝑑𝑑 + 1 𝑒 π‘‘π‘Œ 𝑑𝑑 . By substituting system (3), we get, π‘‘π‘Š 𝑑𝑑 = 𝑅 βˆ’ (𝛿𝑋𝑆 + πœ‚π‘‹πΌ + 𝛾 π‘Œ 𝑒 ). Next, we choose π‘ž = min{𝛿, πœ‚, 𝛾}. Thus, we get π‘‘π‘Š 𝑑𝑑 ≀ 𝑅 βˆ’ π‘žπ‘Š, By using the theory of differential equation, we get π‘Š(𝑑) ≀ 𝑅 π‘ž + πΆπ‘’βˆ’π‘žπ‘‘, where 𝐢 is the arbitrary positive constant. For 𝑑 β†’ ∞, we get lim π‘‘β†’βˆž sup π‘Š(𝑑) ≀ 𝑅 π‘ž . Thus, all solutions of system (2) enter into 𝛺 = {(𝑋𝑆, 𝑋𝐼, π‘Œ) ∈ ℝ+ 3 : π‘Š(𝑑) ≀ 𝑅 π‘ž }. The Equilibrium Points By setting 𝑋�̇� = 𝑋𝐼̇ = οΏ½Μ‡οΏ½ = 0, we get the possible equilibrium points as below. 1. Axial equilibrium 𝐸0 ( 𝑅 𝛿 , 0,0) always existent. Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 231 2. Planar equilibrium 𝐸1 ( πœ‚ 𝛽 , π›½π‘…βˆ’π›Ώπœ‚ πœ‚π›½ , 0) exists when 𝛽𝑅 > π›Ώπœ‚. 3. Interior equilibrium πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—), where 𝑋𝑆 βˆ— = 𝑅 𝛽𝑋𝐼 βˆ—+𝛿 , π‘Œβˆ— = (π‘Žπ‘’βˆ’π›Ύ)(𝑋𝐼 βˆ—βˆ’π‘š) π›Ύπœ‰ , and 𝑋𝐼 βˆ— is the positive root of the quadratic equation (3) as follows. 𝐴1(𝑋𝐼 βˆ—)2 + 𝐴2𝑋𝐼 βˆ— + 𝐴3 = 0, (3) with 𝐴1 = βˆ’π›½(π‘Žπ‘’ βˆ’ 𝛾 + π‘’πœ‚πœ‰), 𝐴2 = π›½π‘š(π‘Žπ‘’ βˆ’ 𝛾) βˆ’ 𝛿(π‘Žπ‘’ βˆ’ 𝛾 + π‘’πœ‚πœ‰) + π‘’π‘…π›½πœ‰, 𝐴3 = π›Ώπ‘š(π‘Žπ‘’ βˆ’ 𝛾). This point πΈβˆ— exists when it satisfies the condition π‘Žπ‘’ > 𝛾 and 𝑋𝐼 βˆ— > π‘š. It is clear to show that 𝐴1 < 0 and 𝐴3 > 0. Thus, the determinant of the equation (3) is 𝐷 = (𝐴2) 2 βˆ’ 𝐴1𝐴3 β‰₯ 0. Therefore, to get the explicit form of the root of the equation (3), we have to check the following two cases: a. For 𝐷 = 0, the equation (3) has a twin positive root where 𝑋𝐼 βˆ— = βˆ’ 𝐴2 2𝐴1 with 𝐴2 > 0. b. For 𝐷 > 0, the probability that equation (3) has positive roots is as follows. ο‚· If 𝐴2 > 0, then the equation (3) has two positive roots where 𝑋𝐼1,2 βˆ— = βˆ’π΄2±√𝐷 2𝐴1 . ο‚· If 𝐴2 < 0, then the equation (3) has a single positive root where 𝑋𝐼 βˆ— = βˆ’π΄2βˆ’βˆšπ· 2𝐴1 . Dynamical Behaviour To investigate the stability in system (2), we have to determine the eigenvalues of the Jacobian matrix. Here, we identify the Jacobian matrix at 𝐸(𝑋𝑆, 𝑋𝐼, π‘Œ) as follows. 𝐽(𝐸) = [ βˆ’π›½π‘‹πΌ βˆ’ 𝛿 βˆ’π›½π‘‹π‘† 0 𝛽𝑋𝐼 𝛽𝑋𝑆 βˆ’ π‘Žπœ‰π‘Œ2 (𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ) 2 βˆ’ πœ‚ βˆ’ π‘Ž(𝑋𝐼 βˆ’ π‘š) 2 (𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ) 2 0 π‘Žπ‘’πœ‰π‘Œ2 (𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ) 2 π‘Žπ‘’(𝑋𝐼 βˆ’ π‘š) 2 (𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ) 2 βˆ’ 𝛾 ] . (4) To check the stability of 𝐸0 ( 𝑅 𝛿 , 0,0), we get the Jacobian matrix by replacing 𝐸(𝑋𝑆, 𝑋𝐼, π‘Œ) in the equation (4) with 𝐸0 ( 𝑅 𝛿 , 0,0). Hence, we obtain the eigenvalues of the Jacobian matrix, namely πœ†1 = βˆ’π›Ώ, πœ†2 = 𝛽𝑅 𝛿 βˆ’ πœ‚, and πœ†3 = π‘Žπ‘’ βˆ’ 𝛾. The point 𝐸0 is locally asymptotically stable when π‘Žπ‘’ < 𝛾 and 𝛽𝑅 < π›Ώπœ‚. To investigate the stability of 𝐸1 ( πœ‚ 𝛽 , π›½π‘…βˆ’π›Ώπœ‚ πœ‚π›½ , 0), we have to identify the Jacobian matrix by replacing 𝐸(𝑋𝑆, 𝑋𝐼, π‘Œ) in the equation (4) with 𝐸1 ( πœ‚ 𝛽 , π›½π‘…βˆ’π›Ώπœ‚ πœ‚π›½ , 0). Thus, we get the eigenvalues πœ†1 = π‘Žπ‘’ βˆ’ 𝛾 and other eigenvalues are the roots of the quadratic equation πœ†2 + πœ‘1πœ† + πœ‘2 = 0 with πœ‘1 = π›½π‘…βˆ’π›Ώπœ‚ πœ‚ + 𝛿 and πœ‘2 = 𝛽𝑅 βˆ’ π›Ώπœ‚. By using the Routh-Hurwitz criteria, the eigenvalues have negative real roots when 𝛽𝑅 > π›Ώπœ‚. Hence, the point 𝐸1 is locally asymptotically stable when π‘Žπ‘’ < 𝛾 and 𝛽𝑅 > π›Ώπœ‚. From the above discussion, we get the following theorem as follows. Theorem 2. The axial equilibrium 𝐸0 is locally asymptotically stable when π‘Žπ‘’ < 𝛾 and 𝛽𝑅 < π›Ώπœ‚ and the planar equilibrium 𝐸1 is locally asymptotically stable when π‘Žπ‘’ < 𝛾 and 𝛽𝑅 > π›Ώπœ‚. Theorem 2 means that if the level of natural mortality in predator is lower than a certain value and the level of infection is lower than a certain value, then the point 𝐸0 Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 232 becomes stable. Meanwhile, if the level of infection is greater than a certain value, then the point 𝐸1 becomes stable. Now, we shall investigate the stability of πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—). By replacing 𝐸(𝑋𝑆, 𝑋𝐼, π‘Œ) in the equation (4) with πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—), we get the Jacobian matrix where 𝑒𝑖𝑗 is the entry of matrix with row 𝑖 and column 𝑗 as follows. 𝐽(πΈβˆ—) = [ 𝑒11 𝑒12 0 𝑒21 𝑒22 𝑒23 0 𝑒32 𝑒33 ], with 𝑒11 = βˆ’π›½π‘‹πΌ βˆ— βˆ’ 𝛿; 𝑒12 = βˆ’π›½π‘‹π‘† βˆ—; 𝑒21 = 𝛽𝑋𝐼 βˆ—, 𝑒22 = 𝛽𝑋𝑆 βˆ— βˆ’ π‘Žπœ‰(π‘Œβˆ—)2 (𝑋𝐼 βˆ— βˆ’ π‘š + πœ‰π‘Œβˆ—)2 βˆ’ πœ‚; 𝑒23 = βˆ’ π‘Ž(𝑋𝐼 βˆ— βˆ’ π‘š)2 (𝑋𝐼 βˆ— βˆ’ π‘š + πœ‰π‘Œβˆ—)2 , 𝑒32 = π‘Žπ‘’πœ‰(π‘Œβˆ—)2 (𝑋𝐼 βˆ— βˆ’ π‘š + πœ‰π‘Œβˆ—)2 ; 𝑒33 = βˆ’ π‘Žπ‘’πœ‰π‘Œβˆ—(𝑋𝐼 βˆ— βˆ’ π‘š) (𝑋𝐼 βˆ— βˆ’ π‘š + πœ‰π‘Œβˆ—)2 . Hence, we obtain the characteristic equation of πΈβˆ—, namely πœ†3 + πœ‡1πœ† 2 + πœ‡2πœ† + πœ‡3 = 0 (5) with πœ‡1 = βˆ’(𝑒11 + 𝑒22 + 𝑒33), πœ‡2 = 𝑒11𝑒22 + 𝑒11𝑒33 + 𝑒22𝑒33 βˆ’ 𝑒12𝑒21 βˆ’ 𝑒23𝑒32, πœ‡3 = βˆ’π‘’11𝑒22𝑒33 + 𝑒11𝑒23𝑒32 + 𝑒12𝑒21𝑒33. By using the Routh-Hurwitz criteria, the eigenvalues have negative real roots when πœ‡1 > 0, πœ‡3 > 0, πœ‡1πœ‡2 > πœ‡3. Thus, the point 𝐸 βˆ— is locally asymptotically stable when πœ‡1 > 0, πœ‡3 > 0, πœ‡1πœ‡2 > πœ‡3. Therefore, we get the following theorem. Theorem 3. The interior equilibrium πΈβˆ— is locally asymptotically stable when it satisfies the condition πœ‡1 > 0, πœ‡3 > 0, πœ‡1πœ‡2 > πœ‡3. In the next theorem, we shall prove that the point 𝐸1 is globally asymptotically stable under a certain condition. Theorem 4. The point 𝐸1 is globally asymptotically stable in the 𝑋𝑆 βˆ’ 𝑋𝐼 plane. Proof: By applying the Dulac function as 𝐻(𝑋𝑆, 𝑋𝐼) = 1 𝑋𝐼 , we have 𝐻(𝑋𝑆, 𝑋𝐼) = 1 𝑋𝐼 , β„Ž1(𝑋𝑆, 𝑋𝐼) = 𝑅 βˆ’ 𝛽𝑋𝑆𝑋𝐼 βˆ’ 𝛿𝑋𝑆, β„Ž2(𝑋𝑆, 𝑋𝐼) = 𝛽𝑋𝑆𝑋𝐼 βˆ’ πœ‚π‘‹πΌ, where 𝐻(𝑋𝑆, 𝑋𝐼) > 0 in the 𝑋𝑆 βˆ’ 𝑋𝐼 plane. Thus, we get βˆ†(𝑋𝑆, 𝑋𝐼) = πœ• πœ•π‘‹π‘† (β„Ž1𝐻) + πœ• πœ•π‘‹πΌ (β„Ž2𝐻) = βˆ’π›½ βˆ’ 𝛿 𝑋𝐼 < 0. Base on Bendixson-Dulac criteria, there is no limit cycle in the 𝑋𝑆 βˆ’ 𝑋𝐼 plane. Thus, the point 𝐸1 is globally asymptotically stable in the 𝑋𝑆 βˆ’ 𝑋𝐼 plane. Hopf Bifurcation In this section, we shall investigate Hopf bifurcation around πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—). Hopf bifurcation guarantees that all solutions of system (2) enter a limit cycle around πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—). Here, we choose the constant π‘š as a bifurcation parameter. Hopf bifurcation around πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—) is presented in the following theorem. Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 233 Theorem 5. The system (2) undergoes a Hopf bifurcation around πΈβˆ—(𝑋𝑆 βˆ—, 𝑋𝐼 βˆ—, π‘Œβˆ—) when π‘š passes through a critical value π‘š = π‘šπ‘. Proof: We consider equation (5) as the characteristic equation at πΈβˆ—. Next, we choose π‘š = π‘šπ‘ such that πœ‡1πœ‡2 = πœ‡3 where πœ‡1, πœ‡2, πœ‡3 > 0. Therefore, we get (πœ†2 + πœ‡2)(πœ† + πœ‡1) = 0 (6) with the roots are πœ†1,2 = Β±π‘–βˆšπœ‡2 and πœ†3 = βˆ’πœ‡1. For all π‘š, the roots become πœ†1,2 = 𝑣1(π‘š) Β± 𝑖𝑣2(π‘š) and πœ†3 = βˆ’πœ‡1(π‘š) where 𝑣1(π‘š) and 𝑣2(π‘š) are real. Next, we shall prove the transversality condition as follows. 𝑑 (𝑅𝑒 πœ†π‘—(π‘š)) π‘‘π‘š | π‘š=π‘šπ‘ β‰  0, 𝑗 = 1,2. By substituting πœ†1 = 𝑣1(π‘š) + 𝑖𝑣2(π‘š) into equation (6) and differentiating to π‘š, we obtain 𝐾𝑣1Μ‡ βˆ’ 𝐿𝑣2Μ‡ + 𝑀 = 0, 𝐿𝑣1Μ‡ + 𝐾𝑣2Μ‡ + 𝑁 = 0, (7) where 𝐾 = 3(𝑣1 2 βˆ’ 𝑣2 2) + πœ‡2 + 2𝑣1πœ‡1, 𝐿 = 6𝑣1𝑣2 + 2𝑣2πœ‡1, 𝑀 = πœ‡1Μ‡(𝑣1 2 βˆ’ 𝑣2 2) + 𝑣1πœ‡2Μ‡ + πœ‡3Μ‡, and 𝑁 = 2𝑣1𝑣2πœ‡1Μ‡ + 𝑣2πœ‡2Μ‡. Next, we solve equation (7). Thus, we have 𝑣1Μ‡ = βˆ’ 𝐾𝑀 + 𝐿𝑁 𝐾2 + 𝐿2 . Since 𝐾𝑀 + 𝐿𝑁 β‰  0 and πœ‡1Μ‡ β‰  0, we obtain 𝑑 (𝑅𝑒 πœ†π‘—(π‘š)) π‘‘π‘š | π‘š=π‘šπ‘ = βˆ’ 𝐾𝑀 + 𝐿𝑁 𝐾2 + 𝐿2 | π‘š=π‘šπ‘ β‰  0, 𝑗 = 1,2, and πœ†3(π‘šπ‘) = βˆ’πœ‡1(π‘šπ‘) < 0. Thus, the system (2) occurs Hopf bifurcation when π‘š passes through a critical value π‘š = π‘šπ‘. Persistence To show that all species are present and are not extinct in the future time, we shall prove that system (2) is uniform persistence. Theorem 6. Let the assumption of Theorem 4 holds. If the inequalities π‘Žπ‘’ > 𝛾 and 𝛽𝑅 > π›Ώπœ‚ hold, then the system (2) is uniform persistence. Proof: We consider average Lyapunov function 𝜎(𝑋) = 𝑋𝑆 π‘Ÿ1𝑋𝐼 π‘Ÿ2π‘Œπ‘Ÿ3 with π‘Ÿ1, π‘Ÿ2, π‘Ÿ3 > 0. Here, 𝜎(𝑋) is nonnegative 𝐢1 in ℝ+ 3 . Thus, we have πœ—(𝑋) = 1 𝜎 π‘‘πœŽ 𝑑𝑑 = π‘Ÿ1 𝑋�̇� 𝑋𝑆 + π‘Ÿ2 𝑋𝐼̇ 𝑋𝐼 + π‘Ÿ3 οΏ½Μ‡οΏ½ π‘Œ , = π‘Ÿ1 ( 𝑅 𝑋𝑆 βˆ’ 𝛽𝑋𝐼 βˆ’ 𝛿) + π‘Ÿ2 (𝛽𝑋𝑆 βˆ’ π‘Ž(𝑋𝐼 βˆ’ π‘š)π‘Œ 𝑋𝐼(𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ) βˆ’ πœ‚) + π‘Ÿ3 ( π‘Žπ‘’(𝑋𝐼 βˆ’ π‘š) 𝑋𝐼 βˆ’ π‘š + πœ‰π‘Œ βˆ’ 𝛾). In the system, the point 𝐸1 is the only equilibrium point which is no limit cycle around the equilibrium point. Therefore, Theorem 4 holds. Hence, it is enough to prove that πœ—(𝑋) > 0 for all equilibrium point 𝑋 ∈ 𝑏𝑑 ℝ+ 3 . Thus, we get πœ—(𝐸0) = π‘Ÿ2 ( 𝛽𝑅 βˆ’ π›Ώπœ‚ 𝛿 ) + π‘Ÿ3(π‘Žπ‘’ βˆ’ 𝛾) > 0, πœ—(𝐸1) = π‘Ÿ3(π‘Žπ‘’ βˆ’ 𝛾) > 0. Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 234 We note that if the inequalities π‘Žπ‘’ > 𝛾 and 𝛽𝑅 > π›Ώπœ‚ hold, then πœ—(𝐸0) > 0 holds. Meanwhile, if the inequalities π‘Žπ‘’ > 𝛾 holds, then πœ—(𝐸1) > 0 holds. Therefore, Theorem 6 expresses the instability of the point 𝐸0 and 𝐸1. Numerical Solutions The analytical results obtained are incomplete without numerical investigation. In this section, we present the numerical solution by using the 5th-order predictor-corrector method at βˆ†π‘‘ = 0.01. Here, we will give four simulations to verify our analytical results and also to demonstrate the effect of the prey refuge. We choose several parameters as in (8). Their units are given as in Table 1. 𝑅 = 2, 𝛽 = 1, 𝛿 = 1, πœ‚ = 0.5, 𝛾 = 0.5, πœ‰ = 1, π‘Ž = 2, 𝑒 = 0.75, π‘š = 0.5. (8) Simulation 1. It confirms that system (2) has an equilibrium πΈβˆ—(1.0685,0.8717,0.7434). On simulation 1, Theorem 3 and Theorem 6 are satisfied. We observe that all solutions converge to the point πΈβˆ—, see figure 1(a). Thus, the point is locally asymptotically stable, see figure 1(b), and the system (2) is uniform persistence. (a) (b) Figure 1. The Dynamics of System (2) with π‘š = 0.5 and Other Parameters as in (8) Simulation 2. We replace π‘š = 0.5 into π‘š = 0.0002. Here, we investigate that system (2) has an equilibrium point πΈβˆ—(1.8305,0.0926,0.1848). On simulation 2, Theorem 3 is not satisfied but Theorem 6 is satisfied. Therefore, system (2) is uniform persistence but the equilibrium point πΈβˆ— is unstable, see figure 2. (a) (b) Figure 2. This Figure Shows The Instability of System (2) Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 235 Simulation 3. To see Hopf bifurcation in the system (2), we choose π‘š = π‘šπ‘ = 0.0004 as a bifurcation parameter. Thus, we identify that the equilibrium point of system (2) is πΈβˆ—(1.8277,0.0943,0.1878). If we choose π‘š = 0.005 > π‘šπ‘ = 0.0004, then all solutions of the system (2) convergent to πΈβˆ—(1.7795,0.1239,0.2378), see figure 3(a). Meanwhile, figure 3(b) shows the phase portrait of the system (2) with π‘š = 0.005 which means that the point πΈβˆ— is stable. Furthermore, if we choose π‘š = 0.0001 < π‘šπ‘ = 0.0004, then the point πΈβˆ—(1.8319,0.0918,0.1833) is unstable, see figure 4(a). Meanwhile, the phase portrait in the system (2) with π‘š = 0.0001 that is presented in figure 4(b) means that the solution of the system (2) enter a limit cycle around πΈβˆ—. Thus, Theorem 5 is satisfied. Simulation 4. To see the effect of prey refuge in the system (2), we use several values of prey refuge, namely π‘š1 = 0.1, π‘š2 = 0.65, and π‘š3 = 1.3. Figure (5) shows the time graph of the system (2) by using π‘š makes different values. Here, we obtain that all populations exist no matter how large π‘š with π‘š < 𝑋𝐼. This prey refuge creates the system (2) to become stable rapidly and no extinction occurs. Here, it is worthy to attention that when we choose the constant π‘š with π‘š < 𝑋𝐼, then the measure of prey refuge doesn’t lead to predator extinction. (a) (b) Figure 3. The Dynamics of System (2) with π‘š = 0.005 and Other Parameters as in (8) (a) (b) Figure 4. The Dynamics of System (2) with π‘š = 0.0001 and Other Parameters as in (8) Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 236 Figure 5. The Effect of Prey Refuge by Using π‘š Makes Different Values CONCLUSIONS In this study, we have observed three species in the eco-epidemiology model with the ratio-dependent function incorporating prey refuge. We obtain three equilibrium points where all points, i.e. axial equilibrium, planar equilibrium, and interior equilibrium, are locally asymptotically stable under certain conditions. Moreover, the planar equilibrium point in the system (2) is globally asymptotically stable. Next, we find that Hopf bifurcation occurs around the interior equilibrium by choosing a bifurcation parameter in the constant π‘š. When π‘š > π‘šπ‘ = 0.0004, the system (2) is stable. However, when π‘š < π‘šπ‘ = 0.0004, the system (2) is unstable. Furthermore, we also find a condition for uniform persistence. If the level of natural mortality in predators is lower than a certain value and the level of infection is greater than a certain value, then all species exist in the future time. Next, by applying the prey refuge, all populations exist no matter how large π‘š with π‘š < 𝑋𝐼. We conclude that system (2) is stable faster and there is no extinction. REFERENCES [1] C. Maji, D. Kesh, and D. Mukherjee, β€œBifurcation and global stability in an eco- epidemic model with refuge,” Energy, Ecol. Environ., vol. 4, no. 3, pp. 103–115, 2019, doi: 10.1007/s40974-019-00117-6. [2] A. K. Pal and G. P. Samanta, β€œStability analysis of an eco-epidemiological model incorporating a prey refuge,” Nonlinear Anal. Model. Control, vol. 2014, 2014, doi: 10.1155/2014/978758. [3] B. Mukhopadhyay and R. Bhattacharyya, β€œEffects of deterministic and random refuge in a prey-predator model with parasite infection,” Math. Biosci., vol. 239, no. 1, pp. 124–130, 2012, doi: 10.1016/j.mbs.2012.04.007. [4] C. Huang, H. Zhang, J. Cao, and H. Hu, β€œStability and Hopf Bifurcation of a Delayed Prey-Predator Model with Disease in the Predator,” Int. J. Bifurc. Chaos, vol. 29, no. 7, p. 23, 2019, doi: 10.1142/S0218127419500913. [5] S. A. Wuhaib and Y. A. B. U. Hasan, β€œPredator-Prey Interactions With Harvesting of Predator With Prey in Refuge,” Commun. Math. Biol. Neurosci., vol. 2013, no. 1, pp. 1–19, 2013. [6] S. P. Bera, A. Maiti, and G. P. Samanta, β€œA prey-predator model with infection in both prey and predator,” Filomat, vol. 29, no. 8, pp. 1753–1767, 2015, doi: Dynamical of Ratio-Dependent Eco-epidemiology Model with Prey Refuge Adin Lazuardy Firdiansyah 237 10.2298/FIL1508753B. [7] R. K. Upadhyay and P. Roy, β€œSpread of a disease and its effect on population dynamics in an eco-epidemiological system,” Commun. Nonlinear Sci. Numer. Simul., vol. 19, no. 12, pp. 4170–4184, 2014, doi: 10.1016/j.cnsns.2014.04.016. [8] D. Mukherjee, β€œHopf bifurcation in an eco-epidemic model,” Appl. Math. Comput., vol. 217, no. 5, pp. 2118–2124, 2010, doi: 10.1016/j.amc.2010.07.010. [9] P. J. Hudson, A. P. Dobson, and D. Newborn, β€œDo Parasites make Prey Vulnerable to Predation? Red Grouse and Parasites,” J. Anim. Ecol., vol. 61, no. 3, p. 681, 1992, doi: 10.2307/5623. [10] N. Apreutesei and G. Dimitriu, β€œOn a prey-predator reactiondiffusion system with Holling type III functional response,” J. Comput. Appl. Math., vol. 235, no. 2, pp. 366– 379, 2010, doi: 10.1016/j.cam.2010.05.040. [11] A. L. Firdiansyah, β€œEffect of Prey Refuge and Harvesting on Dynamics of Eco- epidemiological Model with Holling Type III,” vol. 3, no. 1, pp. 16–25, 2021. [12] A. K. Misra and B. Dubey, β€œA ratio-dependent predator-prey model with delay and harvesting,” J. Biol. Syst., vol. 18, no. 2, pp. 437–453, 2010, doi: 10.1142/S021833901000341X. [13] U. De Lausanne and S. Brook, β€œCoupling in predator-prey dynamics: ratio- dependence,” J. Theor. Biol., vol. 139, pp. 311–326, 1989. [14] M. Verma and A. K. Misra, β€œModeling the Effect of Prey Refuge on a Ratio-Dependent Predator–Prey System with the Allee Effect,” Bull. Math. Biol., vol. 80, no. 3, pp. 626– 656, 2018, doi: 10.1007/s11538-018-0394-6. [15] S. Wang, Z. Ma, and W. Wang, β€œDynamical behavior of a generalized eco- epidemiological system with prey refuge,” Adv. Differ. Equations, vol. 2018, no. 1, pp. 1–20, 2018, doi: 10.1186/s13662-018-1704-x.