Bifurcation Analysis in a Harvested Modiﬁed Leslie–Gower Model Incorporated with the Fear Factor and Prey Refuge

: Fear and prey refuges are two signiﬁcant topics in the ecological community because they are closely associated with the connectivity of natural resources. The effect of fear on prey populations and prey refuges (proportional to both the prey and predator) is investigated in the nonlinear-type predator-harvested Leslie–Gower model. This type of prey refuge is much more sensible and realistic than the constant prey refuge model. Because there is less research on the dynamics of this type of prey refuge, the current study has been considered to strengthen the existing literature. The number and stability properties of all positive equilibria are examined. Since the calculations for the determinant and trace of the Jacobian matrix are quite complicated at these equilibria, the stability of certain positive equilibria is evaluated using a numerical simulation process. Sotomayor’s theorem is used to derive a precise mathematical conﬁrmation of the appearance of saddle-node bifurcation and transcritical bifurcation. Furthermore, numerical simulations are provided to visually demonstrate the dynamics of the system and the stability of the limit cycle is discussed with the help of the ﬁrst Lyapunov number. We perform some sensitivity investigations on our model solutions in relation to three key model parameters: the fear impact, prey refuges, and harvesting. Our ﬁndings could facilitate some biological understanding of the interactions between predators and prey.


Introduction
Several researchers have explored the dynamical intricacy of interacting prey-predator models in depth to comprehend the species' long-term behavior.A wide range of models have been developed to investigate the dynamic relationship between prey and their specialist predators.These models rely on the mathematical formulation by Lotka and Volterra.It should be noted, however, that certain models are classified as Gause-type models [1][2][3][4].The Leslie-Gower form relies on the premise that a predator's population reduction and the per-capita availability of its prey are inversely correlated.As a matter of fact, Leslie [5] established a predator-prey model in which the carrying capacity of the predator's habitat is inversely correlated with the amount of prey.Leslie and Gower [6] and Pielou [7] have also examined an intriguing approach to predator dynamics.Similar to this, the researchers investigated the Leslie-Gower model, taking into account the effects of other ecological aspects, including harvesting, the Allee effect, prey refuges, and cooperative hunting.The functional response, or the interaction term, is another crucial element of population dynamics.Several ecologists and biologists have modified and analyzed various functional responses in the literature, including the Holling, ratio-dependent, Beddington-DeAngelis, and Crowley-Martin responses.For instance, the predator-prey model that states the predator eats its preferred food (prey) in a ratio-dependent manner was taken into consideration by the authors of [8].They carefully examined how the Allee factor and the fear factor affect the size of the prey population.Also, they showed that the model under consideration experiences a number of bifurcations, which include saddle-node bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation.Local and global stability analysis for the Leslie-Gower model and its modified model have been carried out by many researchers [9][10][11].
Despite the fact that some ecologists and biologists agree that prey-predator interactions are inadequately defined solely by regular predation and that a certain amount of fear must be assumed, few mathematical models have been constructed to quantify the possibility that fear factors influence the size of the prey population.This is partly due to the absence of specific experimental data illustrating how populations of terrestrial vertebrates can be impacted by fear [12,13].Zanette et al. [14] just finished manipulating song sparrows throughout the course of a mating season to see if the assumed predation risk may impact reproduction despite the lack of direct killing.Fear may also have an impact on a juvenile prey's physiological state, which might reduce its odds of surviving to adulthood.Birds, for instance, have predator defenses that they activate when they hear a predator, and when they are breeding they will leave their nests as soon as a threat is detected.Such anti-predator behavior can have long-term detrimental impacts on reproduction and population growth, even though it may enhance the probability of survival.Recent mathematical research has investigated the various changes in dynamical behavior that the fear effect might bring about in prey-predator models [15].The modified Leslie-Gower model was employed by the authors in [16] to explore how fear affects the population of prey.They looked into the existence of different bifurcation behaviors and showed that the system experiences a series of dynamic behavior switches as the cost of fear rises, which eventually cause the prey population to go extinct while the predators survive because alternative prey is abundant.In [17], the authors studied the modified Leslie-Gower model with the impact of the fear effect and nonlinear harvesting in both prey and predators.Up until certain thresholds, it is seen that the fear rate stabilizes the system; after that, it causes prey extinction.
The prey refuge helps to safeguard prey to some extent and improves species coexistence by lowering the risk of extinction caused by predation.The refuge lowers preypredator interactions and increases prey survival from extinction, according to a wide range of observational and empirical evidence [18,19].Prey refuges have also been the subject of several theoretical investigations into their impacts, and it has generally been shown that they promote persistence by stabilizing the prey-predator relationship [20].According to some empirical research, refuges can save prey from going extinct by stabilizing the community equilibrium and lowering the predator-prey interactions' tendency to oscillate.Most recently, Al-Salti et al. [15] considered and studied the prey-predator model in regards to a refuge factor with a variable carrying capacity.Further, from the perspective of fulfilling human requirements, fishing, forestry, and animal management frequently practice population harvesting and the exploitation of biological resources.There is considerable interest in using bioeconomic modeling to shed light on the scientific administration of natural resources like fisheries and forests with regard to protection for the future benefit of humankind.Hence, we consider prey refuges in the harvested predator-prey model in this study.
Inspired by the existing literature, this work uses the nonlinear, harvested modified Leslie-Gower model with Holling type II interaction.Further, the study examines the results of including a fear of predation and refuge to balance out predation.To the extent of our knowledge, there is less research concerning the modified Leslie-Gower model with predator harvesting, a refuge factor proportional to both species, and a fear factor on the prey population.The conditions for various local bifurcation occurrences are described.The suggested model is illustrated by extensive numerical computations in terms of phase portraits and bifurcation diagrams.In this article, we studied details of methods followed by [21,22].
In this study, we will address the following three important aims: i To explore how the prey population fear and prey refuge impact the dynamical behavior in the proposed system.The prey refuge considered here depends on both the prey and predator and this paper shows how it impacts long-term survival.ii To investigate the necessary parametric conditions for the existence of equilibrium points, local and global stability, and bifurcation around the coexisting equilibrium point.iii To determine whether the system solution is close to a coexisting equilibrium point or exhibits periodic oscillations by examining the initial sizes of predators and prey.
The paper is structured as follows.In Section 2, we give a deep formulation of a mathematical model to consider in this work.Basic properties such as positive invariance, boundedness, and persistence are provided in Section 3. The conditions for the existence and local and global stability of all possible equilibria are given in Section 4. The conditions for the occurrence of various bifurcation behavior in the proposed model are in Section 5.

Mathematical Formulation
First, the Leslie-Gower models are formulated with coupled, nonlinear ordinary differential equations, which illustrate the interaction between prey and their specialist predators in the form below: It is considered that prey grow logistically and have a linear increase in their intake rate with food density, where a denotes the attack rate.Also, dy dt = r 2 y(1 − y bx ) denotes the growth of the predator population is of the logistic form, but the traditional K 2 , which assesses the carrying capacity determined by the available resources in the ecosystem, is K 2 = bx, proportional to the number of prey (b is the ratio of prey to predators).The Lesie-Gower term in this equation is denoted by y/bx.It assesses the decline in the predator population brought on by the scarcity (per capita y/x) of the prey.If there is an extreme shortage, y can search for other species, but this will restrict its growth because its preferred food (x) is not widely accessible.A positive constant might be added to the denominator to solve this problem.Then, the model becomes the modified Leslie-Gower model [23].Some recent studies on Leslie-Gower models are cited in [16].In this work, we confine ourselves to the modified Leslie-Gower model, which is of the form The modified Leslie-Gower model is obtained by adding a positive constant n to the denominator in the second equation of model (1).

Fear Effect on Prey Population
Model (2) with the assumption that the presence of predators causes fear in the prey population is provided by Here, we multiplied the function f (k, y) = 1 1+ky in the birth rate of the prey population due to fear induced by predators, where k is the amount of anti-predator defense due to fear [16,17].Hence, the function f (k, y) = 1 1+ky meets a biological meaning and needs to satisfy the following conditions: The above biological assumptions are given in Appendix A. For other important results on considering fear levels in prey-predator dynamics, see [24][25][26][27].

Harvesting in Predator Population
In the disciplines of fisheries, forestry, and wildlife management, using biological resources and harvesting species are common practices.It is important to note that harvesting is continuing for a long time before the extinction of the population.Harvesting can be implemented in three ways: (a) a constant-yield H(y) = constant, (b) a constant-effort H(y) = qEy, or (c) a Michaelis-Menten type H(y) = qEy cE+ly .The modified Leslie-Gower model with a nonlinear harvesting rate in predators is demonstrated below: A version of the model (4) without the fear effect has been studied in [28].In constant effort harvesting, y is finite and fixed or as the y → ∞ if E is finite and fixed.Hence, the significance of choosing nonlinear harvesting is that the unrealistic features have been removed such that qEy cE+ly → qy c as E → ∞ and qEy cE+ly → qE l as y → ∞.Nonlinear harvesting is more realistic from an economic and biological point of view than other types of harvesting.

Prey Refuge
Next, taking into account prey refuges proportional to both populations brings our model system closer to reality since, in certain natural systems, prey refuges may be impacted by the size of both predators and prey.In light of this, the current work attempts to examine how refuges and harvesting affect a Holling type II prey-predator model [29].However, research into the Leslie-Gower model with group defense is infrequent.Sokol and Howell [30] proposed a modified Leslie-Gower predator-prey system with group defense.Prey refuges were included in the model investigated in [30], and the spatial component revealed that species distributions are highly susceptible to group defense compared to prey refuges [31].
Model (4) includes a prey refuge that is proportional to both populations, i.e., the refuge size is µxy from the predator, with µ ∈ [0, 1].The incorporation of a prey refuge allows (1 − µy)x of the prey to be accessible for the predator to hunt [10].To guarantee that the allowable range of the refuge is 0 ≤ (1 − µy) ≤ 1 for a realistic biosystem [32], we only permit those tiny values of µ such that (1 − µy) ≥ 0, i.e., y ≤ 1 µ .Next, model ( 4) is provided by where x and y stand for the prey and generalist predator densities, respectively, at time t; r 0 and d 1 are the growth and death rate of the prey population; K 1 is the prey's carrying capacity; the term y/(1 + mx) is the Holling type II function [2,33,34] where m > 0 represents a reduction in the predation rate at high predator densities due to mutual interference among the predators while searching for food or m is the product of the feeding rate and processing time, i.e., processing and searching for food are mutually exclusive events; and r 2 denotes the growth rate of the predator population size.The biological meaning of all parameters are presented in Table 1., a 2 = r 2 b , α = qE l and β = cE l , then the model with fewer parameters takes the following form: where x(0) = x 0 and y(0) = y 0 are the initial conditions for model (6).Hereafter, we consider model (6) for analysis in the below sections.

Some Preliminaries
In this section, we provide some preliminaries, such as the existence, positivity, persistence, and boundedness of the solution for model (6).The persistence in a predator-prey model plays a significant role in mathematical ecology that guarantees the long-term exis-tence of all species.Here, we have demonstrated the persistence criterion using the average Lyapunov function [35].

Existence and Positive Invariance
) t , and model ( 6) is given as where Fi ∈ C ∞ (R) and i = 1, 2. Since F is a smooth function of the variables (x, y) in the positive quadrant , where F is Lipschitzian function of two variables, the local existence and uniqueness of the solution set the hold on the Ω.

Persistence
Let us consider where p 1 and p 2 are positive constants.Now, define the function Φ as: Hence, the solution of model ( 6) is permanent if are satisfied.
Proof.The given initial condition for model (6) ensures that a unique solution (x(t), y(t)) exists and is defined on the interval [0, ξ), where 0 < ξ ≤ +∞.This is due to the fact that the right-hand side of model ( 6) is continuously differentiable and locally Lipschitz in the domain R 2 + , ensuring the existence and uniqueness of the solution.From model (6), we have It can be concluded that model ( 6) is positively invariant for all t ≥ 0.
Theorem 2. All solutions starting in R 2 + are uniformly bounded.
Proof.From model (6), and it is clear that x(t) ≤ r 0 −d 1 a 1 := M 1 as t → ∞, since it is considered that the intrinsic growth rate of the prey size should be positive, i.e., r 0 − d 1 > 0.
Then, using the maximum M 1 in the second equation of model ( 6), we have where µ 1 = 1 − µy > 0.Then, from the above inequality, we have which ensures the boundedness of the solutions.

Equilibria and Their Stability
This section investigates the equilibria points' existence in the proposed model and provides a qualitative assessment of their stability.To obtain the equilibria of model ( 6), we need to solve the prey and predator nullcline equation, which is We can see that model (6) possesses the following equilibria: 1.The trivial equilibrium E 0 (0, 0) always exists.
4. Now, we are interested in the coexistence equilibrium point E * (x * , y * ), where x * is calculated from the second equation of (10), which is where , and y * is calculated from the following equation: where Here, the coefficients B i , i = 1, 2, . . ., 8 in equation ( 13) are influenced by the system parameters.The analytical expression for equilibria to the aforementioned model is extremely challenging.By solving numerically, one can derive the coexistence equilibrium E * .Then, we have the following lemma on the existence of various equilibria.
Remark 1.For model ( 6), E 0 and Ê always exist.Ē exists if r 2 > α β holds.From Descartes's rule of sign changes, it is clear that B 1 > 0 and B 8 < 0 if a 1 A 1 > A 4 (r 0 − d 1 ).Then, (13) has at least one positive root.Furthermore, the number of sign changes in (13) can determine the number of coexistence equilibria of model (6), and it is also necessary that x * > 0. Equation ( 13) has exactly one positive root if there is only one sign change occurring in the coefficients B i .Hereafter, we consider that if model (6) has two coexistence equilibria they are named E 1 and E 2 .If it has only one coexistence equilibrium, it is named E 1 .
Furthermore, the nullcline plot illustrates how the number of equilibria varies by varying the fear parameter with the other parameters held constant.For the choice of fixed parameter values k, for instance, Figure 1 illustrates the existence of two, one, and no interior equilibrium points for model (6).For the choice of fixed parameter values r 0 = 1.01, d = 0.01, a 1 = 1, e = 0.85635, µ = 0.015, m = 0.015, r 2 = 2.60, a 2 = 0.4142, n = 0.01, α = 1.475193, and β = 0.2 and varying the fear parameter k, model ( 6) has the two interior equilibrium points E 1 = (0.28693, 0.79719) and E 2 = (0.60581, 0.44076) for k = 0.05, unique equilibrium E 1 = (0.37126, 0.52976) for k = 0.4093, and no equilibria for k = 0.5.Similar changes in the number of equilibria can occur by varying both the prey refuge µ and harvesting parameter α.

Local Stability
We now examine the local stability properties around the equilibria for model (6).For any given equilibrium (x, y), the Jacobian matrix is where The local stability properties of the equilibria E 0 , Ê, and Ē are stated in the following theorem.
The eigenvalues of the The following theorem provides the conditions for the local asymptotic stability of an arbitrary interior equilibrium point E * .Theorem 4. The coexistence equilibrium E * is locally asymptotically stable if and only if Tr(J E * ) < 0 and Det(J E * ) > 0.
Proof.The Jacobian matrix at E * is given by where The characteristic polynomial of J E * is given by where Tr(J E * ) = a 11 + a 22 and Det(J E * ) = a 11 a 22 − a 12 a 21 .
Furthermore, the dynamical behaviors are summarized in Table 2.
Table 2. Existence and stability conditions of equilibria for model (6).

Equilibrium and Feasibility Stability Coordinate
Condition Status E * (x * , y * ) at least one sign change Tr(J E * ) < 0 and in coefficients B i in ( 13) Det(J E * ) > 0 and x * > 0 Then, we conclude that E 0 and Ê always exist, i.e., the birth rate of the prey population r 0 is always greater than d 1 .Additionally, if Ê exists and is stable it violates the condition for the existence of Ē.Moreover, the relationship between the effects of other parameters has been discussed in the numerical simulations section.

Global Stability
Local stability only guarantees the behavior of the system in the small neighbourhood of an equilibrium point.It does not provide information about the long-term behavior of the system or its stability properties over the entire state space.In contrast, global stability refers to the property of a dynamical system where all trajectories, regardless of their initial conditions, converge to a stable equilibrium point.In this subsection, we concentrate on the global coexistence property of the arbitrary coexistence equilibrium point E * .The global stability properties of model ( 6) are attained by considering a suitable Lyapunov function.The Lyapunov function used in this article has been widely considered in [32,36,37].From Remark 1, if model (6) has exactly one coexistence equilibrium point E 1 (x 1 , y 1 ) then we have the following results for the globally asymptotically stable condition for E 1 (x 1 , y 1 ).
Proof.Let us take the suitable Lyapunov function V(x, y) : R 2 + → R as follows: where ) and V 2 (y) = y − y 1 − y 1 ln(y/y 1 ).The positive constant Ω is defined below.Both functions are well defined and continuous on R 2 + .V(x, y) is positive in R 2 + except at E 1 (x 1 , y 1 ) and V(x, y) = 0 at E 1 (x 1 , y 1 ).Moreover, ∂x < 0 for x < x 1 , and ∂V 2 (y) ∂y > 0 for y > y 1 , ∂V 2 (y) ∂y < 0 when y < y 1 .The time derivative of V 1 and V 2 at the solutions of model (6) after using E 1 is .
We define ).By substituting (18) and (19), we obtain Choosing Ω such that we obtain We can see that the coefficient of (y − y 1 ) 2 < 0 and it is necessary to have a negative in the coefficient of (x − x 1 ) 2 ; for this, we have Clearly, V is positive definite for (x, y) ∈ R 2 + \(x 1 , y 1 ).Since the quadratic form (23) in the previous equation is positively defined, every trajectory in the positive quadrant other than (x 1 , y 1 ) has dV dt < 0. As a result, if the condition is met then E 1 is globally asymptotically steady: dV dt | (x 1 ,y 1 ) = 0.Then, the Lyapunov function V constructed here follows the Lyapunov-Lasalle's invariance principle [38].

Bifurcation Analysis
The study of dynamical systems that undergo abrupt qualitative changes in behavior due to changes in their parameter values is the focus of the mathematical field of bifurcation.For further information on the foundations of local bifurcation analysis, see [21,22,39].In this section, the parametric conditions for the occurrence of various bifurcation behavior for model ( 6) is discussed here, which is saddle-node, transcritical, and Hopf bifurcation near E * by varying the fear, refuge, and harvesting parameter.

Hopf Bifurcation
In this subsection, we provide the condition for the existence of Hopf bifurcation around the arbitrary interior equilibrium E * , and the stability of the Hopf bifurcation is discussed with the help of the value of the first Lyapunov coefficient.Theorem 6.The interior equilibrium E * undergoes Hopf bifurcation if Tr(J E * ) = 0 and Det(J E * ) > 0.
Theorem 7. Assume that the system has an interior equilibrium E * and satisfies Theorem 6; then, E * changes its stability via Hopf bifurcation at some threshold k = k HB and satisfies d dk Tr[E * ]| k=k HB = 0.
Proof.We must confirm the Hopf bifurcation's transversality condition in order to guarantee the stability changes induced by nondegenerate Hopf bifurcation.Clearly, When the parametric limitation Tr[E * ] = 0 and the aforementioned transversality requirement are both met then E * loses stability via Hopf bifurcation.

Stability of Limit Cycle
To analyze the stability direction of the limit cycle, we proceed to compute the first Lyapunov number l 1 at the equilibrium point E * of model (6).
To begin with, we shift the equilibrium E of model ( 6) to (0, 0) by applying x = x − x * and y = ŷ − y * .This yields model ( 6) in a neighborhood of the origin as Since model ( 6) exhibits the Hopf bifurcation, we have ( â10 + b01 = 0) and ∆ = â10 b01 − â01 b10 > 0. The coefficients âij and bji are determinant by Then, the quantity l 1 says the stability of the limit cycle for model ( 6) is given by A supercritical Hopf bifurcation destabilizes the equilibrium E * when l 1 < 0, whereas a subcritical Hopf bifurcation occurs when l 1 > 0. Since the above expression for the first Lyapunov number l 1 is complex, we cannot determine the sign of l 1 from the above expression analytically.As a result, we have found it in the numerical part at the Hopf bifurcation point.

Nonexistence of Periodic Solution
Now, we will express model ( 6 • (H F) = ∂ ∂x According to Bedixson-Dulac's criterion for the nonexistence of periodic orbits [40], if • (H F) < 0 then the present system does not exhibit any periodic orbits.

Saddle-Node Bifurcation
Next, we provide the parametric condition for the occurrence of saddle-node bifurcation in the following.Theorem 8.If Det(J E * ) = 0 at a critical threshold µ = µ SN , then model ( 6) admits saddle-node bifurcation at E * .
Proof.Let us take a bifurcation parameter µ and apply Sotomayor's Theorem [39] to show model ( 6) admits a saddle-node bifurcation at the arbitrary equilibrium point E * .According to [39], since Det(J E * ) = λ 1 λ 2 = 0, then either λ 1 or λ 2 must be zero and the other less than zero.Also, Tr(J E * ) < 0. Let F = ( F1 , F2 ) T ; then, matrix J E * is written as where Let µ SD be the critical value, such that J E * has the eigenvalue zero, such that Det(J E * ) = 0 at µ SD .Also, let the matrices J E * and J E * T have the eigenvectors Θ = (θ 1 , θ 2 ) T and Φ = (φ 1 , φ 2 ) T for the zero eigenvalue, which gives Θ = (− â22 , â21 ) T and Φ = ( â22 , − â12 ) T .Furthermore, we obtain get Therefore, Ω 1 = 0 at µ = µ SD , and also where 3  , and also According to Sotomayor's Theorem, the system exhibits a saddle-node bifurcation around E * (x * , y * ) at µ = µ SD if Ω 1 = 0 and Ω 2 = 0. To confirm Ω 1 , Ω 2 = 0, we calculated it numerically.Thus, it can be inferred that the variation of the parameter µ across the critical threshold µ = µ SD results in a change in the number of interior equilibria of model ( 6) from zero to one to two.

Numerical Simulations
The dynamic nature of the proposed model that corresponds to changes in the fear, prey refuge, and harvesting parameters cannot be expressed explicitly or analyzed analytically.Therefore, numerical simulations are required to comprehend the model's dynamics, which were carried out with the help of XPPAUT [41] and PPLANE [42] software.Let us consider the model parameters as follows: In order to check the effect of fear, the prey refuge, and harvesting in the considered model (6), we consider the fixed parameter values in Table 3 and vary the k, µ, and α.The dynamics of the model are analyzed using phase portraits and one-or two-parameter bifurcation diagrams, as discussed below.For model (6) with the fixed parameters in Table 3, varying the fear parameter k exhibits complex behavior.For instance, the numbers for the existence of interior equilibria are shown using nullcline analysis (see Figure 1).To verify whether the species can coexist for a long period of time, it is crucial to examine the local stability of the equilib-ria.If k = 0.1, the system has four equilibrium points: E 0 (0, 0) is unstable, Ê = (1, 0) is unstable, E 1 (0.2884, 0.75692) is a spiral sink that satisfies Theorem 4 and whose eigenvalues are λ 1,2 = −0.066223± 0.81822i, and E 2 (0.58067, 0.44519) is a saddle point, which is shown in Figure 2a.If k = 0.15, the equilibrium point E 1 (0.29127, 0.72082) is a spiral sink surrounded by an unstable limit cycle, as shown in Figure 2b.On further increasing k HB = 0.1719, the equilibrium E 1 (0.293, 0.70599) is a center, which is shown in Figure 2c.Then, the eigenvalues are λ 1,2 = ±0.75151,which satisfies Theorem 6, and

Effect of Harvesting
Further, it is necessary to verify the effect of harvesting the predator population in the system dynamics.We fix all other parameters as in Table 3 and vary the harvesting parameter α.For the smaller value of α = 0.1, the system has the following equilibrium points: E 0 (0, 0), Ê(1, 0), E 1 (0.14202, 0.91344), and E 2 (1.1536, −0.16237) (see Figure 3a).In E 2 , the predator population size is negative and it is biologically meaningless.And on increasing α TC = 0.52, the equilibrium E 2 and Ê collide with each other and exchange their stability properties, which ensure the existence of transcritical bifurcation (see Figure 3b).We have exactly one coexistence equilibrium point E 1 (0.16494, 0.88891), which satisfies Theorem 5.Then, E 1 is globally asymptotically stable and also clearly shown in Figure 3b.

Effect of Prey Refuge
The phenomena of the existence of one, two, and no equilibria can also occur for the refuge parameter by fixing the parameters in Table 3 and varying the refuge parameter.If µ = 0.2, model ( 6) has four equilibria, E(0, 0), Ê(1, 0), E 1 (0.34004, 0.81828), and is a spiral sink, whose eigenvalues are λ 1,2 = −0.27115± 0.772464, and E 2 (1, 0) is a saddle point (see Figure 4a).If µ SN = 0.9374, the number of interior equilibria is reduced to one and it is a saddle point E 1 (0.71968, 0.56214) whose eigenvalues are λ 1 = −0.69553and λ 2 = 0.021212 (see Figure 4b), which ensures the occurrence of saddle bifurcation for the parameter µ.And the quantities Ω 1 = 0.196311 = 0 and Ω 2 = −1.51299= 0 ensure the existence of saddle-node bifurcation for model (6) by using the Sotomayor theorem.
For a clear view of such a critical transaction, the one-parameter bifurcation by varying k ∈ (0, 0.5), α ∈ (0, 1), and µ ∈ (0, 2) is shown in Figure 5a-c, and all other parameters are fixed as in Table 3.One of the prominent phenomena of the ecological system is extinction criteria and it also necessary to show the dynamics of model (6) with the combination of this parameter.We plotted the two-parameter bifurcation diagrams in Figure 6a-c.The extinction region, the stable and unstable region of the interior equilibria with the particular choice of parameter values, is clearly shown.Further, it is necessary to show how the model studied in this article differs from other earlier works.We compared some works from the existing literature in Table 4.The extinction of the population for the particular choice of k, µ, α, and other values is fixed as in Table 3 and is shown in Figure 7a-c

Reference
Prey Refuge Fear Effect Harvesting Functional Response [16] × × Holling IV [17] Constant Holling II [43] × × Holling IV This paper Both prey and predator Holling II The dynamics of model ( 2), the model with fear, in which the predator consumes the prey in the form of a Holling type II interaction, were examined by the authors in [16].Transcritical bifurcation, Hopf bifurcation, and Bogdanov-Takens bifurcation are all present in the under-considered model.The authors of [17] looked at nonlinear harvesting in both species, prey refuges, and fear in the prey population.They demonstrated how the model experiences transcritical and saddle-node bifurcations.Mukherjee [43] investigated that the modified Leslie-Gower model with Holling type IV interaction, used to study the fear impact on a predator-prey system, experiences similar sorts of bifurcations.In this study, we examined the modified Leslie-Gower model with nonlinear harvesting in the predator population, fear in the prey population, and a prey refuge proportional to both species.We showed that model (6) exhibits saddle-node bifurcation, transcritical bifurcation, Hopf bifurcation, the extinction of species, and cusp-type dynamics.Since the computation approach offers fascinating details about how ecological systems function, there are significant drawbacks because (a) nonlinearity in the system frequently results in highly sensible solutions to variations in system parameters and (b) we frequently lack knowledge of true parameter values.This suggests that the specific decisions we produced regarding the values of the parameters and the precise structure of functional relationships may have significantly influenced the results.Analytical methods are precise and broad but only work with very basic models.But even the simplest model's formal analysis can serve as an initial guide for particular ideas that are easily evaluated through numerical simulations.3 show that the trajectories move towards Ê, i.e., the predator population become extinct.

Conclusions
In this study, we studied the influence of fear in the prey population, the refuge effect proportional to both prey and predators, and nonlinear harvesting in the predator population.We have assumed that prey ensures the admissible range of a refuge 0 ≤ (1 − µy) ≤ 1 in our considered model.We also took into account a consideration of nonlinear harvesting in the predator population.To analyze the effect of this biological situation, we considered the modified Leslie-Gower model, where the prey consumes its favorite food in the form of a Holling type II interaction term.It is important in ecology that the predator can switch over to other food in the abundance of its favorite prey.The existence of all biologically feasible equilibria has been investigated in order to analyze the local dynamics of the considered model.The local stability analysis was carried out in order to show the long-term coexistence of populations, which was analyzed with the help of eigenvalues of the corresponding Jacobian matrices of the equilibria.Further, we took the conventional Lyapunov function to derive the analytical condition for the global stability of the coexisting equilibrium point.The complex dynamical behavior of the model was discussed in terms of bifurcation analysis and extinction of populations, since both play a crucial role in conserving the populations in the ecological system.We showed the occurrence of saddle bifurcation, transcritical bifurcation, and Hopf bifurcation by varying the fear parameter k, refuge parameter µ, and harvesting parameter α.The occurrence of transcritical and saddle-node bifurcation was verified with the help of normal form the-ory [39] and we also showed that the existence of subcritical Hopf bifurcation is confirmed with the help of calculating the first Lyapunov number, i.e., the unstable limit cycle.We also showed that the population becomes extinct by increasing this parameter at a certain critical threshold.
When studying interactions between prey and predators or other ecological systems, achieving sustainability is of utmost importance.This is often accomplished by establishing a stable limit cycle, which helps to maintain the system's long-term stability.But in our proposed model we found the coexisting equilibrium surrounded by an unstable limit cycle.So it is necessary to allow small fluctuations in the populations inside the limit cycle, otherwise it leads to the extinction of species.Enhancing the biological realism of the studied model is one of the challenges for future research.Biological systems are more complex and it is difficult to study their underlying mechanisms.It is necessary to study our model by considering various functional responses and possibly extend it to the food chain model.

Figure 1 .
k = 0.5 Nullcline plots: the orange color line depicts the prey nullcline, the pink color depicts the predator nullcline, and the red points denote the equilibrium points.

Figure 2 .
d dk Tr[E 1 ]| k=k HB = (− r 0 y 1 (k HB y 1 +1) 2 ) = −0.567061= 0. Also, the unstable periodic solution collides with the saddle point.The local amplification near the interior equilibrium point E 1 of Figure 2c is shown in Figure 2d.It is clear there exists a Hopf bifurcation on varying f .It is clear that the eigenvalues of the Jacobian matrix at E 1 (0.293, 0.70599) are λ 1,2 = ±0.70599and it also satisfies the transversality condition of Theorem 7. Further, we found the first Lyapunov coefficient is l 1 = 75.1529π> 0, and it is clear that the existing limit cycle is unstable.(a) represents E 1 is a stable spiral and E 2 is a saddle point.(b) represents the stable E 1 is a stable spiral surrounded by unstable periodic solutions.(c) represents the E 1 is a center.(d) represents the local amplification of figure (c) near E 1 .

Figure 3 .
(a) represents E 1 is stable, E 2 is unstable, and Ê is a saddle point.(b) represents E 1 is stable, and when Ê collides with E 2 it becomes a saddle point.(c) represents an unstable E 1 , and when an unstable limit cycle collides with E 2 it becomes a saddle point and Ê becomes a stable node.(d) represents E 1 is the center.

Figure 4 .Figure 5 .
Figure 4. (a) Locally asymptotically stable E 1 and saddle point E 2 .(b) E 1 and E 2 collide and become an unstable node.

Figure 7 .
Phase portraits of model(6) with particular values of k, µ, and α with other fixed parameter values in Table

Table 1 .
Biological meaning of the parameters.

Table 4 .
Earlier works on Leslie-Gower model.