Effect of Fear in Leslie-Gower Predator-Prey Model with Beddington-DeAngelis Functional Response Incorporating Prey Refuge

—In the present paper, we study the effect of an- tipredator behavior due to fear of predation on a modified Leslie- Gower predator-prey model incorporating prey refuge which predation rate of predators follows Beddington-DeAngelis functional response. The biological justification of the model is demonstrated through non-negativity, boundedness, and permanence. Next, we perform the analysis of equilibrium and local stability. We obtain four equilibrium points where two points are locally asymptotically stable and other points are unstable. Besides, we show the effect of the fear in the model and obtain a conclusion that the increased rate of fear can decrease the density of both populations, and prey populations become extinct. Meanwhile, for the case with a constant rate of fear, the prey refuge helpful to the existence of both populations. However, for the case with the fear effect is large, prey refuge cannot cause the extinction of predators. Several numerical simulations are performed to support our analytical results.


I. INTRODUCTION
I N their habitat, each individual always interacts with another individual of the same type or different types, either individuals in a population or individual from another population. Interaction between individuals that tends to increase or decrease the population densities depend on whether the interactions are beneficial or detrimental. Predation is one of the interactions between individuals, like predators and prey. In predation, predators eat prey as their source of food. Thus, the number of the predator is large, while the number of the prey is lower.
In the past few decades, the existence of interacting predators and prey is described in the predator-prey model. It is used to understand the long-time behavior of individuals. Most of the predator-prey model used by several researchers is the classical Lotka-Volterra model where the prey consumption rate of predator is the growth rate of predator with conversion factors [1]. This paper concentrates on the Leslie-Gower model where both populations are growing logistically with the environmental carrying capacity of predator is proportional to the number of prey, emphasizing the fact that there is an upper limit to the rate of increase in both populations, which is not recognized in the Lotka-Volterra model [2].
Many researchers study the Leslie-Gower model with its various modifications. The model has been modified to confirm A.L. Firdiansyah  actual conditions, such as the Alle effect, harvesting, etc. Aziz and Okiye [3] have modified the Leslie-Gower model and showed global stability of the interior equilibrium point in the modified Leslie-Gower model. Then, Gupta and Chandra [4] modified Aziz and Okiye's model by adding Michaelis-Menten type prey harvesting. Meanwhile, Cai et al. [5] also modified Aziz and Okiye's model by incorporating the Allee effect on prey and investigated the existence of positive equilibriums, stability, and Hopf bifurcation in the model. In the predator-prey model, predators who consume prey are denoted as a functional response which means that the predation rate of predators to prey per capita [6]. It is well that the functional response plays an important component in all predator-prey interactions. Holling [7] introduced three types of functional responses that depend on the density of prey, namely Holling type I, Holling type II, and Holling type III. It has been studied by many researchers. Base on the explanation of [8], the Holling type I means the predation of predator increases linearly with the number of prey but then suddenly until a constant value when predators are satiated, e.g. [9]. The Holling type II means the predation of predator increases when the number of prey decreases, e.g. [3]. Meanwhile, the Holling type III means the predation of predator increases when the number of prey is huge but the predation of predator decreases when the number of prey is less. In this function, predators focus to eat prey in a location where is more abundant, e.g. [10].
It is well known that three Holling types of functional response are modeled as a function of prey density only. However, this function is unrealistic because it neglects the interference among predators [8]. In some situations, predators have to compete or share food. By experiment, predators do interfere with each other's activities causing the competition effects and then the prey changes its behavior when the threat of predators increases [11]. Therefore, the predator densities are important to consider in the predatorprey model [6]. Beddington [12] and DeAngelis et al. [13] suggested a function that can describe the ecology system more reasonably. They proposed a function form that depends on predator species and the environment that protects the prey [6]. The function form is Beddington-DeAngelis functional response. Currently, this function has been studied by many researchers. As in [2], Yu has modified Aziz and Okiye's model by changing Holling type II into Beddington-DeAngelis functional response and investigated the stability of equilibrium points.
Predators can directly disrupt the ecological system by consuming prey and also can indirectly affect the behavior and psychology of prey [14]. Predators affect prey populations indirectly with the fear of predation [15]. Pangle et al. [16] found that indirect effects on prey population growth rate are larger than the direct effect of the predator. For example, the mule deer spends less time foraging because of the predation of mountain lions [17]. The bird populations respond to predatory sound by escaping from their nets in times of danger [15]. Currently, many studies explain that fear has a pretty strong impact on the ecology system. As in [18], their experiment showed that the female birds that experienced frequent nest predation produced fewer eggs in subsequent nests. Therefore, the presence of predators has a stronger impact on prey demographics than direct predation. Zanette et al. [19] have done an experiment on song sparrows which shows that 40 percent of the offspring produced by song sparrows (Melospiza melodia) are reduced due to fear of the predator.
This phenomenon of fear has been adapted by several mathematical studies. As in [20], Wang et al. proposed the predator-prey model by incorporating the fear effect. They obtain a conclusion that the fear effect can stabilize the predator-prey model by negating the existence of periodic solutions. In 2019, Pal et al. [1] observe the fear effect in the modified Leslie-Gower model where predators collaborate during hunting. They obtain a conclusion that fear becomes the stabilization of the model by negating the existence of periodic solutions and creates the more potent model. Pal et al. [11] observe the fear effect in the predator-prey model with Beddington-DeAngelis functional response. They find that fear plays an important component in the stability of an ecosystem. Then, Kundu et al. [21] investigate the fear effect in the discrete-time predator-prey model and obtain that the fear effect increases the stability of the model. Moreover, they find that the fear effect decreases the birth rate of prey.
Recently, hiding behavior has become an interesting topic in mathematical studies. Predation threats and limited resources can affect the migration of the prey population so that hiding behavior is an option for survival [22]. This hiding behavior gives some degree of protection for prey populations which can protect them from extinction. Therefore, prey refuge is a more relevant component of predator-prey dynamics and can prevent the extinction of prey populations [23].
The combination of prey refuge and fear effect has been discussed by many researchers. Wang et al. [24] observe the fear effect in the modified Leslie-Gower predator-prey model by incorporating prey refuge. They find that the fear effect is complex, namely (i) the fear effect can decrease the density of both population and prey populations become extinction, (ii) the effect of fear on the stability is rich and complex, (iii) the existence of prey and predator depends on prey refuge with a constant of fear. Meanwhile, Zhang et al. [25] also investigate the fear effect in the predatorprey model by adding prey refuge. They conclude that the fear effect does not just reduce the density of predator but also can stabilize the model by excluding the existence of periodic solutions. Moreover, the fear effect has a strong impact on the extinction of predators.
Based on the motivation above, we respect a modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response as in [2] by including the effect of fear and prey refuge. The primary discussion of our research is to observe the effect of fear in the model where the rate of fear has an effect on the stability of the model and makes the model is more potent. Several numerical simulations are performed to confirm our analytical results.

II. THE MATHEMATICAL MODEL
We respect the model as in [2] where x 1 (t) is the number of prey at time t and x 2 (t) is the number of predator at time t. The mathematical model is as follows: with x 1 (0) > 0 and x 2 (0) > 0. All parameters are positive values and their definitions are given as follows: r1(r2) is the birth rate of prey (or predators), p(β) is the competition rate of prey (or predators), α is the reduction rate of prey into predators, a(γ) is environmental protects prey (or predators), and b, c are appropriate constants.
In this paper, we assume that prey populations can refuge from predation. Prey populations that are completely protected from predation are denoted by parameter η. Meanwhile, (1 − η)x is the number of prey outside protection where η ∈ [0, 1) is a constant that measures the protection of prey. Thus, model (1) can be rewritten as follows: with x 1 (0) > 0 and x 2 (0) > 0. All parameters are positive values and their definitions are similar to model (1). According to [19], the fear effect can influence prey production. Thus, we incorporate the fear effect by multiplying the birth rate of prey with a factor f (k, x2) = 1 1+kx2 which means the impact of anti-predator response due to fear. Meanwhile, k is the rate of fear which expresses anti-predator behavior in prey. We assume that the frightened prey forages less and makes the young less protected so that decreases the birth rate of prey [26]. Biologically speaking, the fear effect satisfies several conditions as follows: The meaning of conditions can be shown in [22]. Hence, the model transforms into the following model.
with x 1 (0) > 0 and x 2 (0) > 0. All parameters are positive values and their definitions are similar to model (1). Before we begin to analyze the mathematical model, we set model (3) by applying the following scaling: Thus, model (3) becomes where The initial conditions are x(0) > 0 and y(0) > 0. All parameters are positive values.

A. Non-negativity
It is known that the model (4) consists of two interacting populations. Therefore, we have to show that the model solutions have non-negativity values. The non-negative solution is guaranteed by theorem as below.
Proof: From the model (4), we obtain that Because x(0), y(0) > 0 and the exponential function is nonnegative for any real number, then the model solutions remain non-negative, i.e., x(t) ≥ 0 and y(t) ≥ 0.

B. Boundedness
It is known that the populations have limited resources in the nature. So, the populations must be bounded. The boundedness for populations is showed by using the following Lemma as in [27].
Proof: From the first equation of the model (4), we have Assume that µ > δ. By using Lemma 1, we obtain For any sufficiently small ε 1 > 0, there is a T 1 ≥ 0 such that Based on the results above and the second equation, we get Assume that η < 1 and µ > δ. Then from Lemma 1, we obtain

C. Permanence
In this part, we present the permanence of the model (4). In the biological meaning, this guarantees that every population can survive for a long time. Analytically, the model (4)  Theorem 3: The model (4) is permanent when ω > 0.
Proof: From Theorem 2, for any sufficiently small ε 2 > 0, there exists T 2 ≥ 0 such that By using Theorem 2 and the first equation of the model (4), we get where Assume that ω > 0. Based on Lemma 1, we get

IV. THE EQUILIBRIUM AND STABILITY ANALYSIS
In this subsection, we investigate the dynamics of model (4) and recall several results as follows: A. The equilibrium points By setting dx dt = dy dt = 0, we get the equilibrium points of the model (4) as follows: • The trivial point of model (4) with It is well known that U 1 and U 2 are always positive. However, the signs of U 3 and U 4 are not obvious. Therefore, by using Descartes's rule of signs, we have several conditions for U 3 and U 4 as follows: • When U 3 > 0 and U 4 > 0, the equation (5) does not contain positive root because there is no change of sign. • When U 3 > 0 and U 4 < 0, the equation (5) contains one positive root.
• When U 3 < 0 and U 4 < 0, the equation (5) contains one positive root. • When U 3 < 0 and U 4 > 0, the equation (5) contains two positive roots. Based on the results above, if U 4 < 0, then the equation (5) contains at least one positive root. Thus, the necessary condition that has at least one positive root for the existence of E * is U 4 must be negative. Meanwhile, the explicit form of the root of equation (5) can be obtained by using Cardano's approach as in [28].

B. The Stability Analysis
Now, we shall analyze the stability of model (4) by investigating eigenvalues from the Jacobian matrix at E 0 (0, 0), E 1 (µ − δ, 0), and E 2 (0, ν θ ) which are respectively given as follows: Based on the above analytical results, we obtain several conditions for the stability of each equilibrium point. The equilibrium point E 0 and E 1 are unstable. The equilibrium point E 2 is locally asymptotically stable when θµ θ+ρν < δνξ+νϕ(1−η)+δθ θ+ξν . Meanwhile, the point E 2 is unstable when θµ θ+ρν > δνξ+νϕ(1−η)+δθ θ+ξν . To investigate the stability of the interior point E * (x * , y * ), we identify the Jacobian matrix at E * where φ ij is the entry of matrix with row i and column j as follows: The trace and determinant values of the Jacobian matrix at E * are given as follows:  (6) By using the Routh-Hurwitz criteria, the point E * is locally asymptotically stable if it satisfies condition tr(J E * ) < 0 and det(J E * ) > 0. Therefore, the condition of stability for E * is as below:
In Figure 1b, we set ρ = 0, η = 0.25 which means that model (4) does not include fear but only prey refuge. In this case, the result is the same as in Figure 1a. The difference of the previous figure is the point E * changes from (1.0894, 1.6991) to (1.1976, 1.426) which is affected by prey refuge.
In Figure 2a, we set ρ = 1, η = 0.25. In this case, model (4) includes the fear effect and prey refuge. The result of this case is similar to Figure 1b. However, the difference of the  (6) previous figure is the point E * changes from (1.1976, 1.426) to (0.4587, 0.6344). If we set ρ = 10, η = 0.25 as in Figure 2b, then the point E * changes from (0.4587, 0.6344) to (0.037, 0.1825). It is noted that the density of both populations decreases which is caused by the fear effect. Therefore, the rate of fear affects the density of both populations.
In Figure 3a, we set ρ = 10, η = 0.85 which means that model (4) includes the effect of fear and prey refuge where the number of prey outside prey refuge is less than the previous case. By comparing Figure 2b with Figure 3a, we obtain that E * changes from (0.037, 0.1825) to (0.1778, 0.181). It is noted that the value of x * of E * increases from 0.037 to 0.1778. Meanwhile, the value of y * of E * reduce from 0.1825 to 0.181. Therefore, the existence of prey and predator with a constant rate of fear needs large enough prey refuge. Thus, prey refuge helpful to the existence of both populations.
In Figure 3b, we set ρ = 100, η = 0.85 which means that model (4) includes prey refuge and the fear effect that is greater than the previous case. The result of this case shows that the point E 2 (0, 0.14286) is locally asymptotically stable and model (4) has no interior point E * . It is noted that the rate of the fear reduces the number of both populations where x * = 0 and y * > 0 forever. Therefore, the fear factor cannot cause the extinction of predators but only the extinction of the prey.
In Figure 4, we set ρ = 0, η = 0.85 which means that model (4) does not include fear but only prey refuge. The result of this case is the point E * is locally asymptotically stable and this is different from Figure 3b. It shows that the prey refuge cannot cause the extinction of the predator.  (6) VI. CONCLUSIONS In this paper, we have investigated the fear effect in the modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response including prey refuge. We find four equilibrium points exist under certain conditions where two equilibrium points are locally asymptotically stable, namely E 2 and E * and other points are unstable. Biologically, we conclude that the effect of fear and prey refuge influences the model (4). The fear effect decreases the number of both populations. If the fear effect is getting large, then the number of both populations is lower, and prey populations become extinct. These results are similar to [24], [29]. For the case with a constant of fear, the existence of prey and predator needs large enough prey refuge. This result is similar to [24]. However, for the case with the fear effect is large, prey refuge cannot cause the extinction of predators. For future work, we can observe the global stability of the interior point on the model (4).