Dynamical Behaviors of a Modified Leslie-Gower Predator-Prey System with Fear Effect and Prey Refuge

Ke Yuan

School of Mathematics and Statistics, Fujian Normal University, Fuzhou, Fujian, China.

DOI: 10.4236/ojmsi.2024.124011                

Abstract

In this paper, the dynamical behaviors of a modified Leslie-Gower predator-prey model incorporating fear effect and prey refuge are investigated. We delve into the construction of the model and its biological significance, with preliminary results encompassing positivity, boundedness, and persistence. The stability of the system’s boundary and positive equilibrium points is proven by calculating the real part of the eigenvalues of the Jacobian matrix. At the positive equilibrium point, we demonstrate that the system’s unique positive equilibrium is globally asymptotically stable by using the Dulac criterion. Furthermore, at this equilibrium point, we employ the Implicit Function Theorem to discuss how fear effects and prey refuges influence the population densities of both prey and predators. Finally, numerical simulations are conducted to validate the above-mentioned conclusions and explored the impact of Predator-taxis sensitivity $\alpha$ on dynamics of the system.

Keywords

Fear Effect, Prey Refuge, Predator-Taxis Sensitivity, Population Density, Stability

1. Introduction

In ecology, the study on the dynamic behavior of predator-prey models has always been a hot topic. The primary means to advance the theoretical work of predator-prey models is to construct mathematical models that are more practically meaningful, which represents a significant application of mathematics in ecology. In this field, Lotka [1] and Volterra [2] proposed the renowned Lotka-Volterra model. This model illustrates that under conditions of abundant food, the population densities of predators and prey will grow exponentially. Considering the limited environmental resources, the number of species cannot grow indefinitely. Leslie [3] and Gower [4] introduced the classic Leslie-Gower model, which emphasizes that both predators and prey have upper limits to their growth, and the carrying capacity of predators is proportional to the number of prey. However, in nature, many predators do not solely rely on a single food source; in the absence of their preferred food, predators will hunt other food sources for survival. Thus, the carrying capacity of predators cannot fully reflect this reality. In particular, Aziz-Alaoui and Okiye [5] proposed a improved Leslie-Gower predator-prey model. The modified model better reflects biological facts, as predators can survive due to the presence of substitute prey even if their primary prey becomes extinct. Recently, numerous scholars had investigated predator-prey models with Leslie-Gower functional response functions [6]-[8].

In the intricate dynamics of predator-prey interactions, two crucial factors significantly influence the population size of prey species. Firstly, predators directly impact prey populations through predation. Secondly, the mere presence and indirect behaviors of predators also exert an influence on prey populations. In 1998, theoretical biologist Lima [9] argued that the presence of predators elicits fear in prey, leading to alterations in their behavioral and physiological traits, which subsequently affects prey population sizes. This effect, under certain conditions, can far outweigh the direct consequences of predation. When confronted with the risk of predators, prey exhibit a series of behavioral changes, such as modifications in habitat selection and reductions in foraging and breeding activities [10]-[12]. Notably, Zanette et al. [12] demonstrated the profound impact of predator-induced fear on prey populations during the breeding season of song sparrows. By manipulating predator sounds without direct predation, they observed a 40% decline in the fecundity of prey populations, thereby validating the significant influence of predator fear on prey populations. However, the implications of fear extend beyond these findings; at heightened levels, fear can significantly compromise the physical condition of juvenile prey and have detrimental effects on the survival of adult prey [13] [14].

To investigate the influence of fear effects, Wang et al. [15] initially proposed a predator-prey model incorporating fear effects, demonstrating that a high level of fear can stabilize the predator-prey system by excluding the existence of periodic solutions, whereas a relatively low level of fear can induce multiple limit cycles, leading to a bistable phenomenon. Furthermore, various models incorporating fear effects had been extensively studied [16]-[20]. In [16], Das et al. analyzed a predator-prey model with a Holling II functional response, which simultaneously considers the influence of fear effects on both the birth rate and mortality rate of the prey population. Waqas et al. [17] introduced a fear factor into a discrete-time predator-prey model to investigate its impact on the complexity and dynamic behavior of the system. This study primarily focused on assessing the local stability of the trivial and boundary equilibrium points, and exploring the criteria for Neimark-Sacker and period-doubling bifurcations. Sarkar et al. [18] improved the expression of fear functions by incorporating a minimum fear cost. Based on [18], Dong et al. [19] proposed a fear function with saturated fear cost and studied a predator-prey model that incorporates both saturated fear cost and Predator-taxis sensitivity to prey attraction. Their model aimed to analyze how these factors jointly affect the dynamics and stability of the predator-prey interactions.

In nature, not all prey are captured by predators, as prey seek refuges to evade predation [21] [22]. For prey, finding shelters constitutes a low-cost anti-predation behavior that directly reduces encounters with predators, thereby mitigating the risk of extinction and contributing to the stability of the system [23] [24]. Two common forms of refuge are proportional refuge and constant refuge. Research by J. Ghosh et al. [25] has demonstrated that a certain degree of refuge can facilitate the coexistence of all populations, while a refuge proportion exceeding a certain threshold may lead to the extinction of predators. Li et al. [26] have investigated the dynamical behaviors of a discrete predator-prey system incorporating the fear effect and refuges. By employing the Center Manifold Theorem and bifurcation theory, they calculated the critical coefficients for Flip and Hopf bifurcations, as well as the properties of the bifurcating periodic solutions. To cope with the risk of predation, prey populations defend themselves against predators by seeking shelters. Nevertheless, prey’s energy for foraging and refuge-seeking is limited. Devoting more energy to finding shelters would inevitably lead to a reduction in reproduction due to such investment. Consequently, prey refuge entails both costs and benefits. Here, the cost refers to the decrease in the prey population’s birth rate caused by seeking refuge, while the benefit stems from the reduced risk of being killed by predators. In most predator-prey models, the costs and benefits of refuges are treated as independent factors. However, given the finitude of time and energy, they might not be independent of each other [27]. Hence, we correlate these two factors by introducing a Predator-taxis sensitivity parameter $\alpha$ . As the Predator-taxis sensitivity increases, prey populations allocate more energy to seeking shelters, resulting in a decreased predation success rate. Concurrently, with less time for foraging, the reproduction rate of the prey population also declines.

Based on the aforementioned research backdrop, we incorporate the fear effect and proportional refuge into a modified Leslie-Gower predator-prey model, and intertwine these two independent factors through the Predator-taxis sensitivity parameter $\alpha$ . The subsequent sections of this paper are organized as follows: In Section 2, we meticulously construct the model and delve into its positivity, boundedness, and permanence. In Section 3, we present the existence and stability of various equilibrium points within the model. In Section 4, we investigate the influence of fear effect and prey refuge on population density. In Section 5, we leverage MATLAB software to conduct numerical simulations that validate our preceding analytical findings. These simulations provide a visual representation of the model’s dynamical behaviors, reinforcing the theoretical results. The conclusion is given in Section 6.

2. Formulation of Mathematical Model

Leslie and Gower proposed the following Leslie-Gower [28] predator-prey model.

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx-dx-b{x}^2-p\left(x\right)y,\hfill \\ \frac{\text{d}y}{\text{d}t}=sy\left(1-\frac{y}{ax}\right),\hfill \end{array}$ (2.1)

where $x\left(t\right)$ , $y\left(t\right)$ represent the density of prey and predator at the time $t$ , respectively; The parameter $r$ represents the birth rate of the prey, $d$ denotes the death rate of the prey, $b$ is the intraspecific competition coefficient among prey individuals, and $s$ signifies the intrinsic growth rate of the predator population. The predator preys upon the prey according to a functional response $p\left(x\right)$ and an environmental carrying capacity $ax$ , where $a$ captures the efficiency of converting prey food quality into predator births. All parameters $r,d,b,s,a$ are assumed to be strictly positive constants. Additionally, we postulate that the intrinsic growth rate of the prey population is greater than zero, i.e., $r>d$ , as this condition is necessary for the persistence of the prey population and to avoid its extinction.

In the presence of predators, prey populations engage in anti-predator behaviors such as refuge-seeking to cope with the risk of predation. In this paper, we consider Holling I functional response with proportional refuges, given by $p\left(x\right)=\beta \left(1-m\right)x$ , where the parameter $\beta$ represents the capture rate of the predator, and $m\in \left(0,1\right)$ indicating the refuge coefficient. The fear of predators can reduce the birth rate of prey populations and also affect their physical conditions, yet it does not lead to the extinction of prey populations. Even if the fear level is sufficiently high, prey can survive under saturated fear costs [19]. Therefore, we incorporate the fear function $f\left(k,y\right)=\eta +\frac{1-\eta }{1+ky}$ , where $\eta \in \left(0,1\right)$ represents the saturated fear cost and $k$ denotes the fear level. We consider that the fear effect simultaneously impacts the birth rate and mortality rate of prey. In times of food scarcity, predators will turn to capturing other prey to survive. Thus, our model accounts for the predators having substitute prey. Consequently, we consider the following model:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(\eta +\frac{1-\eta }{1+ky}\right)-d{x}^2-\beta \left(1-m\right)xy,\hfill \\ \frac{\text{d}y}{\text{d}t}=sy\left(1-\frac{y}{a\left(1-m\right)x+c}\right),\hfill \end{array}$ (2.2)

where $r$ is the intrinsic growth rate of prey, $d$ represents death due to intra-prey competion, and $c$ represents the number of other substitute prey. In the aforementioned system (2.2), the costs and benefits of the refuge are treated as mutually independent. However, as mentioned in the introduction, they may not be independent of each other. This is because the prey’s energy allocated for seeking refuge and reproduction is limited. If too much energy is invested in seeking refuge, reproduction may subsequently decrease. We introduce the Predator-taxis sensitivity $\alpha$ to account for this relationship. Therefore, we modify the function

$p\left(x\right)$ and $f\left(k,y\right)$ in (2.2), as $f\left(k,y\right)\equiv f\left(k,\alpha ,y\right)=\eta +\frac{1-\eta }{1+k\alpha y}$ , and $p\left(x\right)\equiv p\left(\alpha ,x\right)=\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x$ . Consequently, the form of the system considered in this paper is as following:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(\eta +\frac{1-\eta }{1+k\alpha y}\right)-d{x}^2-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)xy=F\left(x,y\right),\hfill \\ \frac{\text{d}y}{\text{d}t}=sy\left(1-\frac{y}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c}\right)=G\left(x,y\right).\hfill \end{array}$ (2.3)

A description of the parameters is given in Table 1.

Table 1. The model parameters and variables are described as follows.

Parameter	Description
$x$	Density of prey population
$y$	Density of predator population
$r$	Growth rate of the prey population
$d$	Intraspecific competition coefficient among prey population
$\beta$	Rate of predation
$a$	Food quality of prey for conversion into predator’s growth
$m$	Refuge coefficient
$\eta$	Saturated fear cost
$k$	Level of fear
$c$	Substitute prey of predator
$s$	Growth rate of the predator population
$\alpha$	Predator-taxis sensitivity

Here, the function $f\left(k,\alpha ,y\right)=\eta +\frac{1-\eta }{1+k\alpha y}$ and $p\left(\alpha ,x\right)=\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x$ have the following properties.

1) $f\left(0,\alpha ,y\right)=1$ , $f\left(k,0,y\right)=1$ , $f\left(k,\alpha ,0\right)=1$ : implies that when prey are insensitive to predators or predators are absent, the prey population will not exhibit any anti-predator behavior, and the fear function will have no impact on the birth rate of the prey.

2) ${\lim }_{k\to \infty }f\left(k,\alpha ,y\right)=\eta$ , ${\lim }_{\alpha \to \infty }f\left(k,\alpha ,y\right)=\eta$ , ${\lim }_{y\to \infty }f\left(k,\alpha ,y\right)=\eta$ : means that the fear effect will not lead to the extinction of the prey population. Even in high-fear environments and with abundant predator numbers, prey can survive through various anti-predator behaviors, such as seeking refuge. At this point, the growth rate of the prey population tends to be $r\eta$ .

3) $\frac{\partial f\left(k,\alpha ,y\right)}{\partial k}<0$ , $\frac{\partial f\left(k,\alpha ,y\right)}{\partial \alpha }<0$ , $\frac{\partial f\left(k,\alpha ,y\right)}{\partial y}<0$ : means that the growth rate of the prey population decreases with increasing levels of fear, Predator-taxis sensitivity, and the number of predators.

4) $p\left(0,x\right)=\beta x$ : implies that when prey are insensitive to predators, predators consume prey at their maximum rate, which is denoted by the efficiency parameter $\beta x$ .

5) ${\lim }_{\alpha \to \infty }p\left(\alpha ,x\right)=\beta \left(1-m\right)x$ : implies that when prey become highly sensitive to predators, they maximize their efforts to seek refuge. Under such circumstances, the availability of refuge in the environment reaches its saturation point.

6) $\frac{\partial p\left(\alpha ,x\right)}{\partial \alpha }<0$ : indicates that the predator’s capture rate decreases as the Predator-taxis sensitivity increases.

Well-Posedness

In this subsection, we delve deeper into the positivity, boundedness, and permanence of the system (2.3). This analysis serves as the cornerstone for understanding the fundamental dynamical properties of the system.

Theorem 2.1. Every solution of system (2.3) with initial condition $x\left(0\right)>0$ , $y\left(0\right)>0$ is positive and ultimately bounded, and the bounded positive invariant set is $\Omega =\left\{\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2|0<x<\frac{r}{d},0<y<K\right\}$ , where $K=\frac{ra}{d}\left(1-\frac{\alpha m}{1+\alpha }\right)+c$ .

Proof. Since the right-hand side of system (2.3), $F\left(x,y\right)$ and $G\left(x,y\right)$ are completely continuous and locally Lipschitzian on ${ℝ}_{+}^2$ , the solution $\left(x\left(t\right),y\left(t\right)\right)$ of system (2.3) with initial conition $\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ exists and is unique on ${ℝ}_{+}^2$ .

From system (2.3), we can easily obtain that:

$\begin{array}{l}x\left(t\right)=x\left(0\right)\text{exp}\left\{{\displaystyle {\int }_0^t\left[r\left(\eta +\frac{1-\eta }{1+k\alpha y\left(u\right)}\right)-dx\left(u\right)-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y\left(u\right)\right]\text{d}u}\right\},\\ y\left(t\right)=y\left(0\right)\text{exp}\left\{{\displaystyle {\int }_0^t\left[s\left(1-\frac{y\left(u\right)}{a\left(1-\frac{\alpha m}{1+\alpha }x\left(u\right)\right)+c}\right)\right]\text{d}u}\right\}.\end{array}$

This shows that all solutions of system (2.3) with positive initial values remain positive.

Next, we will prove the boundedness of the system (2.3).

From the first equation of (2.3) we can get

$\frac{\text{d}x}{\text{d}t}\le rx-d{x}^2.$

Applying Lemma 1 in [29] to the differential inequality, then

$\underset{t\to \infty }{\lim }\sup x\left(t\right)\le \frac{r}{d}.$

Therefore, for any $ϵ>0$ , there exists $T>0$ such that $x\left(t\right)\le \frac{r}{d}+ϵ$ for $t>T$ . When $t>T$ , From the second equation of system (2.3), for large $t>0$ , we have

$\frac{\text{d}y}{\text{d}t}\le sy\left(1-\frac{y}{a\left(1-\frac{\alpha m}{1+\alpha }\right)\left(\frac{r}{d}+ϵ\right)+c}\right),t>T.$

Then

$\underset{t\to \infty }{\lim }\sup y\left(t\right)\le \frac{ra}{d}\left(1-\frac{\alpha m}{1+\alpha }\right)+c.$

Thus, system (2.3) is bounded and the bounded set is

$\Omega =\left\{\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2|0<x<\frac{r}{d},0<y<K\right\}$ .

The investigation of the system’s permanence holds significant importance, as it signifies the ability of both species to sustainably coexist and avoid extinction in their mutual interactions over time.

System (2.3) is uniformly permanent is equivalent to there exists $ϵ>0$ , for every solution of system is non-negative such that ${\lim }_{t\to \infty }\inf x\left(t\right)\ge ϵ$ , ${\lim }_{t\to \infty }\inf y\left(t\right)\ge ϵ$ .

Theorem 2.2. If (H1) $r\left(1+\eta k\alpha K\right)-\beta K\left(1+k\alpha K\right)>0$ holds, then system (2.3) is permanent.

Proof. From Theorem 2.1 we can get for sufficiently large $t>0$

$\begin{array}{c}\frac{\text{d}x}{\text{d}t}=rx\left(\eta +\frac{1-\eta }{1+k\alpha y}\right)-d{x}^2-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)xy\\ \ge rx\left(\eta +\frac{1-\eta }{1+k\alpha K}\right)-d{x}^2-\beta xK\\ =\left[r\left(\eta +\frac{1-\eta }{1+k\alpha K}\right)-\beta K\right]x-d{x}^2,\end{array}$

and

$\frac{\text{d}y}{\text{d}t}=sy\left(1-\frac{y}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c}\right)\ge sy\left(1-\frac{y}{c}\right).$

Applying Lemma 1 in [29] to the above differential inequality, combine with (H1) then

$\begin{array}{l}\underset{t\to \infty }{\lim }\inf x\left(t\right)\ge \frac{r\left(1+\eta k\alpha K\right)-\beta K\left(1+k\alpha K\right)}{d\left(1+k\alpha K\right)},\underset{t\to \infty }{\lim }\inf y\left(t\right)\ge c.\hfill \end{array}$

Choose $ϵ=\min \left\{\frac{r\left(1+\eta k\alpha K\right)-\beta K\left(1+k\alpha K\right)}{d\left(1+k\alpha K\right)},c\right\}$ , combine with Theorem 2.1, we have completed the proof of the permanence of system (2.3).

3. Equilibria and Their Stability

In this Section, we investigate the existence and stability of various equilibrium points within the system. This examination is pivotal in predicting the long-term behavior of the predator-prey interaction under different ecological scenarios.

3.1. Existence Analysis of Equilibria

It is easy to get that system (2.3) has the following three equilibrium points:

i) The trivial equilibrium $E_0\left(0,0\right)$ .

ii) The predator-free equilibrium $E_1\left(\frac{r}{d},0\right),$ the prey-free equilibrium $E_2\left(0,c\right)$ .

Now, we analyze the existence of the positive equilibrium $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ , where the expressions of $x^{\text{*}}$ and $y^{\text{*}}$ can be obtained by the following equation.

$\{\begin{array}{l}r\left(\eta +\frac{1-\eta }{1+k\alpha y}\right)-dx-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y=0,\hfill \\ s\left(1-\frac{y}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c}\right)=0.\hfill \end{array}$ (3.1)

From the second equation of Equation (3.1), it is easy to know

$y^{\text{*}}=a\left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}+c.$

Substituting the expression of $y^{\text{*}}$ into the first equation of Equation (3.1) yields the following equation.

$A{x}^2+Bx+C=0$ (3.2)

where

$\begin{array}{l}A=ka\alpha \theta \left(d+a\beta {\theta }^2\right)>0,\\ B=\left[d\left(1+kc\alpha \right)+a\beta {\theta }^2\right]-ka\alpha \theta \left(r\eta -2c\beta \theta \right),\\ C=c\beta \theta \left(1+kc\alpha \right)-r\left(1+kc\alpha \eta \right).\end{array}$

For the convenience of calculation, denote $\theta =1-\frac{\alpha m}{1+\alpha }$ , according to the expression of $\theta$ , it can be inferred that $\theta \in \left(0,1\right)$ .

Next, we can apply the Descartes’ rule of signs to determine the positive solution of Equation (3.2). If $C>0$ , we can get

$r<\frac{c\beta \theta \left(1+kc\alpha \right)}{1+kc\alpha \eta }.$ (3.3)

Note that $\eta \in \left(0,1\right)$ , we can multiply both sides of Equation (3.3) by $\eta$ and then reduce $2c\beta \theta$ , which results in

$\begin{array}{l}r\eta -2c\beta \theta <\frac{c\beta \theta \eta \left(1+kc\alpha \right)}{1+kc\alpha \eta }-2c\beta \theta \\ \iff r\eta -2c\beta \theta <\frac{c\beta \theta \left(\eta -2\right)-k\alpha \beta \theta \eta c^2}{1+kc\alpha \eta }<0\\ \iff B>0.\end{array}$

In this case, Equation (3.2) has no positive roots. To have a root for Equation (3.2), assume that $C<0$ , we can get the Equation (3.2) has a unique non-negative root

$x^{\text{*}}=\frac{-B+\sqrt{B^2-4AC}}{2A}.$

Theorem 3.1. If $C<0$ holds, the system (2.3) has a unique coexistence equilibrium $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ .

3.2. Stability Analysis

The stability of the equilibrium point is determined by the real parts of eigenvalues of the Jacobian matrix. The Jacobian matrix of system (2.3) is given by

$J_E=\left[\begin{array}{ll}{F}_x\hfill & F_y\hfill \\ G_x\hfill & G_y\hfill \end{array}\right]$ , where

$F_x=r\left(\eta +\frac{1-\eta }{1+k\alpha y}\right)-2dx-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y,$

$F_y=\frac{-kr\alpha \left(1-\eta \right)x}{{\left(1+k\alpha y\right)}^2}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x,$

$G_x=\frac{as{y}^2\left(1-\frac{\alpha m}{1+\alpha }\right)}{{\left[a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c\right]}^2},$

$G_y=s-\frac{2sy}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c}.$

Theorem 3.2. The trivial equilibrium $E_0\left(0,0\right)$ is an unstable node, $E_1\left(\frac{r}{d},0\right)$ is a saddle point, $E_2\left(0,c\right)$ is also a saddle point if $r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)-c\beta \left(1-\frac{\alpha m}{1+\alpha }\right)>0,$ otherwise, $E_2$ is a stable node.

Proof. The Jacobian matrix of system (2.3) at points $E_0\left(0,0\right)$ , $E_1\left(\frac{r}{d},0\right)$ , $E_2\left(0,c\right)$ are respectively

$J\left(E_0\right)=\left[\begin{array}{cc}r& 0\\ 0& s\end{array}\right],J_{\left(E_1\right)}=\left[\begin{array}{cc}-r& \frac{-r}{d}\left[kr\alpha \left(1-\eta \right)+\beta \left(1-\frac{\alpha m}{1+\alpha }\right)\right]\\ 0& s\end{array}\right],$

$J_{\left(E_2\right)}=\left[\begin{array}{cc}r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)-c\beta \left(1-\frac{\alpha m}{1+\alpha }\right)& 0\\ as\left(1-\frac{\alpha m}{1+\alpha }\right)& -s\end{array}\right].$

It is obvious that $E_0\left(0,0\right)$ is an unstable node, $E_1\left(\frac{r}{d},0\right)$ is a saddle point. The eigenvalues of $E_2\left(0,c\right)$ are $r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)-c\beta \left(1-\frac{\alpha m}{1+\alpha }\right)$ and $-s$ , when $r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)-c\beta \left(1-\frac{\alpha m}{1+\alpha }\right)>0$ , $E_2\left(0,c\right)$ is a saddle point, otherwise, $E_2\left(0,c\right)$ is stable node.

Similarly, the Jacobian matrix at the coexistence equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ as follow:

$J\left(E^{\text{*}}\right)=\left[\begin{array}{cc}-d{x}^{\text{*}}& \frac{-kr\alpha \left(1-\eta \right)x^{\text{*}}}{{\left(1+k\alpha y^{\text{*}}\right)}^2}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}\\ as\left(1-\frac{\alpha m}{1+\alpha }\right)& -s\end{array}\right].$

The characteristic equation is

${\lambda }^2-Tr\left(J_{E^{\text{*}}}\right)\lambda +D\left(J_{E^{\text{*}}}\right)=0,$

where

$Tr\left(J_{E^{\text{*}}}\right)=-d{x}^{\text{*}}-s<0,$

$D\left(J_{E^{\text{*}}}\right)=ds{x}^{\text{*}}+as\left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}\left(\frac{kr\alpha \left(1-\eta \right)}{{\left(1+k\alpha y^{\text{*}}\right)}^2}+\beta \left(1-\frac{\alpha m}{1+\alpha }\right)\right)>0.$

Hence, the coexistence equilibrium $E^{\text{*}}$ is local asymptotic stable.

Next, to demonstrate the global asymptotic stability of $E^{\text{*}}$ ,it is important to prove the absence of periodic orbits in the region Ω for system (2.3) by Bendison-Dulac criteria [30].

We assume Dulac function $B\left(x,y\right)=\frac{1}{xy},x>0,y>0$ , then

$\begin{array}{c}\text{Δ}=\frac{\partial \left(B\left(x,y\right)F_1\left(x,y\right)\right)}{\partial x}+\frac{\partial \left(B\left(x,y\right)G_1\left(x,y\right)\right)}{\partial y}\\ =-\frac{d}{y}-\frac{s}{x}\frac{1}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x+c}.\end{array}$

It is obvious that $\text{Δ}<0$ holds within Ω, under the condition of $C<0$ , we can obtain the following inequality

$\begin{array}{l}C<0\iff c\beta \theta \left(1+kc\alpha \right)<r\left(1+kc\alpha \eta \right)\\ \iff c\beta \theta <\frac{r\left(1+kc\alpha \eta \right)}{\left(1+kc\alpha \right)}=\frac{r\eta \left(\frac{1}{\eta }+kc\alpha \right)}{\left(1+kc\alpha \right)}<r\eta <r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)\\ \iff r\left(\eta +\frac{1-\eta }{1+kc\alpha }\right)-c\beta \left(1-\frac{\alpha m}{1+\alpha }\right)>0.\end{array}$

In other words, when $C<0$ , the $E_2\left(0,c\right)$ is a saddle point, and the system does not admit any closed orbits. Consequently, the unique coexistence equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ is globally asymptotically stable.

In conclusion, we have the following result.

Theorem 3.3. If $C<0$ holds, the system (2.3) has a unique coexistence equilibrium $E^{\text{*}}$ , and global asymptotic stable.

The global stability of the unique coexistence equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ in system (2.3) will be intuitively presented in Section 5 through numerical simulation.

4. Population Density Analysis

In this section, we study the effects of fear level $k$ and refuge coefficient $m$ on the coexistence equilibrium $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ where the expressions of $x^{\text{*}}$ and $y^{\text{*}}$ is satisfy the following equation:

$\{\begin{array}{l}r{x}^{\text{*}}\left(\eta +\frac{1-\eta }{1+k\alpha y^{\text{*}}}\right)-d{x}^{\text{*}2}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}{y}^{\text{*}}=0,\hfill \\ s{y}^{\text{*}}\left(1-\frac{y^{\text{*}}}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}+c}\right)=0.\hfill \end{array}$

Denote $P=r{x}^{\text{*}}\left(\eta +\frac{1-\eta }{1+k\alpha y^{\text{*}}}\right)-d{x}^{\text{*}2}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}{y}^{\text{*}}$ , $Q=s{y}^{\text{*}}\left(1-\frac{y^{\text{*}}}{a\left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}+c}\right)$ , the determinant of Jacobian matrix at $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ of system (2.3) is given by

$\begin{array}{c}J=\frac{D\left(P,Q\right)}{D\left(x^{\text{*}},y^{\text{*}}\right)}=\left|\begin{array}{cc}{P}_{x^{\text{*}}}& P_{y^{\text{*}}}\\ Q_{x^{\text{*}}}& Q_{y^{\text{*}}}\end{array}\right|\\ =\left[\begin{array}{cc}-d{x}^{\text{*}}& \frac{-kr\alpha \left(1-\eta \right)x^{\text{*}}}{{\left(1+k\alpha y^{\text{*}}\right)}^2}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}\\ as\left(1-\frac{\alpha m}{1+\alpha }\right)& -s\end{array}\right]\\ =ds{x}^{\text{*}}+as\left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}\left(\frac{kr\alpha \left(1-\eta \right)}{{\left(1+k\alpha y^{\text{*}}\right)}^2}+\beta \left(1-\frac{\alpha m}{1+\alpha }\right)\right)>0.\end{array}$

The proof here is based on the following implicit function existence theorem.

Theorem 4.1 (implicit function existence theorem [31])

Suppose the functions $F\left(x,y,\cdots ,z,u,v,\cdots ,w\right)=0$ and $G\left(x,y,\cdots ,z,u,v,\cdots ,w\right)=0$ (where $x,y,\cdots ,z$ are independent variables, and $u,v,\cdots ,w$ are dependent variables or implicit functions) are continuous in some neighborhood of the point $\left(x_0,y_0,\cdots ,z_0,u_0,v_0,\cdots ,w_0\right)$ , and their partial derivatives with respect to $u,v,\cdots ,w$ exist. Furthermore, assume that at this point, the determinant of the Jacobian matrix.

$J=\left(\begin{array}{cccc}\frac{\partial F}{\partial u}& \frac{\partial F}{\partial v}& \cdots & \frac{\partial F}{\partial w}\\ \frac{\partial G}{\partial u}& \frac{\partial G}{\partial v}& \cdots & \frac{\partial G}{\partial w}\end{array}\right)$

is non-zero, i.e., $\text{det}\left(J\right)\ne 0$ .

Then, in some neighborhood of the point $\left(x_0,y_0,\cdots ,z_0\right)$ , the equations $F=0$ and $G=0$ uniquely determine a set of implicit functions $u=u\left(x,y,\cdots ,z\right)$ , $v=v\left(x,y,\cdots ,z\right)$ , $\cdots$ , $w=w\left(x,y,\cdots ,z\right)$ , such that these implicit functions are continuous at the point $\left(x_0,y_0,\cdots ,z_0\right)$ and their partial derivatives with respect to $x,y,\cdots ,z$ exist.

Furthermore, these partial derivatives can be calculated using the following formulas:

$\frac{\partial u}{\partial x}=-\frac{J_1}{\text{det}\left(J\right)},\frac{\partial u}{\partial y}=-\frac{J_2}{\text{det}\left(J\right)},\cdots ,\frac{\partial w}{\partial z}=-\frac{J_n}{\text{det}\left(J\right)},$

where $J_1,J_2,\cdots ,J_n$ are the determinants obtained by replacing the corresponding columns of the Jacobian matrix $J$ with $\left(\frac{\partial F}{\partial x},\frac{\partial F}{\partial y},\cdots ,\frac{\partial F}{\partial z}\right)$ and $\left(\frac{\partial G}{\partial x},\frac{\partial G}{\partial y},\cdots ,\frac{\partial G}{\partial z}\right)$ , respectively.

By applying the theorem of implicit function, we can derive the derivatives of $x^{\text{*}}$ and $y^{\text{*}}$ with fear level $k$ :

$\begin{array}{c}\frac{\text{d}{x}^{\text{*}}}{\text{d}k}=-\frac{1}{J}\frac{D\left(P,Q\right)}{D\left(k,y^{\text{*}}\right)}=-\frac{1}{J}\left|\begin{array}{cc}{P}_k& P_{y^{\text{*}}}\\ Q_k& Q_{y^{\text{*}}}\end{array}\right|\\ =-\frac{1}{J}\frac{sr\alpha \left(1-\eta \right)x^{\text{*}}{y}^{\text{*}}}{{\left(1+k\alpha y^{\text{*}}\right)}^2}<0,\end{array}$

$\begin{array}{c}\frac{\text{d}{y}^{\text{*}}}{\text{d}k}=-\frac{1}{J}\frac{D\left(P,Q\right)}{D\left(x^{\text{*}},k\right)}=-\frac{1}{J}\left|\begin{array}{cc}{P}_{x^{\text{*}}}& P_k\\ Q_{x^{\text{*}}}& Q_k\end{array}\right|\\ =-\frac{1}{J}\frac{r\alpha \left(1-\eta \right)x^{\text{*}}{y}^{\text{*}}}{{\left(1+k\alpha y^{\text{*}}\right)}^2}\cdot as\left(1-\frac{\alpha m}{1+\alpha }\right)<0.\end{array}$

Thus, it follows that, the population densities of both prey and predators are inversely related to the parameter $k$ , suggesting that the fear effect negatively impacts their respective population sizes.

Next, we delve into the impact of prey refuges on the population densities of both prey and predators. Similarly, it can be concluded that:

$\begin{array}{c}\frac{\text{d}{x}^{\text{*}}}{\text{d}m}=-\frac{1}{J}\frac{D\left(P,Q\right)}{D\left(m,y^{\text{*}}\right)}=-\frac{1}{J}\left|\begin{array}{cc}{P}_m& P_{y^{\text{*}}}\\ Q_m& Q_{y^{\text{*}}}\end{array}\right|\\ =\frac{1}{J}\left(\frac{S\alpha \beta x^{\text{*}}{y}^{\text{*}}}{1+\alpha }+\frac{as\alpha x^{\text{*}}}{1+\alpha }\left(\frac{kr\alpha \left(1-\eta \right)x^{\text{*}}}{{\left(1+k\alpha y^{\text{*}}\right)}^2}+\beta \left(1-\frac{\alpha m}{1+\alpha }\right)x^{\text{*}}\right)\right)>0,\end{array}$

$\begin{array}{c}\frac{\text{d}{y}^{\text{*}}}{\text{d}m}=-\frac{1}{J}\frac{D\left(P,Q\right)}{D\left(x^{\text{*}},m\right)}=-\frac{1}{J}\left|\begin{array}{cc}{P}_{x^{\text{*}}}& P_m\\ Q_{x^{\text{*}}}& Q_m\end{array}\right|\\ =-\frac{1}{J}\left(\frac{as\alpha x^{\text{*}}}{1+\alpha }\left(d{x}^{\text{*}}-\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y^{\text{*}}\right)\right).\end{array}$

From the above that, $\frac{d{x}^{\text{*}}}{dm}>0$ can be inferred that the prey population density is an increasing function of $m$ , indicating that the prey refuge contributes positively to the growth of prey population density.

When $d{x}^{\text{*}}<\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y^{\text{*}}$ , $\frac{d{y}^{\text{*}}}{dm}>0$ , suggesting that the prey refuge is beneficial for increasing the predator population density under these conditions. Conversely, when $d{x}^{\text{*}}>\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y^{\text{*}}$ , $\frac{d{y}^{\text{*}}}{dm}<0$ , which means that the prey refuge is not conducive to the growth predator population density in this scenario. Specially, when $d{x}^{\text{*}}=\beta \left(1-\frac{\alpha m}{1+\alpha }\right)y^{\text{*}}$ , the predator population density attains its maximum value.

5. Numerical Simulations

In this section, we will validate our previous analysis by numerical simulations. All parameter values in this section are based on their biological significance in reality.

Example 5.1. To visualize the influence of fear level $k$ on the population densitives of both prey and predators, we denote the following parametric values:

$r=0.5,\eta =0.2,\alpha =0.25,d=0.1,\beta =0.5,m=0.3,s=0.2,a=0.5,c=0.1,$ (5.1)

and vary the parameter $k$ .

Under the condition of parameter (5.1), the inequality $C<0$ always holds at all times. Indicating that the system maintains a unique coexistence equilibrium consistenly and global asymptotic stable. The Figure 1(a), depicts a unique interior equilibrium point $E^{\text{*}}\left(1.2309,0.6785\right)$ for the model (2.3) with $k=1$ ; Figure 1(b) depicts a unique interior equilibrium point $E^{\text{*}}\left(1.0316,0.58485\right)$ for the model (2.3) with $k=3$ ; Figure 1(c) depicts a unique interior equilibrium point $E^{\text{*}}\left(0.7501,0.4525\right)$ for the model (2.3) with $k=10$ . Figure 1(d) depicts a unique interior equilibrium point $E^{\text{*}}\left(0.3333,0.2566\right)$ for the model (2.3) with $k=100$ . As can be seen from Figure 1, as the fear level K increases, the population densities of both prey and predator decrease. This indicates that the fear level is detrimental to the increase in population densities of both prey and predator. Especially, in Figure 1(d), even with a high fear level $k$ , the prey population still exists. In other words, having a fear effect with saturated fear cost in the system (2.3), redator and prey populations can coexist with minimal fear cost.

Figure 1. Population densities of prey and predator under varying levels of fear $k$ . The values of the other parameters are provided in (5.1).

Example 5.2. Considering the following parameters, let’s explore the impact of prey refuges on population density in ecology.

$r=0.5,\eta =0.2,\alpha =1,d=0.1,\beta =0.5,k=2,s=0.2,a=0.5,c=0.1,$ (5.2)

and vary the parameter $m$ .

Figure 2. Population densities of prey and predator under varying prey refuge m. The values of the other parameters are provided in (5.2).

As can be seen from Figure 2(a), here exists a positive correlation between prey population density and prey refuge, suggesting that an increase in prey refuge will elevate the prey population density. In Figure 2(b), when $0<m<0.843$ , the predator population density rises with the augmentation of prey refuge; conversely, when $0.843<m<1$ , the predator population density diminishes as the prey refuge expands. Notably, at the point where $m=0.843$ , the predator population density attains its maximum. This phenomenon can be ecologically interpreted as follows: When prey refuges are scarce, the number of sheltered prey is also limited, and at this juncture, the prey population continues to grow, enabling predators to capture sufficient prey to sustain their existence. However, as the number of prey refuges increases, more prey seek shelter within, making it challenging for predators to hunt enough prey to maintain a high population density, thereby leading to a decline in the predator population density.

Example 5.3. To visualize the impact of the Predator-taxis sensitivity $\alpha$ on the dynamic behavior of system (2.3), the following parameters are selected:

$r=0.2,\eta =0.25,m=0.3,d=0.2,\beta =0.5,k=5,s=0.2,a=0.8,c=0.2,$ (5.3)

and vary the parameter $\alpha$ .

Under the condition of parameter (5.3), in Figure 3(a), $\alpha =0$ , system (2.3) has a unique globally asymptotically stable equilibrium $E^{\text{*}}\left(0.1667,0.3333\right)$ . In Figure 3(b), $\alpha =1$ , the dynamics is similar to that in Figure 3(a), and the greates difference between these two figures is that the value of $E^{\text{*}}$ decreases from (0.1667, 0.3333) to (0.0663, 0.2451). With further increase of parameter $\alpha$ , the value of $E^{\text{*}}$ decreases from (0.0663,0.2451) in Figure 3(b) to (0.0191, 0.2118) in Figure 3(c). In Figure 3(d), $\alpha =9$ , system (2.3) has no unique globally asymptotically stable equilibrium $E^{\text{*}}$ , the prey-free equilibrium $E_2\left(0,0.2\right)$ of system (2.3) is globally asymptotically stable. This indicates that an increase in Predator-taxis sensitivity $\alpha$ is not conducive to the population density of prey and predators. When the Predator-taxis sensitivity is higher than a certain threshold, the unique coexistence equilibrium point of system (2.3) disappears, and system (2.3) tends to stabilize at the prey-free equilibrium.

Figure 3. Impact of the Predator-taxis sensitivity $\alpha$ on the dynamic behavior of system (2.3). The values of the other parameters are provided in (5.3).

6. Conclusions

In this paper, the dynamical behaviors of a modified Leslie-Gower predator-prey model incorporating fear effect and prey refuge are investigated. We consider that the predator population has substitute food source, enabling the predators to survive even if the prey population becomes extinct. Taking into account the limited energy of prey populations for both foraging and seeking refuge, we hypothesize that an excessive investment in refuge-seeking behaviors will necessarily detract from time spent foraging. We bridge these two independent factors through the Predator-taxis sensitivity $\alpha$ , yielding System, which holds profound biological significance. When examining the impact of fear effects on predator populations, we contend that predator fear does not necessarily lead to prey extinction. Rather, at high levels of fear, prey populations can persist through various anti-predator behaviors. This assertion is corroborated by the validation presented in Figure 1(d).

Next, we demonstrate the non-negativity, boundedness, and persistence of the model under the condition of non-negative initial values. The stability of equilibrium points is investigated by computing the real parts of the eigenvalues of the Jacobian matrix evaluated at these points. Our findings reveal that the trivial equilibrium and the predator-free equilibrium are always saddle points, whereas the prey-free equilibrium is a stable node under certain conditions, otherwise, it is also a saddle point. In the presence of a coexistence equilibrium, we prove that it is always locally asymptotically stable. Under such circumstances, the prey-free equilibrium becomes unstable, and according to the Bendixson-Dulac criterion, the system (2.3) exhibits no closed orbits, thus confirming the global asymptotic stability of the coexistence equilibrium. In the scenario where a coexistence equilibrium exists, we employ the implicit function theorem to derive the derivatives of prey and predator population densities with respect to the fear level $k$ and the prey refuge coefficient $m$ . The computational results indicate that the fear effect is detrimental to the increase in both prey and predator population densities. Conversely, the prey refuge positively contributes to the enhancement of prey population density. When the prey refuge is relatively small, it acts as a facilitator for predator population density. However, as the prey refuge exceeds a certain threshold, it exerts an inhibitory effect on the predator population density.

Finally, we validate the aforementioned conclusions through numerical simulations. Under the condition of parameter (5.1), in Figure 1, it illustrates the impact of fear effect on prey and predator population densities. As the fear level increases, both prey and predator population densities decrease. Notably, in Figure 1(d), when the fear level is exceedingly high, the prey population persists, which aligns with observations in the real world. Under the condition of parameter (5.2), in Figure 2, it demonstrates the influence of prey refuge on prey and predator population densities. Specifically, prey refuge fosters the growth of prey population density. Furthermore, a low level of prey refuge stimulates predator population density, whereas a high level of prey refuge inhibits the growth of predator population density. Under the condition of parameter (5.3), we numerically validate the effect of the Predator-taxis sensitivity $\alpha$ on the dynamical behavior of system (2.3). As evident in Figure 3, as $\alpha$ progressively increases, the population densities of both prey and predator diminish. Notably, upon exceeding a specific threshold value of $\alpha$ , the unique positive equilibrium point of system (2.3) vanishes, and the system converges to stability at the predator-free equilibrium.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Lotka, A.J. (1925) Elements of Physical Biology. Williams & Wilkins.
[2]	Volterra, V. (1928) Variations and Fluctuations of the Number of Individuals in Animal Species Living Together. ICES Journal of Marine Science, 3, 3-51. https://doi.org/10.1093/icesjms/3.1.3
[3]	Leslie, P.H. (1948) Some Further Notes on the Use of Matrices in Population Mathematics. Biometrika, 35, 213-245. https://doi.org/10.1093/biomet/35.3-4.213
[4]	Leslie, P.H. and Gower, J.C. (1960) The Properties of a Stochastic Model for the Predator-Prey Type of Interaction between Two Species. Biometrika, 47, 219-234. https://doi.org/10.1093/biomet/47.3-4.219
[5]	Aziz-Alaoui, M.A. and Daher Okiye, M. (2003) Boundedness and Global Stability for a Predator-Prey Model with Modified Leslie-Gower and Holling-Type II Schemes. Applied Mathematics Letters, 16, 1069-1075. https://doi.org/10.1016/s0893-9659(03)90096-6
[6]	Shi, T. and Wen, Z. (2024) Canard Cycles and Homoclinic Orbit of a Leslie-Gower Predator-Prey Model with Allee Effect and Holling Type II Functional Response. Qualitative Theory of Dynamical Systems, 23, Article No. 197. https://doi.org/10.1007/s12346-024-01059-z
[7]	Wang, K., Xu, X. and Liu, M. (2024) Global Hopf Bifurcation of a Diffusive Modified Leslie-Gower Predator-Prey Model with Delay and Michaelis-Menten Type Prey Harvesting. Qualitative Theory of Dynamical Systems, 23, Article No. 81. https://doi.org/10.1007/s12346-023-00939-0
[8]	Wang, D. and Ma, Y. (2024) Bifurcation and Turing Pattern Formation in a Diffusion Modified Leslie-Gower Predator-Prey Model with Crowley-Martin Functional Response. Journal of Applied Mathematics and Physics, 12, 2190-2211. https://doi.org/10.4236/jamp.2024.126133
[9]	Lima, S.L. (1998) Nonlethal Effects in the Ecology of Predator-Prey Interactions. BioScience, 48, 25-34. https://doi.org/10.2307/1313225
[10]	Maerz, J.C., Panebianco, N.L. and Madison, D.M. (2001) Effects of Predator Chemical Cues and Behavioral Biorhythms on Foraging, Activity of Terrestrial Salamanders. Journal of Chemical Ecology, 27, 1333-1344. https://doi.org/10.1023/a:1010309108210
[11]	Sheriff, M.J., Krebs, C.J. and Boonstra, R. (2009) The Sensitive Hare: Sublethal Effects of Predator Stress on Reproduction in Snowshoe Hares. Journal of Animal Ecology, 78, 1249-1258. https://doi.org/10.1111/j.1365-2656.2009.01552.x
[12]	Zanette, L.Y., White, A.F., Allen, M.C. and Clinchy, M. (2011) Perceived Predation Risk Reduces the Number of Offspring Songbirds Produce per Year. Science, 334, 1398-1401. https://doi.org/10.1126/science.1210908
[13]	Creel, S. and Christianson, D. (2008) Relationships between Direct Predation and Risk Effects. Trends in Ecology & Evolution, 23, 194-201. https://doi.org/10.1016/j.tree.2007.12.004
[14]	Cresswell, W. (2010) Predation in Bird Populations. Journal of Ornithology, 152, 251-263. https://doi.org/10.1007/s10336-010-0638-1
[15]	Wang, X., Zanette, L. and Zou, X. (2016) Modelling the Fear Effect in Predator-Prey Interactions. Journal of Mathematical Biology, 73, 1179-1204. https://doi.org/10.1007/s00285-016-0989-1
[16]	Das, B.K., Sahoo, D. and Samanta, G.P. (2022) Impact of Fear in a Delay-Induced Predator-Prey System with Intraspecific Competition within Predator Species. Mathematics and Computers in Simulation, 191, 134-156. https://doi.org/10.1016/j.matcom.2021.08.005
[17]	Ishaque, W., Din, Q., Khan, K.A. and Mabela, R.M. (2023) Dynamics of Predator-Prey Model Based on Fear Effect with Bifurcation Analysis and Chaos Control. Qualitative Theory of Dynamical Systems, 23, Article No. 26. https://doi.org/10.1007/s12346-023-00878-w
[18]	Sarkar, K. and Khajanchi, S. (2020) Impact of Fear Effect on the Growth of Prey in a Predator-Prey Interaction Model. Ecological Complexity, 42, Article ID: 100826. https://doi.org/10.1016/j.ecocom.2020.100826
[19]	Dong, Y., Wu, D., Shen, C. and Ye, L. (2021) Influence of Fear Effect and Predator-Taxis Sensitivity on Dynamical Behavior of a Predator-Prey Model. Zeitschrift für angewandte Mathematik und Physik, 73, Article No. 25. https://doi.org/10.1007/s00033-021-01659-8
[20]	Hamdallah, S.A.A. and Arafa, A.A. (2023) Stability Analysis of Filippov Prey–predator Model with Fear Effect and Prey Refuge. Journal of Applied Mathematics and Computing, 70, 73-102. https://doi.org/10.1007/s12190-023-01934-z
[21]	Sih, A. (1987) Prey Refuges and Predator-Prey Stability. Theoretical Population Biology, 31, 1-12. https://doi.org/10.1016/0040-5809(87)90019-0
[22]	Huang, Y., Chen, F. and Zhong, L. (2006) Stability Analysis of a Prey-Predator Model with Holling Type III Response Function Incorporating a Prey Refuge. Applied Mathematics and Computation, 182, 672-683. https://doi.org/10.1016/j.amc.2006.04.030
[23]	Ma, Z., Li, Z. and Wang, S. (2008) Effect of Prey Refuges on a Predator-Prey System. Journal of Lanzhou University, 44, 103-106.
[24]	Mandal, S. and Tiwari, P.K. (2024) Schooling Behavior in a Generalist Predator-Prey System: Exploring Fear, Refuge and Cooperative Strategies in a Stochastic Environment. The European Physical Journal Plus, 139, Article No. 29. https://doi.org/10.1140/epjp/s13360-023-04787-4
[25]	Ghosh, J., Sahoo, B. and Poria, S. (2017) Prey-Predator Dynamics with Prey Refuge Providing Additional Food to Predator. Chaos, Solitons & Fractals, 96, 110-119. https://doi.org/10.1016/j.chaos.2017.01.010
[26]	Li, W., Zhang, C. and Wang, M. (2024) Analysis of the Dynamical Properties of Discrete Predator-Prey Systems with Fear Effects and Refuges. Discrete Dynamics in Nature and Society, 2024, Article ID: 9185585. https://doi.org/10.1155/2024/9185585
[27]	Sasmal, S.K. and Takeuchi, Y. (2020) Dynamics of a Predator-Prey System with Fear and Group Defense. Journal of Mathematical Analysis and Applications, 481, Article ID: 123471. https://doi.org/10.1016/j.jmaa.2019.123471
[28]	Pielou, E.C. (1969) An Introduction to Mathematical Ecology. Wiley.
[29]	Chen, F. (2005) On a Nonlinear Nonautonomous Predator-Prey Model with Diffusion and Distributed Delay. Journal of Computational and Applied Mathematics, 180, 33-49. https://doi.org/10.1016/j.cam.2004.10.001
[30]	Murray, J.D. (2002) Mathematical Biology I. Springer.
[31]	Krantz, S.G. and Parks, H.R. (2003) The Implicit Function Theorem: History, Theory and Applications. Birkhauser.