1. 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
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 delay-induced 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.
2. 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;
(
) is the conversion rate of the predator;
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
[
34,
35,
36].
Now, let
, where
b and
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
by
and the carrying capacity
k by
, where
and
satisfy the following conditions [
18]:
Here, the parameter
, where
, describes the level of fear that drives the prey’s anti-predation behavior. In particular, let
and
. Thus, by incorporating the fear effects in model (
2), we have obtained the following predator–prey system:
where
,
and
. Throughout the analysis we will take
, if not stated otherwise. The description and unit of the model parameters are presented in
Table 1.
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].
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.
In [
38], the authors presumed and analysed the dynamics of a predator–prey interplay by including a nonlinear prey refuge function
, where
(>0) is a half-saturation constant for refuge function and
(>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.
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
,
, the fear function
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
, where
, 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.
7. 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 non-delayed 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,
(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 (
) that governs the reproduction of predator population following prey hunting. To acquire the rich dynamics of (
3), we introduce both these delays
and
in system (
3). Therefore, after incorporating delays, the modified form of system (
3) is:
Let
C denote the Banach space of continuous functions
equipped with supremum norm
with
.
By linearizing system (
14) about
, we have:
where
are
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
, 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
will be derived using the same analytic arguments mentioned in [
42,
43], in the presence of both the time delays.
.
In absence of both the time delays, system (
14) is reduced to system (
3). We have already described the condition of the stability of
in Theorem 5 in absence of
and
.
, .
Therefore, the characteristic equation, Equation (
16), becomes
where
,
,
,
.
For
, Theorem 5 gives the criteria for which all the zeros of Equation (
17) are negative or have negative real parts, whereas for
, it has an infinite number of roots. According to Rouche’s theorem and continuity in
, the signs of the roots of Equation (
17) will change if they cross the imaginary axis. Thus, by putting
in (
17) and separating the real and imaginary parts, we get
By squaring and adding Equations (
18) and (19), we get
By substituting
, we get the following equation in
:
If there is no positive real root in Equation (
21), real
does not exist. Thus, for any
,
is locally asymptotically stable. If Equation (
21) has a positive real root, say
such that
. Thus,
are two purely imaginary roots of (
17).
From (
18) and (19), we get the values of
as:
for
. Thus, if Theorem 5 holds, system (
14) will be LAS around the interior equilibrium
for
. Therefore, by Butler’s lemma,
is stable for
, where
, and
is unstable for
, provided the transversality condition is satisfied. We now verify the transversality condition
. By differentiating Equation (
17) with respect to
, we get:
This gives
Now,
Therefore, the transversality condition is satisfied, and so Hopf bifurcation occurs at
if
. The following theorem can now be stated.
Theorem 15. Suppose that equilibrium exists and is LAS for . Additionally, let be a positive root of Equation (21). Then, such that equilibrium of system (14) is LAS when and unstable for , wherefor . Furthermore, system (14) will exhibit a Hopf bifurcation around when provided . , .
Here, we consider the case when breeding delay (
) is absent but gestation delay (
) is present. Therefore, the characteristic equation, Equation (
16), becomes
where
,
,
and
. We can state the following theorem by relying on the same analysis as in Case II.
Theorem 16. Suppose that equilibrium exists and is LAS for . Additionally, let when . Then, such that equilibrium of system (14) is LAS when and unstable for , wherefor . Furthermore, system (14) will exhibit a Hopf bifurcation around when , provided . , .
Assume
as a fixed number in the interval
and
. Let
be the root of the characteristic equation, Equation (
16). By putting this value into Equation (
16), we have
comparing real and imaginary parts,
and
The root of the characteristic equation, Equation (
16), must be purely imaginary for equilibrium
to change in stability. Therefore, by putting
in Equations (
25) and (
26), we get
By squaring and adding Equations (
27) and (
28), we have
If
, Equation (
29) has a positive root, say,
. Consequently, Equation (
16) has a pair of imaginary roots
for some fixed
. Therefore, by solving Equations (
27) and (
28) for
at
, we get:
for
, where
,
,
,
and
.
Now, if equilibrium
is LAS in absence of both the time delays, then by Butler’s lemma, it will remain stable for
, such that
and system (
14) undergoes Hopf bifurcation (around
) at
, provided the transversality condition,
, is satisfied.
By differentiating Equations (
25) and (
26) with respect to
, we get:
where
By solving Equations (
31) and (
32), we get:
Therefore, Hopf bifurcation (around
) occurs at
if
.
Theorem 17. If equilibrium exists and is LAS in absence of both the time delays, then for if , there exists such that equilibrium is LAS when and unstable when , wherefor . Further, system (14) will exhibit a Hopf bifurcation around for , provided . , .
We can state the following theorem relying on the same analysis as in Case IV.
Theorem 18. If equilibrium exists and is LAS in absence of both the time delays, then for if , there exists such that equilibrium is LAS when and unstable when , wherefor . Further, system (14) will exhibit a Hopf bifurcation around equilibrium for , provided where 8. 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:
,
,
,
,
,
,
,
and
. Here, we see that the death rate,
, 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
from
, while leaving all other parameters to be the same as in
Figure 2, and we can see that at
, i.e., when
, one eigenvalue is zero, and the stability of trivial equilibrium point,
, exchanges with that of the axial equilibrium point,
, i.e., transcritical bifurcation occurs. If we further decrease
, at
, the axial equilibrium point exchanges its stability with that of the positive interior equilibrium point
, and finally at
, system (
3) undergoes stable Hopf bifurcation around the positive interior equilibrium point (see
Figure 3).
To study the stability nature of the predator-free (axial) equilibrium point, we took the parameter set:
,
,
,
,
,
,
,
and
. Here,
and the stability condition
is satisfied. Therefore,
is stable. 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
,
exchanges its stability with that of
(see
Figure 5a). We have seen a similar type of behavior when we increase the parameter value of
from zero onward up to
; the stability of
exchanges with that of
(see
Figure 5b). In the case of parameter
h, the stability of
exchanges with that of the equilibrium point
at
(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
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
,
, depends on
and
. For both the parameters,
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
, we chose the parametric set:
,
,
,
,
,
,
,
and
, which satisfies the existence and stability criteria of the coexistence steady state
. We depict the corresponding time series and phase portrait in
Figure 8. From here, it is stated that the coexistence steady state
is LAS. We check the population sizes of both the prey and predator numerically by varying the fear parameter
and taking other parameters to be the same as in
Figure 8 and find that model (
3) exhibits a Hopf bifurcation (around
) at
(see
Figure 9). Both the prey and predator populations move from a stable nature to an oscillatory nature as the fear level
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
as we increase it from
(see
Figure 11).
Figure 12 illustrates the growth of the predator population with respect to
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
.
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
; therefore, it is appropriate to consider
b as the enrichment parameter while keeping other parametric values unchanged.
Considering a parametric set as {
,
,
,
,
,
and
}, the bifurcation diagram with respect to
is depicted in
Figure 14. The figure shows that for a high level of fear
, the environmental carrying capacity of prey changes the system’s oscillatory behavior to be stable around coexistence equilibrium
.
Figure 15 represents two parametric bifurcation diagrams in the
plane, where the solid red line indicates the Hopf curve. In the region above the Hopf-curve, one coexistence steady state (
) exists and it is locally asymptotically stable and
loses its stability through a super-critical Hopf bifurcation when the parametric value passes the Hopf-curve.
8.1. 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:
and it is seen that when
and
, all the stability conditions (as mention in Theorem 5) of coexistence equilibrium points are satisfied; i.e., the coexistence equilibrium point
is stable (see
Figure 16).
For the delay (Case II), when
and
, we have analytically found that delayed system, (
14), exhibits a Hopf bifurcation at some critical value of
. After choosing the same parameter set as in
Figure 16, we have drawn the corresponding Hopf bifurcation diagram (see
Figure 17). For
Figure 18, we have taken
, and we can see that the delayed system, (
14), is stable, whereas in
Figure 19, we can observe the oscillatory nature of delayed system (
14) around
for
. Therefore, the delay parameter
can make system (
14) go from stable to unstable as we increase the value of
from zero onward while taking
.
For Case III: when
and
, we found that delayed system (
14) undergoes a Hopf bifurcation for some critical value of
, as depicted in
Figure 20. Again, we took
and
to draw the time series. We can observe the stability and oscillatory nature of system (
14), as depicted in
Figure 21 and
Figure 22, respectively.
Figure 23 is the Hopf bifurcation diagram with the delay parameter
when
and Hopf bifurcation occurs at
. Additionally, we have drawn the time series of the delayed system, (
14), with
(see
Figure 24),
(see
Figure 25) and fixed
.
It can be observed in
Figure 24 and
Figure 25 that when
, the system is stable for
and unstable otherwise. However, numerically it can be noticed that when
, for any positive value of
, the system is always unstable and forms a limit cycle around the coexistence state
(see
Figure 26). As a consequence, system (
14) does not exhibit any switching behavior for any
, regardless of the value of
.
Finally, in Case V: for
and
, we have drawn the bifurcation diagram with
. It can be observed that system (
14) exhibits a Hopf bifurcation at some critical point
(see
Figure 27). We have drawn the time series for system (
14) with
(see
Figure 28) and
(see
Figure 29), which is consistent with the bifurcation diagram in
Figure 27.
Figure 28 and
Figure 29 illustrate that the system is stable for
provided
, and unstable otherwise. However, it is evident from numerical simulations that an unstable limit cycle is generated when
for any positive value of
(see
Figure 30). This leads to the conclusion that system (
14) has no switching behavior for any
, irrespective of the values of
.
8.2. 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
and
with Gaussian white noises
and
in system (
3), where
and
are independent Gaussian white noises with the following characteristics:
Here,
is the intensity or strength of
, 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
and
, where
denotes standard Brownian motion in two-dimension. It was assumed that
is positive and bounded for
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
Figure 31 and
Figure 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
, respectively.
Figure 31 and
Figure 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
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
. Since, for this parameter set, the death rate of the prey population,
, is greater than the birth rate,
, 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,
, is greater than the birth rate, the predator population cannot persist in the ecosystem, yet the prey species survive properly without predation.
9. 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
, where
, are incorporated into the model to account for prey’s anti-predation response with positive parameters
, where
, measuring the level of fear. Clearly,
is decreasing with increasing
,
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,
, and conversion rate,
, each has an effective impact on the proposed model because the feasibility criteria
and
of the coexistence equilibrium depend on
and
as mentioned in Equation (
5). Consequently, as we vary
and
individually, both transcritical and Hopf bifurcation can be found (see
Figure 7a,b). Theorem 11 states that at some threshold value
, the system experiences transcritical bifurcation. The predator-free equilibrium point exchanges its stability with an interior equilibrium point at this threshold value,
. The effect of the death rate of the predators,
, has been numerically studied, as shown in
Figure 7a. Whenever the value of death rate
exceeds its threshold value,
, 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
, the coexistence equilibrium point is stable up to some threshold value of
. The prey and predator populations oscillate periodically if
diminishes further. Again, Theorem 12 states that at
, 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
. A numerical simulation depicted that as the fear level,
, rises, prey and predator populations exhibit a more stable nature (see
Figure 9). Additionally, for fear level
, incorporated into carrying capacity of prey species, it can be observed that a higher value of
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
,
is independent of the parameters
and
. However,
depends on both parameters. Moreover,
, so at an interior steady state, successive incrementation of fear level (
) can decrease the prosperity of the predator population (see
Figure 9b). Again,
; therefore, at
, the growth of the predator species also decreases because of continuous incrementation of
(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
. Therefore, the predator population cannot get enough prey for their survival.
- (4)
Numerical simulations show that a Hopf-bifurcation is exhibited around
if we increase the birth rate
of the prey species. From our study, it can be observed that coexistence equilibrium
of system (
3) is stable when
and is unstable with oscillatory nature when
(
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
; 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
and
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) oscillate around the deterministic equilibrium values of 2 and
, respectively (see
Figure 31 and
Figure 32). Additionally, these two figures (
Figure 31 and
Figure 32) depict that for both systems (
35) (stochastic system) and (
3) (deterministic model), neither predator nor prey goes extinct.
- (b)
In both systems, (
3) and (
35), if we consider
, i.e., if the birth rate
of the prey species is less than the death rate
, the prey population is not able to survive in ecosystem, and the predator population eventually becomes extinct due to a shortage of food (see
Figure 2 and
Figure 34).
- (c)
For both models, (
3) and (
35), if the death rate of the predator
is greater than its birth rate, then the predator species will move towards extinction, but the prey species could persist in the ecosystem (see
Figure 4 and
Figure 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.