Complex Dynamics of a Predator–Prey Interaction with Fear Effect in Deterministic and Fluctuating Environments

: Many ecological models have received much attention in the past few years. In particular, predator–prey interactions have been examined from many angles to capture and explain various environmental phenomena meaningfully. Although the consumption of prey directly by the predator is a well-known ecological phenomenon, theoretical biologists suggest that the impact of anti-predator behavior due to the fear of predators (felt by prey) can be even more crucial in shaping prey demography. In this article, we develop a predator–prey model that considers the effects of fear on prey reproduction and on environmental carrying capacity of prey species. We also include two delays: prey species birth delay inﬂuenced by fear of the predator and predator gestation delay. The global stability of each equilibrium point and its basic dynamical features have been investigated. Furthermore, the “paradox of enrichment” is shown to exist in our system. By analysing our system of nonlinear delay differential equations, we gain some insights into how fear and delays affect on population dynamics. To demonstrate our ﬁndings, we also perform some numerical computations and simulations. Finally, to evaluate the inﬂuence of a ﬂuctuating environment, we compare our proposed system to a stochastic model with Gaussian white noise terms.


Introduction
The dynamical nature of predator-prey interactions is genuinely influential in mathematical biology, particularly in ecosystems, where the prophecy of predator-prey interactions and the predation process play a crucial role in maintaining balance and biodiversity.Mathematical modeling is an effective method for studying the aforementioned biological processes [1,2].As a result, various models have been constructed and investigated.We formally categorize them as deterministic models, especially those containing systems governed by ordinary differential equations [2,3], partial differential equations [1,4,5], fractional differential equations [6], stochastic systems [7], systems including time lags [8], network models [9], etc.
Malthus, near the end of the eighteenth century, proposed a theory on population dynamics based on the idea that population growth is proportionate to the population density present in the ecosystem [10].In this case, the populace will increase (or decrease) exponentially, but the situation of exponentially increasing growth does not match the real-world conditions.Verhulst [11] proposed the logistic growth model which incorporates a constant per-capita population growth rate and intra-species competition proportional to the current population size.In this circumstance, the population increases or decreases as a function of its initial size, and in the long run, approaches the environment carrying capacity from, respectively, above or below that level.The traditional Lotka-Volterra model [12], which was formulated for the predator-prey interaction, describes predation as a linear, and hence unlimited, function of the prey population density.This original theoretical form was later replaced with saturated curves called functional responses (of predator to the prey density) and improved by including a logistic growth term for the prey species.To simulate realistic behavior among predator and prey species, a number of functional response forms were proposed.For the last few decades, the Lotka-Volterra model has been evolving, incorporating diverse assumptions, such as the Allee effect, fear effect, gestation delay, habitat complexity, additional food and switching feeding, to explore how the system behaves locally and globally [13][14][15].
As predation is easy to witness in nature, most studies on predator-prey systems only investigate the direct consumption of prey by the predator [16][17][18].However, it is recognised that the sheer presence of a predator could have psychological and physiological impacts on the prey that are more harmful than direct consumption to the population.In general, animals defend themselves against predators in the wild by altering habitat use, attentiveness, foraging behavior and physiological changes [19,20].Prey may move from a higher-to a lower-risk zone to reduce predation and thereby lose energy, specifically if the conditions in the low-risk zone is poor.Furthermore, scared prey may eat less, resulting in hunger and decreased reproduction [16,19].Prey species are always under psychological stress due to the threat of assault, and in some circumstances, they die solely through fear rather than direct consumption.Several studies have found that many prey species spend a lot of time for watching predators and less time gathering food, which leads to fewer eggs being produced [21][22][23].Zanette et al. [23] found that female song sparrows (Melospiza melodia) exposed to predatory noises produce 40% fewer fledglings than birds exposed to non-predatory sounds.Creel et al. [24] used numerous data to demonstrate that wolves implicitly affect elk reproductive physiology and population dynamics through the cost of fear.Wang et al. [18] first developed a model that included the fear effect and investigated whether large amounts of fear may stabilize the predator-prey model by eliminating oscillatory behavior.In 2019, Zhang et al. [25] introduced the dynamic behavior of a predator-prey interaction model, including both fear and prey refuge.In 2020, Wang and Zou [26] studied a model incorporating a fear effect in a system of ordinary differential equations as a cost; however, in this research article, they also considered an anti-predation strategy and digestion delay.In 2022, Das et al. [27] constructed and analysed a predator-prey system that introduces the cost of fear into the birth and death rates of the prey population and a gestation delay.
In biological processes, a time delay is a common element.Time delays may occur in population biology due to food digestion, maturation, newborn predators' gestation and other factors.We can get several extensive explanations regarding the usefulness of time delays in practical systems from the classical books [28,29].After consuming prey, the predator must wait a certain amount of time to reproduce, mediated by the gestation period.As a result, in model construction and biological elucidation, the time delay between prey capture and contribution to the predator's growth is crucial [30].Such time-delay elements add realism and complexity to the system.In general, delays exert destabilizing effects, which may occur as oscillations via Hopf bifurcations [31][32][33].Ruan [30] explored various forms of delays and the dynamics of the associated systems in a review of works on predator-prey models with discrete delay.
Motivated by the preceding discussion, we have developed a predator-prey model that considers the influences of fear on prey reproduction and the environment carrying capacity of prey species.It also takes into account two delays: a prey birth delay and a predator gestation delay.The effects of fear on the proposed dynamic system, and the impacts of the delay factors, are the key research topics.This work is organized as follows.The mathematical model is formulated in the next section.Section 3 is devoted to the system's well-posedness.The non-delayed system's equilibrium points and their local stability, along with the global stability, are examined in Section 4. Conditions for uniform persistence are determined in Section 5. We look at local bifurcations, such as transcritical and Hopf bifurcations, in Section 6. Section 7 deals with the stability analysis of the delayinduced model.To show our theoretical findings, we have performed several numerical simulations, which are reported in Section 8. Additionally, the mortality rates of both species are perturbed by Gaussian white noise terms to compare the deterministic system to the stochastic model.Finally, Section 9 deals with the discussion and conclusions of this work.

Model Formation
A typical predator-prey population model in which prey logistic growth is considered and prey is the only food source of predator (i.e., in the absence of prey, a predator will go extinct) can be written as: where x and y denote the prey and predator populations, respectively; r and k are the intrinsic growth rate and the carrying capacity of the prey, respectively, (in absence of predator); α is the intake rate; γ (0 < γ < 1) is the conversion rate of the predator; d 2 is the natural mortality rate of the predator; h is the prey handling time.Here, the interaction between prey and predator is considered as a Holling type-II functional response αx 1+αhx [34][35][36].Now, let r = b − d 1 , where b and d 1 denote the birth rate and natural death rate of the prey, respectively.Therefore, Equation (1) can be rewritten as: ( Recent field observations and the results of the evidence suggest that the mere presence of a predator could alter the natural behavior of prey, and that it could affect population size.Defensive behaviors such as prevention, alertness, alarm calls and gathering against the predator might reduce direct mortality from predation temporarily.Lifelong fitness of the prey species can be damaged by reducing the growth rate and fecundity because of reduced intake and mating opportunities.Thus, we multiply the reproduction term bx by f 1 (k 1 , y) and the carrying capacity k by f 2 (k 2 , y), where f 1 (k 1 , y) and f 2 (k 2 , y) satisfy the following conditions [18]: Here, the parameter k i , where i = 1, 2, describes the level of fear that drives the prey's antipredation behavior.In particular, let f 1 (k 1 , y) = 1 1+k 1 y and f 2 (k 2 , y) = 1 1+k 2 y .Thus, by incorporating the fear effects in model (2), we have obtained the following predator-prey system: where Throughout the analysis we will take b > d 1 , if not stated otherwise.The description and unit of the model parameters are presented in Table 1.
Table 1.List of model parameters with units and their biological significance [37].

Parameter
Description Unit b intrinsic birth rate of the prey Time Recently, Mondal and Samanta [38] explored a predator-prey interaction incorporating nonlinear prey refuge under the influence of the fear effect and additional food.The following describes how Model (3) differs from the system investigated by Mondal and Samanta [38].
1.In [38], the authors considered that fear of predators only suppresses the birth rate of the prey species.However, in this work, we make the assumption that fear of predators not only decreases the birth rate of prey, but also reduces the environmental carrying capacity of prey species.2. In [38], the authors presumed and analysed the dynamics of a predator-prey interplay by including a nonlinear prey refuge function φ(x, y) = m 1 xy a 0 +y , where a 0 (>0) is a halfsaturation constant for refuge function and m 1 (>0) is the coefficient of prey refuge.On the other hand, the hypothesis used in this work is that the prey population does not take refuge.3.In [38], the authors did not examine the implications of environmental stochasticity for the suggested predator-prey system.However, in this work, the influences of environmental fluctuations have been taken into account too.
Furthermore, Sarkar and Khajanchi [39] proposed and analysed a predator-prey system introducing the cost of fear as an indirect impact of predators on prey.The fundamental difference between Model (3) and the model studied by Sarkar and Khajanchi [39] is the form of a fear function.In the present work, it is assumed that when level of fear k i = 0, i = 1, 2, the fear function f i (k i , y) seems to have no effect on birth, but in [39], they considered the fear function in such a manner that when fear level is zero, there is always some minimal impact on the birth rate of prey; even if the predator population rises infinitely, prey species will be under minimum fear because of the physiological impact of the prey populations being habituated to fear from predator species.However, in this scenario, an increase in predator population would have a significant impact on birth rate of prey species, which might ultimately lead to the extinction of prey individuals.Since there are not enough experimental data to confirm this theory employed by [39], we assume fear functions as f i (k i , y) = 1 1+k i y , where i = 1, 2, as considered by Wang et al. [18].A predator-prey interaction model has been studied by Halder et al. [40] which takes into consideration the fear effect and multiple foraging techniques.They have included a harvesting term for both species and assumed a modified version of the Leslie-Gower functional response.They took the foraging rate of the predator as linear with Holling type-II or Holling type-III foraging.On the other hand, in this work, a logistic predator-prey model with a Holling type-II functional response has been explored in the context of the cost of fear due to predator.Studying the consequences of fear, as it is incorporated into the carrying capacity of prey, and the implications of breeding delay and predator gestation delay, are key interests of this work.

Positivity and Uniform Boundedness
Here, we check the existence, positivity and boundedness of the solution of System (3).Theorem 1.The solution of system (3) exists and is positive for all values of t ≥ 0.
Then, system (3) reduces to: Clearly, xΨ 1 (x, y) and yΨ 2 (x, y) are continuous functions of x and y, and also locally Lipschitzian in R 2 + .Therefore, a solution of (3) exists in [0, ζ), where 0 < ζ < ∞.Now, Equation (4) gives Hence, the theorem is proved.Theorem 2. All solutions of system (3) are uniformly bounded in the region B = {(x, y) : 0 < Proof.By taking the first equation of Model (3), we get It should be mentioned that a greater death rate is always harmful to all species and may cause extinction, which is not desirable from an ecological perspective; therefore, b > d 1 is assumed throughout this work.Now, let Ω(x(t), y(t)) = x(t) + y(t) γ .Then, Therefore, by Gronwall inequality, we have Therefore, all solutions of (3) enter in the region: This completes the proof.

Equilibria and Their Stability Analysis
In this section, we evaluate the proposed system's behavior analytically.

Equilibrium Points
System (3) has three feasible equilibrium points or steady states in R 2 + , namely: h and y * is the unique positive root of the equation: Here, Thus, we have These three equilibrium points are represented in Figure 1 using the parameter set mentioned in the figure caption.Proof.The Jacobian matrix J(0, 0) at E 0 (0, 0) is Therefore, the eigenvalues are (b − d 1 ) and −d 2 .Clearly, for the stability of steady state E 0 (0, 0), we must have b < d 1 and for instability b > d 1 .

Theorem 4. The axial steady state E
and unstable if Proof.The Jacobian matrix corresponding to β , 0 is given by: Thus, the eigenvalues of the Jacobian matrix and Therefore, Theorem 5.The coexistence steady state E Proof.The Jacobian matrix J(x * , y * ) at the coexistence steady state E * (x * , y * ) is given by: , where Therefore, all the eigenvalues of J(x * , y * ) will have negative real parts if j 11 < 0. Thus,
Proof.Consider a Lyapunov function: Then, Thus, dV dt < 0 if b < d 1 and dV dt = 0 at E 0 (0, 0).Therefore, using Lyapunov theorem, we can conclude that the steady state Proof.From the first equation of system (3), Now, second equation of ( 3) gives: Let us take Dulac function as: where m is to be specified later.Let us now calculate the divergence of the vector : where Therefore, it is sufficient to find a m so that condition (6) is satisfied.However, condition (6) holds if For convenience, let m + 1 = m.Then, (7) becomes: The existence of a m satisfying (8) means: Thus, there exists m such that D < 0 for (x, y) ∈ R 2 + .By the Bendixson-Dulac theorem, system (3) has no closed orbits in the first quadrant, provided 0 Hence, E * is globally asymptotically stable.

Uniform Persistence
Proof.Let us consider the following average Lyapunov function: where i > 0, for i = 1, 2. Here, (x, y) is a continuously differentiable non-negative function defined on R 2 + .Differentiating Equation ( 9) with respect to t, we have: β , 0 , we get: We choose 1 and 2 in such a way that Λ(0, Remark 3. Persistency implies instability of the boundary equilibria.

Transcritical Bifurcations
In bifurcation theory, a transcritical bifurcation is a particular kind of local bifurcation in which a stationary point exchanges its stability with another stationary point due to continuous incrementation of the bifurcation parameter.
Theorem 10.System (3) experiences a transcritical bifurcation around the trivial steady state E 0 (0, 0) with bifurcating parameter d 1 at d 1 be the critical value of d 1 for which exactly one eigenvalue of J(0, 0) is zero.Thus, d Therefore, by Sotomayor's theorem, a transcritical bifurcation occurs around E 0 at d Theorem 11.System (3) exhibits a transcritical bifurcation around the axial equilibrium point 2 be the critical value of d 2 for which J b−d 1 β , 0 has exactly one zero eigenvalue at 2 .Thus, d . Now, the eigenvectors of for the zero eigenvalue are U = (u 1 , u 2 ) T and W = (0, 1) T , respectively.
Here, u 2 = 1 and Therefore, by Sotomayor's theorem, a transcritical bifurcation is exhibited around E 1 by taking d 2 as the bifurcating parameter.
Proof.The proof is similar to the proof of Theorem 11.

Hopf Bifurcation
A Hopf bifurcation takes place whenever a periodic solution or limit cycle surrounding an equilibrium point appears or disappears as a parameter's value changes.We shall show that the proposed model exhibits a Hopf bifurcation as we change the value of k 1 .Now, the characteristic equation at the steady state E * (x * , y * ) is where We know that in the case of complex conjugate roots of Equation ( 10), the equilibrium point E * (x * , y * ) is a stable spiral if the real part of the roots of Equation ( 10) are negative and an unstable spiral for the positive real part.Therefore, a stability switch takes place when Equation (10) contains purely imaginary roots.
Theorem 13.System (3) exhibits a Hopf bifurcation around the coexistence equilibrium point E * (x * , y * ) when the bifurcation parameter k 1 crosses a threshold value k 1 is the critical value of k 1 at which Equation (10) gives purely imaginary roots, so T(k 1 , the zeros of Equation ( 10) take the following form: , where p(k 1 ) and q(k 1 ) are real valued functions of k 1 .By putting λ(k 1 ) = p(k 1 ) + iq(k 1 ) in Equation ( 10) and differentiating with respect to k 1 , we have By separating the real and imaginary parts, we get where .
By solving Equations ( 11) and ( 12) for ṗ(k 1 ), we get: Now, for 1 , we consider the following cases: Case-1 Hence, from Equation (13), we have Hence, from Equation ( 13), we have: According to the Hopf bifurcation theorem [3], system (3) switches its stability behavior through Hopf bifurcation, provided ) Thus, for the existence of Hopf bifurcation at 1 , we must have This completes the proof.
Theorem 14. System (3) exhibits a Hopf bifurcation around the coexistence equilibrium point E * (x * , y * ) when the bifurcation parameter k 2 crosses a threshold value k Proof.The proof is similar to the proof of Theorem 13.

Delayed Dynamical System
A time delay is present in many biological processes, both natural and artificial.Exploring the time lags renders the mathematical model more authentic than the nondelayed model.A delay differential equation also explains significantly better intricate dynamics compared to a simple differential equation.
The predator's fear diminishes the birth rate of prey biomass, resulting from a combination of psychological changes, foraging behavioral changes and other factors.The cost of fear cannot be instantaneous because prey requires time to assess the predation risk.As a result, we cannot employ the cost of fear diminishing the prey's birth rate right away.A certain time delay, τ 1 (say), is needed to perform the complete process.Again, in reality, the conversion of consumed prey into predator reproduction does not occur instantly; instead, there is a time delay for predator biomass gestation.Thus, we assume that a constant time delay called gestation delay (τ 2 ) that governs the reproduction of predator population following prey hunting.To acquire the rich dynamics of (3), we introduce both these delays τ 1 and τ 2 in system (3).Therefore, after incorporating delays, the modified form of system (3) is: Let C denote the Banach space of continuous functions ψ : with ξ = max {τ 1 , τ 2 }.
By linearizing system ( 14) about E * (x * , y * ), we have: where A, B, C are 2 × 2 matrices given by where The characteristic equation of the delay system is where As Equation ( 16) is transcendental, it would have an infinite number of complex roots.
To study the local stability behavior of E * , we should investigate the signs of real parts of the zeros of Equation ( 16), which is hard in the presence of both time delays.As a result, we first examine Equation ( 16) in the absence of delay, then in the presence of a single time delay.After that, the local asymptotic stability behavior conditions of equilibrium E * will be derived using the same analytic arguments mentioned in [42,43], in the presence of both the time delays.
Case I In absence of both the time delays, system ( 14) is reduced to system (3).We have already described the condition of the stability of E * in Theorem 5 in absence of τ 1 and τ 2 .
Case II Therefore, the characteristic equation, Equation ( 16), becomes where For τ 1 = 0, Theorem 5 gives the criteria for which all the zeros of Equation ( 17) are negative or have negative real parts, whereas for τ 1 > 0, it has an infinite number of roots.According to Rouche's theorem and continuity in τ 1 , the signs of the roots of Equation (17) will change if they cross the imaginary axis.Thus, by putting λ = iω (ω > 0) in (17) and separating the real and imaginary parts, we get By squaring and adding Equations ( 18) and ( 19), we get By substituting ω 2 = σ, we get the following equation in σ: If there is no positive real root in Equation ( 21), real ω does not exist.Thus, for any τ 1 > 0, E * is locally asymptotically stable.If Equation ( 21) has a positive real root, say σ * such that ω 10 = ± √ σ * .Thus, λ = ±ω 10 are two purely imaginary roots of (17).From ( 18) and ( 19), we get the values of τ 1 as: for j = 0, 1, 2, . . . .Thus, if Theorem 5 holds, system (14) will be LAS around the interior equilibrium E * for τ 1 = τ 2 = 0. Therefore, by Butler's lemma, E * is stable for , and E * is unstable for τ 1 ≥ τ 0 1 , provided the transversality condition is satisfied.We now verify the transversality condition Equation ( 17) with respect to τ 1 , we get: This gives dλ dτ 1 Therefore, the transversality condition is satisfied, and so Hopf bifurcation occurs at 2 ) = 0.The following theorem can now be stated.

Numerical Analysis
In this section, we show numerical simulations to reveal the dynamical behavior of the proposed system (3).First, we have selected the parameter set: {b = 5.5, 2, h = 1 and γ = 0.8}.Here, we see that the death rate, d 1 , of prey is greater than the birth rate, b, of the prey; therefore, as time goes by, the prey species will become extinct, and in the absence of prey, predator species will also become extinct.We have depicted the corresponding time series in Figure 2.Then, we have decreased the parameter d 1 from 5.6, while leaving all other parameters to be the same as in Figure 2, and we can see that at d    To study the stability nature of the predator-free (axial) equilibrium point, we took the parameter set: {b = 5.5, Here, b > d 1 and the stability condition We have drawn the corresponding time series in Figure 4.In the next stage, taking all other parameters to be the same as in Figure 4, we have decreased the value of b and noticed that at b [TC] = 0.1, E 1 exchanges its stability with that of E 0 (see Figure 5a).We have seen a similar type of behavior when we increase the parameter value of d 1 from zero onward up to d [TC] 1 = 5.5; the stability of E 1 exchanges with that of E 0 (see Figure 5b).In the case of parameter h, the stability of E 1 exchanges with that of the equilibrium point E * at h [TC] = 0.8112 (see Figure 6), as we diminish h, while taking all other parameter values to be the same as in Figure 4.However, the death rate d 2 of the predator and the conversion rate γ of the predator have effective impacts on system (3), as the existence of a positive interior equilibrium condition 0 < d 2 < γ h , C < 0, depends on γ and d 2 .For both the parameters, d 2 and γ, we have found transcritical and Hopf bifurcation as we varied any one of the parameters at a time, while taking other parameter values to be the same as in Figure 4 (see Figure 7).
To study the stability nature of E * (x * , y * ), we chose the parametric set: {b = 10, 2 and γ = 0.7}, which satisfies the existence and stability criteria of the coexistence steady state E * .We depict the corresponding time series and phase portrait in Figure 8. From here, it is stated that the coexistence steady state E * (x * , y * ) ≡ E * (2, 3.169) is LAS.We check the population sizes of both the prey and predator numerically by varying the fear parameter k 1 and taking other parameters to be the same as in Figure 8 and find that model (3) exhibits a Hopf bifurcation (around E * ) at k [H 1 ] 1 = 0.337 (see Figure 9).Both the prey and predator populations move from a stable nature to an oscillatory nature as the fear level k 1 increases.This outcome is biologically significant, since the prey species is alert and shows signs of habituation after a certain fear threshold.Numerical simulation of the proposed system indicates that the parameter β has the property of taking system (3) from being unstable to stable by destroying Hopf bifurcation around the coexistence steady state and creating a transcritical bifurcation for a comparatively large value of it.We increase the value of β from zero while all other parameters' values remain same as in Figure 8 (see Figure 10).However, α has the opposite nature in the sense that it makes system (3) go from stable to unstable by creating a Hopf bifurcation around E * as we increase it from α [TC] = 0.353 (see Figure 11).Figure 12 illustrates the growth of the predator population with respect to k 2 when the other parameters are same as those in Figure 8.This choice of parameter set suggests that the predator population will decline with the rising value of k 2 .Figure 13 shows that as parameter b increases, the coexistence equilibrium becomes unstable through a Hopf bifurcation, and predator-prey oscillation ensues.We know this phenomenon as the "paradox of enrichment".The parameter b (birth rate of prey) may be visualized as the enrichment parameter, since after simplification of Model (3), we get the expression for carrying capacity as b−d 1 β ; therefore, it is appropriate to consider b as the enrichment parameter while keeping other parametric values unchanged.x Unstable E * b [TC] =2.70 Considering a parametric set as {b = 5.5, 4 and γ = 0.8}, the bifurcation diagram with respect to k 2 is depicted in Figure 14.The figure shows that for a high level of fear (k 2 ), the environmental carrying capacity of prey changes the system's oscillatory behavior to be stable around coexistence equilibrium E * .Figure 15 represents two parametric bifurcation diagrams in the (k 2 − k 1 ) plane, where the solid red line indicates the Hopf curve.In the region above the Hopf-curve, one coexistence steady state (E * ) exists and it is locally asymptotically stable and E * loses its stability through a super-critical Hopf bifurcation when the parametric value passes the Hopf-curve.

Effects of Time Delays on Predator-Prey Population
In this section, we discuss numerically the analytical findings of the delayed system, ( 14).Here, we take the parameter set as: {b = 8.7, 6} and it is seen that when τ 1 = 0 and τ 2 = 0, all the stability conditions (as mention in Theorem 5) of coexistence equilibrium points are satisfied; i.e., the coexistence equilibrium point E * (x * , y * ) ≡ E * (2.5, 1.564) is stable (see Figure 16).

Study of System (3) with Environmental Stochasticity
Environmental fluctuations are not incorporated in deterministic models.However, these are only ecologically beneficial if the dynamical patterns revealed are still in evidence when stochastic influences are incorporated.Environmental stochasticity is generally considered to cause uncertain population birth and mortality rates.Temperature, humidity, parasites and pathogens, environmental pollution, food quality and other factors influence reproduction, growth and death.As these phenomena are difficult to predict, it is preferable to use a stochastic approach instead of a deterministic one.Assume that the environmental fluctuations will manifest themselves primarily as fluctuations in the natural mortality rate of each species, since these are the main parameters subject to coupling of a prey-predator pair with its environment.This parameters are perturbed by Gaussian white noise, which is one of the most useful forms of noise for modeling rapidly fluctuating phenomena.
Therefore, we perturbed the parameters d 1 and d 2 with Gaussian white noises γ 1 and γ 2 in system (3), where γ 1 and γ 2 are independent Gaussian white noises with the following characteristics: Here, µ i > 0 is the intensity or strength of γ i , and the Dirac delta function δ is defined as follows: where • is the ensemble average of the stochastic process under consideration.Thus, system (3) was modified as follows: We took dt , where w = {w 1 (t), w 2 (t) | t ≥ 0} denotes standard Brownian motion in two-dimension.It was assumed that d i + γ i (t) is positive and bounded for i = 1, 2. Thus, from (34), we have the following stochastic system: We already defined the parameters in Section 2. The Euler Maruyama method was used in MATLAB to determine the dynamical behavior of system (35).
In Figures 31 and 32, the impacts of environmental fluctuation on the species are depicted for the parameter set in Figure 8.It is noticeable that, following some initial transients, prey and predator oscillate around the deterministic steady-state values 2 and 3.169, respectively.Figures 31 and 32 indicate that prey and predator species will never go extinct with this parameter set.Hence, system (35) will persist.We have further depicted the stationary distributions of prey and predator populations at t = 300 in Figure 33.By choosing the same parameter values as in Figure 2, we have drawn the stochastic trajectories of prey and predator populations in Figure 34, while taking the intensities of the fluctuations µ 1 = µ 2 = 0.01.Since, for this parameter set, the death rate of the prey population, d 1 = 5.6, is greater than the birth rate, b = 5.5, after certain time, the prey population goes extinct.Consequently, in the absence of food (prey), the predator population also goes extinct (see Figure 34).
In Figure 35, we depict the stochastic trajectories of both populations with the parameter set of Figure 4.As the death rate, d 2 = 0.95, is greater than the birth rate, the predator population cannot persist in the ecosystem, yet the prey species survive properly without predation.

Discussion and Conclusions
One of the key themes in ecology and evolutionary biology is the study of the dynamics of predator-prey interactions, in which the predator consumes the prey population directly.Current field observations of the fear effect in predator-prey dynamics have highlighted the necessity of improving existing systems that do not consider the fear effect.Over the last few years, researchers have been introducing an anti-predation mechanism in mathematical model formulations to take into account the effect of fear, which results in a cost in reproduction.By studying these models, they were able to acquire some conclusions on the impact of such an anti-predation reaction [13,18,26,44,45].Wang et al. [18] first developed a mathematical model that included the fear effect, and they investigated the model using two different functional responses: (i) Holling type I (linear) and (ii) Holling type II (rectangular hyperbola).In this study, we have conducted an analysis of a Holling type-II functional response mediated by the prey's anti-predation reaction.The functions f (k i , y) = 1 1+k i y , where i = 1, 2, are incorporated into the model to account for prey's anti-predation response with positive parameters k i , where i = 1, 2, measuring the level of fear.Clearly, f (k i , y) is decreasing with increasing k i , i = 1, 2 and y, respectively.The analytical results reveal how the cost of fear of prey's anti-predation reaction affects the population dynamics.We also include two different types of delays: the breeding delay of the prey species influenced by fear and the gestation delay of the predator.The dynamical characteristics, feasibility conditions and stability (local and global) criteria for each steady state of the proposed system have been demonstrated.Some numerical simulations are explored to confirm analytical results showing how fear and biomass transfer delay will affect the population dynamics.Fear effects can affect our formulated system in a variety of ways, as indicated by our analytical and numerical results, which can be stated in the following manner: 1.
The predator's death rate, d 2 , and conversion rate, γ, each has an effective impact on the proposed model because the feasibility criteria 0 < d 2 < γ h and C < 0 of the coexistence equilibrium depend on γ and d 2 as mentioned in Equation (5).Consequently, as we vary d 2 and γ individually, both transcritical and Hopf bifurcation can be found (see Figure 7a,b).Theorem 11 states that at some threshold value d 2 , the predator-free equilibrium point is stable; that is, a higher mortality rate of predators leads to the extinction of the species.However, if it is less than the critical point d

[tc]
2 , the coexistence equilibrium point is stable up to some threshold value of d [H] 2 .The prey and predator populations oscillate periodically if d 2 diminishes further.Again, Theorem 12 states that at γ [tc] , the system exhibits transcritical bifurcation.As the value of γ crosses this critical value, the predator-free equilibrium point dies out, and the coexistence steady state becomes stable up to some threshold value, as shown in Figure 7b.The populations of prey and predator oscillate periodically if γ increases further.

2.
Theorem 13 states that system (3) exhibits Hopf bifurcation around coexistence state at some critical value of fear level k 1 .A numerical simulation depicted that as the fear level, k 1 , rises, prey and predator populations exhibit a more stable nature (see Figure 9).Additionally, for fear level k 2 , incorporated into carrying capacity of prey species, it can be observed that a higher value of k 2 can stabilize the proposed model by excluding the existence of periodic solutions, which is also analytically shown in Theorem 14.These studies demonstrate that the prey population does not stop reproduction due to fear.These types of results have also been obtained in the prey-predator models studied by Mondal et al. (2018) [44], Mondal and Samanta (2021) [13] and Wang et al. ( 2016) [18].These findings are significant from a biological perspective due to the fact that the prey species is cautious and shows signs of habituation when a specific threshold of fear has been reached.That is, once the prey population reaches a certain level, the fear no longer has an effect, since the prey are now aware of the predator and exhibit signs of habituation.

3.
At E * (x * , y * ), x * is independent of the parameters k 1 and k 2 .However, y * depends on both parameters.Moreover, dy * dk 1 < 0, so at an interior steady state, successive incrementation of fear level (k 1 ) can decrease the prosperity of the predator population (see Figure 9b).Again, dy * dk 2 < 0; therefore, at E * (x * , y * ), the growth of the predator species also decreases because of continuous incrementation of k 2 (see Figure 12).Such phenomena are ecologically significant because the prey's anti-predation behavior continuously improves as a result of the ongoing improvement in the cost of predation fear (k i , i = 1, 2).Therefore, the predator population cannot get enough prey for their survival.

4.
Numerical simulations show that a Hopf-bifurcation is exhibited around E * (x * , y * ) if we increase the birth rate (b) of the prey species.From our study, it can be observed that coexistence equilibrium E * of system (3) is stable when b [TC] < b < b [H] and is unstable with oscillatory nature when b > b [H] (Figure 13).Thus, the proposed model supports oscillation with a higher birth rate of prey population.That is, the "paradox of enrichment" is visible in the proposed system (3) with parameter b (birth rate of prey population) as the enrichment parameter because when we simplify the system (3), we determine carrying capacity as b−d 1 β ; thus, use of b as the enrichment parameter while keeping other parametric values as constants is appropriate.Figure 13 expresses this phenomenon in various circumstances.5.
Among the analytical findings of delayed system (14), it can be observed that there are critical values of delay parameters τ 1 and τ 2 in all four cases, such that model (14) experiences Hopf bifurcation under some conditions, as stated in Theorems 15-18.Numerical simulations of the delayed system also support this with analytical outcomes.It has been shown that in all the delay cases, as the value of the delayed parameter increases gradually, there is some threshold value such that the model exhibits periodic solution through supercritical Hopf bifurcation around an interior equilibrium point, which has already been discussed in Section 8.1.

6.
Finally, the proposed model ( 3) is compared numerically to a stochastic system (35) incorporating Gaussian white noise terms in the death rates of prey and predator species because of environmental fluctuations.
(a) It can be noticed that, following some initial transients, the biomass values of prey and predator populations for stochastic model (35)  In both systems, (3) and ( 35), if we consider b < d 1 , i.e., if the birth rate (b) of the prey species is less than the death rate (d 1 ), the prey population is not able to survive in ecosystem, and the predator population eventually becomes extinct due to a shortage of food (see Figures 2 and 34).(c) For both models, (3) and (35), if the death rate of the predator d 2 is greater than its birth rate, then the predator species will move towards extinction, but the prey species could persist in the ecosystem (see Figures 4 and 35).
The proposed model has been constructed under certain limitations, such as (i) prey being the only food source for predators, such that the predators cannot sustain themselves without it and will eventually go extinct.In spite of this, when there is a lack of prey, predators will always search for other sources of sustenance in real life.It is plausible that the alternative food source will not be as nutritious or will not appeal to predator as much.Therefore, providing the predator with an extra source of food would represent a step toward a more realistic scenario.Das and Samanta [45] formulated and analysed a stochastic system that includes additional food sources for the predator species and introduced Gaussian white noise to the death rates of prey and predator.As a result, one could upgrade the considered system (3), by incorporating additional food supplies for the predator.(ii) In the context of an ecological environment, carryover effects will occur in any scenario in which a person's previous experiences and history may explain how they are now performing in a particular situation.Here, in the considered system (3), carry-over effects due to fear are not included.Mondal and Samanta [46] investigated the dynamics of predator-prey interplay and the impact of fear (felt by prey) of the predator and its carry-over effects.(iii) Model (3) neglects to take into consideration the advantageous effects of the anti-predation response.Wang and Zou [26] explored the dynamical behavior of a predator-prey system incorporating the cost of fear through adaptive avoidance of predators.Incorporating these facts (additional food sources for the predator, carry-over effects, the benefit of the anti-predation response, adaptive avoidance of predators, etc.) into our model (3) will make it much more realistic.Furthermore, through the use of diffusion-driven instability, it would be possible to develop a plausible mathematical model that could be used to investigate the impact of spatial diffusion on pattern formation.

Figure 3 .
Figure 3. Bifurcation diagram regarding parameter d 1 , and all other parameters are the same as in Figure 2.

Figure 6 . 4 Figure 7 .
Figure 6.Bifurcation diagram regarding parameter h, and all other parameters are the same as in Figure 4.

Figure 11 .Figure 12 .
Figure 11.Bifurcation diagram regarding parameter α when the other parameters are the same as those in Figure 8.(a) Prey population; (b) Predator population.

Figure 15 .
Figure 15.Two parametric bifurcation diagrams in the (k 2 − k 1 ) parametric plane with all other parameters the same as in Figure 14.

Figure 17 .
Figure 17.Bifurcation diagram regarding delay parameter τ 1 when τ 2 = 0 and all other parameter values are the same as in Figure 16.

Figure 20 .
Figure 20.Bifurcation diagram regarding delay parameter τ 2 when τ 1 = 0 and all other parameter values are the same as in Figure 16.

Figure 31 .
Figure 31.Prey population's stochastic trajectory.Here, the prey population fluctuates around the deterministic equilibrium value 2. Parameters µ 1 = µ 2 = 0.001 and all other parameter values are the same as in Figure 8.

Figure 32 .
Figure 32.Predator population's stochastic trajectory.Here, the predator population fluctuates around the deterministic equilibrium value 3.169.Parameters µ 1 = µ 2 = 0.001 and all other parameter values are the same as in Figure 8.

Figure 33 .
Figure 33.Histograms of prey and predator populations corresponding to system (34).

Figure 34 .
Figure 34.Extinction of both species when µ 1 = µ 2 = 0.01 and all other parameter values are the same as in Figure 2.

Figure 35 .
Figure 35.Extinction of the predators when µ 1 = µ 2 = 0.01 and all other parameter values are the same as in Figure 4.
experiences transcritical bifurcation.The predator-free equilibrium point exchanges its stability with an interior equilibrium point at this threshold value,d [tc]2 .The effect of the death rate of the predators, d 2 , has been numerically studied, as shown in Figure7a.Whenever the value of death rate d 2 exceeds its threshold value, d[tc] ) goes to predator-free equilibrium.Thus, the existence criterion of y * is biologically significant.), x * is independent of the parameters k 1 and k 2 .Thus, dx * dk 1 = dx * dk 2 = 0. + αγd 1 (γ − hd 2 ))y * + βγd 2 k 2 + α 2 (γ − hd 2 ) 2 y * 2 B + 2Ay * < 0. The trivial steady state E 0 (0, 0) is locally asymptotically stable (LAS) if b < d 1 and unstable if b > d 1 .