Next Article in Journal
Risk Evaluation Model of Coal Spontaneous Combustion Based on AEM-AHP-LSTM
Next Article in Special Issue
COVID-19: From Limit Cycle to Stable Focus
Previous Article in Journal
Deep Spatial-Temporal Neural Network for Dense Non-Rigid Structure from Motion
Previous Article in Special Issue
Mathematical Model of Pancreatic Cancer Cell Dynamics Considering the Set of Sequential Mutations and Interaction with the Immune System
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India
*
Author to whom correspondence should be addressed.
Mathematics 2022, 10(20), 3795; https://doi.org/10.3390/math10203795
Submission received: 14 August 2022 / Revised: 27 September 2022 / Accepted: 4 October 2022 / Published: 14 October 2022
(This article belongs to the Special Issue Advances in the Mathematics of Ecological Modelling)

Abstract

:
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 influenced 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 findings, we also perform some numerical computations and simulations. Finally, to evaluate the influence of a fluctuating environment, we compare our proposed system to a stochastic model with Gaussian white noise terms.

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 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 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:
d x d t = r x ( 1 x k ) α x y 1 + α h x , x ( 0 ) > 0 d y d t = γ α x y 1 + α h x d 2 y , y ( 0 ) > 0
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 + α h x [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:
d x d t = b x d 1 x r x 2 k α x y 1 + α h x , x ( 0 ) > 0 d y d t = γ α x y 1 + α h x d 2 y , y ( 0 ) > 0 .
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 b x 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]:
f i ( k i , 0 ) = 1 , f i ( 0 , y ) = 1 , f i ( k i , y ) k i < 0 , f i ( k i , y ) y < 0 , for i = 1 , 2 .
Here, the parameter k i , where i = 1 , 2 , describes the level of fear that drives the prey’s anti-predation 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:
d x d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x , x ( 0 ) > 0 d y d t = γ α x y 1 + α h x d 2 y , y ( 0 ) > 0
where β = r k , b 1 + k 1 y b and k 1 + k 2 y k . 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.
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 ϕ ( x , y ) = m 1 x y a 0 + y , where a 0 (>0) is a half-saturation 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.
  • 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.

3. 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 .
Proof. 
Let
Ψ 1 ( x , y ) = b 1 + k 1 y d 1 β ( 1 + k 2 y ) x α y 1 + α h x , Ψ 2 ( x , y ) = γ α x 1 + α h x d 2 .
Then, system (3) reduces to:
d x d t = x Ψ 1 ( x , y ) , x ( 0 ) > 0 d y d t = y Ψ 2 ( x , y ) , y ( 0 ) > 0
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
x ( t ) = x ( 0 ) exp 0 t Ψ 1 ( x ( s ) , y ( s ) ) d s > 0 , [ x ( 0 ) > 0 ] y ( t ) = y ( 0 ) exp 0 t Ψ 2 ( x ( s ) , y ( s ) ) d s > 0 , [ y ( 0 ) > 0 ]
Hence, the theorem is proved. □
Theorem 2.
All solutions of system (3) are uniformly bounded in the region B = { ( x , y ) : 0 < x ( t ) + y ( t ) γ b 2 4 β μ + δ , f o r a n y δ > 0 } if b > d 1 .
Proof. 
By taking the first equation of Model (3), we get
d x d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x b x d 1 x β x 2 α x y 1 + α h x , b 1 + k 1 y b ( b d 1 ) x β x 2 .
lim sup t x ( t ) b d 1 β , p r o v i d e d b > d 1 . 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,
d Ω d t = d x d t + 1 γ d y d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 d 2 y γ b x d 1 x β x 2 d 2 y γ , b 1 + k 1 y b = β x b 2 β 2 + b 2 4 β d 1 x d 2 y γ b 2 4 β d 1 x d 2 y γ b 2 4 β μ Ω , where   μ = min { d 1 , d 2 } . d Ω d t + μ Ω b 2 4 β .
Therefore, by Gronwall inequality, we have
0 < Ω ( x ( t ) , y ( t ) ) b 2 4 β μ 1 e μ t + Ω ( x ( 0 ) , y ( 0 ) ) e μ t
lim t x ( t ) + y ( t ) γ b 2 4 β μ + δ ,   for   any δ > 0 .
Therefore, all solutions of (3) enter in the region:
B = ( x , y ) : 0 < x ( t ) + y ( t ) γ b 2 4 β μ + δ , for   any δ > 0 .
This completes the proof. □

4. Equilibria and Their Stability Analysis

In this section, we evaluate the proposed system’s behavior analytically.

4.1. Equilibrium Points

System (3) has three feasible equilibrium points or steady states in R + 2 , namely:
(i)
E 0 ( 0 , 0 ) (trivial equilibrium),
(ii)
E 1 b d 1 β , 0 (predator-free equilibrium) exists if b > d 1 ,
(iii)
E * ( x * , y * ) (coexistence equilibrium) where x * = d 2 α ( γ h d 2 ) provided 0 < d 2 < γ h and y * is the unique positive root of the equation:
A y 2 + B y + C = 0 , provided   C < 0 .
Here,
A = k 1 k 2 β γ d 2 + α 2 k 1 ( γ h d 2 ) 2 ( > 0 ) , B = β γ d 2 ( k 1 + k 2 ) + d 1 k 1 α γ ( γ h d 2 ) + α 2 ( γ h d 2 ) 2 ( > 0 ) , C = β γ d 2 γ α ( b d 1 ) ( γ h d 2 ) .
Thus, we have
y * = B + B 2 4 A C 2 A .
These three equilibrium points are represented in Figure 1 using the parameter set mentioned in the figure caption.
Remark 1.
If C > 0 , i.e., if d 2 > γ α b d 1 β 1 + α h b d 1 β , the coexisting equilibrium E * ( x * , y * ) goes to predator-free equilibrium. Thus, the existence criterion of y * is biologically significant.
Remark 2.
At E * ( x * , y * ) , x * is independent of the parameters k 1 and k 2 . Thus, d x * d k 1 = d x * d k 2 = 0 . However, y * depends on both parameters. Moreover,
d y * d k 1 = β γ d 2 + α γ d 1 ( γ h d 2 ) y * + β γ d 2 k 2 + α 2 ( γ h d 2 ) 2 y * 2 B + 2 A y * < 0 .
Therefore, at the coexistence equilibrium point E * , increasing the level of fear ( k 1 ) can lower the size of the predator population.
Again,
d y * d k 2 = β γ d 2 y * + β γ d 2 k 1 y * 2 B + 2 A y * < 0 .
Therefore, at E * ( x * , y * ) , the growth of predator population decreases due to continuous increment in k 2 .

4.2. Local Stability Analysis

Theorem 3.
The trivial steady state E 0 ( 0 , 0 ) is locally asymptotically stable (LAS) if b < d 1 and unstable if b > d 1 .
Proof. 
The Jacobian matrix J 0 , 0 at E 0 ( 0 , 0 ) is
J 0 , 0 = b d 1 0 0 d 2 .
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 1 b d 1 β , 0 is locally asymptotically stable if γ α b d 1 β 1 + α h b d 1 β < d 2 and unstable if γ α b d 1 β 1 + α h b d 1 β > d 2 .
Proof. 
The Jacobian matrix corresponding to E 1 b d 1 β , 0 is given by:
J b d 1 β , 0 = ( b d 1 ) b k 1 ( b d 1 ) β k 2 ( b d 1 ) 2 β α ( b d 1 ) β + α h ( b d 1 ) 0 γ α b d 1 β 1 + α h b d 1 β d 2 .
Thus, the eigenvalues of the Jacobian matrix J b d 1 β , 0 are ( b d 1 ) and γ α b d 1 β 1 + α h b d 1 β d 2 . Therefore, E 1 b d 1 β , 0 is LAS if γ α b d 1 β 1 + α h b d 1 β < d 2 and unstable if γ α b d 1 β 1 + α h b d 1 β > d 2 . □
Theorem 5.
The coexistence steady state E * ( x * , y * ) is LAS if β x * ( 1 + k 2 y * ) + α 2 h x * y * ( 1 + α h x * ) 2 < 0 .
Proof. 
The Jacobian matrix J x * , y * at the coexistence steady state E * ( x * , y * ) is given by:
J x * , y * = j 11 j 12 j 21 j 22 ,
where
j 11 = β x * ( 1 + k 2 y * ) + α 2 h x * y * ( 1 + α h x * ) 2 , j 12 = x * b k 1 ( 1 + k 1 y * ) 2 + β k 2 x * + α ( 1 + α h x ) ( < 0 ) , j 21 = γ α y * ( 1 + α h x * ) 2 ( > 0 ) , j 22 = 0 .
Therefore, all the eigenvalues of J ( x * , y * ) will have negative real parts if j 11 < 0 . Thus, E * ( x * , y * ) is LAS if j 11 = β x * ( 1 + k 2 y * ) + α 2 h x * y * ( 1 + α h x * ) 2 < 0 . □

4.3. Global Stability Analysis

Theorem 6.
The trivial steady state E 0 ( 0 , 0 ) is globally asymptotically stable (GAS) if b < d 1 .
Proof. 
Consider a Lyapunov function:
V ( x , y ) = x ( t ) + y ( t ) γ .
Then,
d V d t = d x d t + 1 γ d y d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 d 2 y γ b x d 1 x β x 2 d 2 y γ , b 1 + k 1 y b = ( d 1 b ) x β x 2 d 2 y γ .
Thus, d V d t < 0 if b < d 1 and d V d t = 0 at E 0 ( 0 , 0 ) . Therefore, using Lyapunov theorem, we can conclude that the steady state E 0 ( 0 , 0 ) is GAS if b < d 1 . □
Theorem 7.
The axial steady state E 1 b d 1 β , 0 is GAS if d 2 > γ h .
Proof. 
From the first equation of system (3),
d x d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 b x d 1 x β x 2 , b 1 + k 1 y b = ( b d 1 ) x 1 x b d 1 β lim sup t x ( t ) b d 1 β .
Now, second equation of (3) gives:
d y d t = γ α x y 1 + α h x d 2 y d 2 γ h y
lim t y ( t ) = 0 i f d 2 > γ h .
Additionally, x ( t ) b d 1 β as y ( t ) 0 .
Theorem 8.
The steady state E * ( x * , y * ) of system (3) is GAS if 0 < γ h d 2 h b and 2 h β d 2 ( α h b β ) ( γ h d 2 ) 0 .
Proof. 
Let
Φ 1 ( x , y ) = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x Φ 2 ( x , y ) = γ α x y 1 + α h x d 2 y .
Let us take Dulac function as:
χ ( x , y ) = x 1 y m ( 1 + k 1 y ) ( 1 + α h x ) ,
where m is to be specified later. Let us now calculate the divergence of the vector:
D ( x , y ) = x χ ( x , y ) Φ 1 ( x , y ) + y χ ( x , y ) Φ 2 ( x , y ) x 1 y m { f 1 ( x , m ) k 1 y + f 2 ( x , m ) } ,
where
f 1 ( x , m ) = 2 α h β x 2 + { α ( m + 2 ) ( γ h d 2 ) β } x d 2 ( m + 2 ) , f 2 ( x , m ) = 2 α h β x 2 + { α ( m + 1 ) ( γ h d 2 ) β + α h b } x d 2 ( m + 1 ) .
Now,
f 2 ( x , m ) f 1 ( x , m ) = α ( γ h d 2 ) + α h b x + d 2 > 0 f o r x [ 0 , )
if 0 < γ h d 2 h b . Therefore, D ( x , y ) < 0 , ( x , y ) R + 2 if
f 2 ( x , m ) 0   for   x [ 0 , ) .
Therefore, it is sufficient to find a m so that condition (6) is satisfied. However, condition (6) holds if
g ( m + 1 ) = : α ( m + 1 ) ( γ h d 2 ) β + α h b 2 8 α h β d 2 ( m + 1 ) 0 .
For convenience, let m + 1 = m ˜ . Then, (7) becomes:
g ( m ˜ ) = α m ˜ ( γ h d 2 ) β + α h b 2 8 α h β d 2 m ˜ < 0 ,
g ( m ˜ ) = α 2 ( γ h d 2 ) 2 m ˜ 2 + 2 ( α h b β ) ( γ h d 2 ) 4 h β d 2 α m ˜ + ( α h b β ) 2 < 0 .
The existence of a m ˜ satisfying (8) means:
g α ( α h b β ) ( γ h d 2 ) 4 α h β d 2 α 2 ( γ h d 2 ) 2 0 ,
which is equivalent to
α ( α h b β ) ( γ h d 2 ) 4 α h β d 2 2 α 2 ( γ h d 2 ) 2 ( α h b β ) 2 0 , i . e . , 2 h β d 2 ( α h b β ) ( γ h d 2 ) 0 .
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 < γ h d 2 h b and 2 h β d 2 ( α h b β ) ( γ h d 2 ) 0 . Hence, E * is globally asymptotically stable. □

5. Uniform Persistence

Theorem 9.
System (3) is uniformly persistent if b > d 1 and γ α b d 1 β 1 + α h b d 1 β > d 2 .
Proof. 
Let us consider the following average Lyapunov function:
ϱ ( x , y ) = x ϱ 1 y ϱ 2 , ( x , y ) R + 2
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:
ϱ ˙ ϱ = ϱ 1 x ˙ x + ϱ 2 y ˙ y
Λ = ϱ 1 b 1 + k 1 y d 1 β ( 1 + k 2 y ) x α y 1 + α h x + ϱ 2 γ α x 1 + α h x d 2
where Λ = ϱ ˙ ϱ . By computing Λ at E 0 ( 0 , 0 ) and E 1 b d 1 β , 0 , we get:
Λ ( 0 , 0 ) = ϱ 1 ( b d 1 ) ϱ 2 d 2 ,
Λ b d 1 β , 0 = ϱ 2 γ α b d 1 β 1 + α h b d 1 β d 2 .
We choose ϱ 1 and ϱ 2 in such a way that Λ ( 0 , 0 ) > 0 provided b > d 1 . Again, Λ b d 1 β , 0 > 0 , if γ α b d 1 β 1 + α h b d 1 β > d 2 . Therefore, system (3) is uniformly persistent [41] if b > d 1 and γ α b d 1 β 1 + α h b d 1 β > d 2 . □
Remark 3.
Persistency implies instability of the boundary equilibria.

6. Local Bifurcations

6.1. 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 [ t c ] = b .
Proof. 
Let d 1 = d 1 [ t c ] be the critical value of d 1 for which exactly one eigenvalue of J ( 0 , 0 ) is zero. Thus, d 1 [ t c ] = b . Now, the eigenvectors of J ( 0 , 0 ) and [ J ( 0 , 0 ) ] T for the zero eigenvalue are U = ( u 1 , u 2 ) T = ( 1 , 0 ) T and W = ( w 1 , w 2 ) T = ( 1 , 0 ) T , respectively. Let us calculate
Δ 1 = W T . G d 1 0 , 0 ; d 1 [ t c ] = ( 1 , 0 ) . G 1 d 1 G 2 d 1 E 0 = ( 1 , 0 ) . x 0 E 0 = 0 ,
where G = ( G 1 , G 2 ) T ; G 1 = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x , G 2 = γ α x y 1 + α h x d 2 y .
Δ 2 = W T . D G d 1 0 , 0 ; d 1 [ t c ] U = ( 1 , 0 ) . 2 G 1 x d 1 2 G 1 y d 1 2 G 2 x d 1 2 G 2 y d 1 E 0 . 1 0 = 1 0 .
Δ 3 = W T . D 2 G E 0 ; d 1 [ t c ] U , U = ( 1 , 0 ) . D G 1 x u 1 + G 1 y u 2 G 2 x u 1 + G 2 y u 2 E 0 . u 1 u 2 = 2 β < 0 .
Therefore, by Sotomayor’s theorem, a transcritical bifurcation occurs around E 0 at d 1 [ t c ] = b . □
Theorem 11.
System (3) exhibits a transcritical bifurcation around the axial equilibrium point E 1 b d 1 β , 0 with a bifurcating parameter d 2 at d 2 [ t c ] = γ α b d 1 β 1 + α h b d 1 β .
Proof. 
Let d 2 = d 2 [ t c ] be the critical value of d 2 for which J b d 1 β , 0 has exactly one zero eigenvalue at d 2 = d 2 [ t c ] . Thus, d 2 [ t c ] = γ α b d 1 β 1 + α h b d 1 β . Now, the eigenvectors of J b d 1 β , 0 and J b d 1 β , 0 T for the zero eigenvalue are U = ( u 1 , u 2 ) T and W = ( 0 , 1 ) T , respectively. Here, u 2 = 1 and u 1 = b k 1 ( b d 1 ) β k 2 ( b d 1 ) 2 β α ( b d 1 ) β + α h ( b d 1 ) b d 1 . Let us calculate Δ 1 , Δ 2 , Δ 3 :
Δ 1 = W T . G d 2 0 , 0 ; d 2 [ t c ] = ( 0 , 1 ) . G 1 d 2 G 2 d 2 E 1 = ( 0 , 1 ) . 0 y E 1 = 0 ,
where G = ( G 1 , G 2 ) T ; G 1 = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x , G 2 = γ α x y 1 + α h x d 2 y .
Δ 2 = W T . D G d 2 0 , 0 ; d 2 [ t c ] U = ( 0 , 1 ) . 2 G 1 x d 2 2 G 1 y d 2 2 G 2 x d 2 2 G 2 y d 2 E 1 . u 1 u 2 = 1 0 . Δ 3 = W T . D 2 G E 0 ; d 2 [ t c ] U , U = ( 0 , 1 ) . D G 1 x u 1 + G 1 y u 2 G 2 x u 1 + G 2 y u 2 E 1 . u 1 u 2 = 2 γ α 1 + α h b d 1 β u 1 < 0 , u 1 = b k 1 ( b d 1 ) β k 2 ( b d 1 ) 2 β α ( b d 1 ) β + α h ( b d 1 ) b d 1 < 0 .
Therefore, by Sotomayor’s theorem, a transcritical bifurcation is exhibited around E 1 by taking d 2 as the bifurcating parameter. □
Theorem 12.
System (3) exhibits a transcritical bifurcation around the axial equilibrium point E 1 b d 1 β , 0 with bifurcating parameter γ at γ [ t c ] = d 2 1 + β α ( b d 1 ) .
Proof. 
The proof is similar to the proof of Theorem 11. □

6.2. 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
λ 2 T λ + D = 0 ,
where T = trace of J ( x * , y * ) = β x * ( 1 + k 2 y * ) + α 2 h x * y * ( 1 + α h x * ) 2 and D = determinant of J ( x * , y * ) = γ α x * y * ( 1 + α h x * ) 2 b k 1 ( 1 + k 1 y * ) 2 + β k 2 x * + α ( 1 + α h x ) .
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 [ H ] , provided T ( k 1 [ H ] ) = 0 , D ( k 1 [ H ] ) > 0 and d T d k 1 k 1 = k 1 [ H ] 0 .
Proof. 
Suppose k 1 = k 1 [ H ] is the critical value of k 1 at which Equation (10) gives purely imaginary roots, so T ( k 1 [ H ] ) = 0 and D ( k 1 [ H ] ) > 0 . Therefore, in an open neighborhood of k 1 [ H ] , the zeros of Equation (10) take the following form: λ 1 = p ( k 1 ) + i q ( k 1 ) and λ 2 = p ( k 1 ) i q ( k 1 ) , where p ( k 1 ) and q ( k 1 ) are real valued functions of k 1 . By putting λ ( k 1 ) = p ( k 1 ) + i q ( k 1 ) in Equation (10) and differentiating with respect to k 1 , we have
2 p ( k 1 ) + i q ( k 1 ) p ˙ ( k 1 ) + i q ˙ ( k 1 ) T ˙ ( k 1 ) p ( k 1 ) + i q ( k 1 ) T ( k 1 ) p ˙ ( k 1 ) + i q ˙ ( k 1 ) + D ˙ ( k 1 ) = 0 .
By separating the real and imaginary parts, we get
X 1 p ˙ ( k 1 ) X 2 q ˙ ( k 1 ) + X 3 = 0
X 2 p ˙ ( k 1 ) + X 1 q ˙ ( k 1 ) + X 4 = 0 ,
where
X 1 = 2 p ( k 1 ) T ( k 1 ) X 2 = 2 q ( k 1 ) X 3 = p ( k 1 ) T ˙ ( k 1 ) + D ˙ ( k 1 ) X 4 = T ˙ ( k 1 ) q ( k 1 ) .
By solving Equations (11) and (12) for p ˙ ( k 1 ) , we get:
p ˙ ( k 1 ) = X 1 X 3 + X 2 X 4 X 1 2 + X 2 2 .
Now, for k 1 = k 1 [ H ] , we consider the following cases:
  • Case-1
p = 0 , q = D . Therefore, X 1 = 0 , X 2 = 2 D , X 3 = D ˙ , X 4 = T ˙ D . Hence, from Equation (13), we have
p ˙ k 1 = k 1 [ H ] = d p ( k 1 ) d k 1 k 1 = k 1 [ H ] = 1 2 d T d k 1 k 1 = k 1 [ H ] .
  • Case-2
p = 0 , q = D . Therefore, X 1 = 0 , X 2 = 2 D , X 3 = D ˙ , X 4 = T ˙ D . Hence, from Equation (13), we have:
p ˙ k 1 = k 1 [ H ] = d p ( k 1 ) d k 1 k 1 = k 1 [ H ] = 1 2 d T d k 1 k 1 = k 1 [ H ] .
According to the Hopf bifurcation theorem [3], system (3) switches its stability behavior through Hopf bifurcation, provided
d d k 1 ( R e λ ( k 1 ) ) k 1 = k 1 [ H ] = d p ( k 1 ) d k 1 k 1 = k 1 [ H ] 0
is satisfied.
Thus, for the existence of Hopf bifurcation at k 1 = k 1 [ H ] , we must have
d T d k 1 k 1 = k 1 [ H ] 0 .
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 2 [ H ] , provided T ( k 2 [ H ] ) = 0 , D ( k 2 [ H ] ) > 0 and d T d k 2 k 2 = k 2 [ H ] 0 .
Proof. 
The proof is similar to the proof of Theorem 13. □

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, τ 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:
d x d t = b x ( t τ 1 ) 1 + k 1 y ( t τ 1 ) d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x , d y d t = γ α x ( t τ 2 ) y ( t τ 2 ) 1 + α h x ( t τ 2 ) d 2 y ,
Let C denote the Banach space of continuous functions ψ : ξ , 0 R + 2 equipped with supremum norm
ψ = sup ξ s 0 ψ 1 s , ψ 2 s x ( s ) = ψ 1 ( s ) , y ( s ) = ψ 2 ( s ) , ψ i ( s ) > 0 , s ξ , 0 , ψ i ( 0 ) > 0 , i = 1 , 2 .
with ξ = max { τ 1 , τ 2 } .
By linearizing system (14) about E * ( x * , y * ) , we have:
d X d t = A X ( t ) + B X ( t τ 1 ) + C X ( t τ 2 )
where A , B , C are 2 × 2 matrices given by
A = a 11 a 12 0 a 22 , B = b 11 b 12 0 0 , C = 0 0 c 21 c 22 ,
where
a 11 = d 1 2 β x * ( 1 + k 2 y * ) α y * ( 1 + α h x * ) 2 , a 12 = β k 2 x * 2 α x * 1 + α h x * , a 22 = d 2 , b 11 = b 1 + k 1 y * , b 12 = k 1 b x * ( 1 + k 1 y * ) 2 , c 21 = γ α y * ( 1 + α h x * ) 2 , c 22 = γ α x * 1 + α h x * .
The characteristic equation of the delay system is
λ 2 + A 1 λ + A 2 + ( B 1 λ + B 2 ) e λ τ 1 + ( C 1 λ + C 2 ) e λ τ 2 + D 1 e λ ( τ 1 + τ 2 ) = 0
where
A 1 = ( a 11 + a 22 ) , A 2 = a 11 a 22 , B 1 = b 11 , B 2 = a 22 b 11 , C 1 = c 22 , C 2 = a 11 c 22 a 12 c 21 , D 1 = b 11 c 22 b 12 c 21 .
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
τ 1 = τ 2 = 0 .
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
τ 1 > 0 , τ 2 = 0 .
Therefore, the characteristic equation, Equation (16), becomes
λ 2 + M 1 λ + M 2 + ( N 1 λ + N 2 ) e λ τ 1 = 0 ,
where M 1 = A 1 + C 1 , M 2 = A 2 + C 2 , N 1 = B 1 , N 2 = B 2 + D 1 .
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
N 2 cos ω τ 1 + N 1 ω sin ω τ 1 = ω 2 M 2
N 1 ω cos ω τ 1 N 2 sin ω τ 1 = M 1 ω .
By squaring and adding Equations (18) and (19), we get
ω 4 + ( M 1 2 N 1 2 2 M 2 ) ω 2 + M 2 2 N 2 2 = 0 .
By substituting ω 2 = σ , we get the following equation in σ :
Ω ( σ ) = σ 2 + ( M 1 2 N 1 2 2 M 2 ) σ + M 2 2 N 2 2 = 0 .
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:
τ 1 j = 1 ω 10 cos 1 ( N 2 M 1 N 1 ) ω 10 2 N 2 M 2 N 2 2 + N 1 2 ω 10 2 + 2 j π ω 10 ,
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 τ 1 < τ 1 0 , where τ 1 0 = min j τ 1 j , and E * is unstable for τ 1 τ 1 0 , provided the transversality condition is satisfied. We now verify the transversality condition d R e ( λ ) d τ 1 τ 1 = τ 1 0 0 . By differentiating Equation (17) with respect to τ 1 , we get:
d λ d τ 1 = λ ( N 1 λ + N 2 ) e λ τ 1 2 λ + M 1 + N 1 e λ τ 1 τ 1 ( N 1 λ + N 2 ) e λ τ 1 .
This gives
d λ d τ 1 1 = 2 λ + M 1 + N 1 e λ τ 1 λ ( N 1 λ + N 2 ) e λ τ 1 τ 1 λ .
Now,
sgn d R e ( λ ) d τ 1 τ 1 = τ 1 0 = sgn d R e ( λ ) d τ 1 τ 1 = τ 1 0 1 = sgn R e d λ d τ 1 1 λ = i ω 10 = sgn 2 ω 10 2 + ( M 1 2 N 1 2 2 M 2 ) N 1 2 ω 10 2 + N 2 2 = sgn Ω ( ω 10 2 ) N 1 2 ω 10 2 + N 2 2 .
Therefore, the transversality condition is satisfied, and so Hopf bifurcation occurs at τ 1 = τ 1 0 if Ω ( ω 10 2 ) 0 . The following theorem can now be stated.
Theorem 15.
Suppose that equilibrium E * exists and is LAS for τ 1 = τ 2 = 0 . Additionally, let σ * = ω 10 2 be a positive root of Equation (21). Then, τ 1 = τ 1 0 such that equilibrium E * of system (14) is LAS when 0 < τ 1 < τ 1 0 and unstable for τ 1 > τ 1 0 , where
τ 1 j = 1 ω 10 cos 1 ( N 2 M 1 N 1 ) ω 10 2 N 2 M 2 N 2 2 + N 1 2 ω 10 2 + 2 j π ω 10 ,
for j = 0 , 1 , 2 , . Furthermore, system (14) will exhibit a Hopf bifurcation around E * when τ 1 = τ 1 0 provided Ω ( ω 10 2 ) 0 .
  • Case III
τ 1 = 0 , τ 2 > 0 .
Here, we consider the case when breeding delay ( τ 1 ) is absent but gestation delay ( τ 2 ) is present. Therefore, the characteristic equation, Equation (16), becomes
λ 2 + a 1 λ + a 2 + ( b 1 λ + b 2 ) e λ τ 2 = 0 ,
where a 1 = A 1 + B 1 , a 2 = A 2 + B 2 , b 1 = C 1 and b 2 = C 2 + D 1 . We can state the following theorem by relying on the same analysis as in Case II.
Theorem 16.
Suppose that equilibrium E * exists and is LAS for τ 1 = τ 2 = 0 . Additionally, let ( A 2 + B 2 ) 2 < ( C 2 + D 1 ) 2 when τ 1 = 0 . Then, τ 2 = τ 2 0 such that equilibrium E * of system (14) is LAS when 0 < τ 2 < τ 2 0 and unstable for τ 2 > τ 2 0 , where
τ 2 j = 1 ω 20 cos 1 ( b 2 a 1 b 1 ) ω 20 2 b 2 a 2 b 2 2 + b 1 2 ω 20 2 + 2 j π ω 20 ,
for j = 0 , 1 , 2 , . Furthermore, system (14) will exhibit a Hopf bifurcation around E * when τ 2 = τ 2 0 , provided d R e ( λ ) d τ 2 τ 2 = τ 2 0 0 .
  • Case IV
τ 2 ( 0 , τ 2 0 ) , τ 1 > 0 .
Assume τ 2 as a fixed number in the interval ( 0 , τ 2 0 ) and τ 1 > 0 . Let λ = u + i v be the root of the characteristic equation, Equation (16). By putting this value into Equation (16), we have
u 2 + 2 u v i v 2 + A 1 u + A 1 v i + A 2 + { ( B 1 u + B 2 ) + B 1 v i } ( cos v τ 1 i sin v τ 1 ) e u τ 1 + { ( C 1 u + C 2 ) + i C 1 v } ( cos v τ 2 i sin v τ 2 ) e u τ 2 + D 1 cos ( τ 1 + τ 2 ) v i sin ( τ 1 + τ 2 ) v e ( τ 1 + τ 2 ) u = 0 ;
comparing real and imaginary parts,
u 2 v 2 + A 1 u + A 2 + { ( B 1 u + B 2 ) cos v τ 1 + B 1 v sin v τ 1 ) } e u τ 1 + { ( C 1 u + C 2 ) cos v τ 2 + C 1 v sin v τ 2 } e u τ 2 + D 1 e ( τ 1 + τ 2 ) u cos ( v ( τ 1 + τ 2 ) ) = 0
and
2 u v + A 1 v + B 1 v cos v τ 1 ( B 1 u + B 2 ) sin v τ 1 e u τ 1 + { C 1 v cos v τ 2 ( C 1 u + C 2 ) sin v τ 2 } e u τ 2 D 1 e ( τ 1 + τ 2 ) u sin ( v ( τ 1 + τ 2 ) ) = 0 .
The root of the characteristic equation, Equation (16), must be purely imaginary for equilibrium E * to change in stability. Therefore, by putting u = 0 in Equations (25) and (26), we get
( B 2 + D 1 cos v τ 2 ) cos v τ 1 + ( B 1 v D 1 sin v τ 2 ) sin v τ 1 v 2 + A 2 + C 2 cos v τ 2 + C 1 v sin v τ 2 = 0 ,
( B 1 v D 1 sin v τ 2 ) cos v τ 1 ( B 2 + D 1 cos v τ 2 ) sin v τ 1 + A 1 v + C 1 v cos v τ 2 C 2 sin v τ 2 = 0 .
By squaring and adding Equations (27) and (28), we have
v 4 + δ 1 v 3 + δ 2 v 2 + δ 3 v + δ 4 = 0 ,
δ 1 = 2 C 1 sin v τ 2 δ 2 = A 1 2 + C 1 2 B 1 2 2 A 2 + 2 ( A 1 C 1 C 2 ) cos v τ 2 δ 3 = 2 ( A 2 C 1 A 1 C 2 + B 1 D 1 ) sin v τ 2 δ 4 = A 2 2 B 2 2 + C 2 2 D 1 2 + 2 ( A 2 C 2 B 2 D 1 ) cos v τ 2 .
If δ 4 < 0 , Equation (29) has a positive root, say, v = ω 10 * . Consequently, Equation (16) has a pair of imaginary roots λ = ± ω 10 * for some fixed τ 2 ( 0 , τ 2 0 ) . Therefore, by solving Equations (27) and (28) for τ 1 at v = ω 10 * , we get:
τ 1 * j = 1 ω 10 * c o s 1 ( η 1 ω 10 * 2 η 2 ) + ( η 3 ω 10 * 2 η 4 ) cos ( ω 10 * τ 2 ) + η 5 ω 10 * sin ( ω 10 * τ 2 ) B 2 2 + D 1 2 + B 1 2 ω 10 * 2 + 2 B 2 D 1 cos ( ω 10 * τ 2 ) 2 B 1 D 1 ω 10 * sin ( ω 10 * τ 2 ) + 2 π j ω 10 *
for j = 0 , 1 , 2 , , where η 1 = ( B 2 + A 1 B 1 ) , η 2 = ( D 1 C 2 + A 2 B 2 ) , η 3 = ( D 1 C 1 B 1 ) , η 4 = ( A 2 D 1 + B 2 C 2 ) and η 5 = ( A 1 D 1 + B 1 C 2 B 2 C 1 ) .
Now, if equilibrium E * is LAS in absence of both the time delays, then by Butler’s lemma, it will remain stable for 0 < τ 1 < τ 1 * 0 , such that τ 1 * 0 = min j 0 τ 1 * j and system (14) undergoes Hopf bifurcation (around E * ) at τ 1 * 0 , provided the transversality condition, d R e ( λ ) d τ 1 τ 1 = τ 1 * 0 0 , is satisfied.
By differentiating Equations (25) and (26) with respect to τ 1 , we get:
P 1 d R e ( λ ) d τ 1 τ 1 = τ 1 * 0 + P 2 d I m ( λ ) d τ 1 τ 1 = τ 1 * 0 = R 1
P 2 d R e ( λ ) d τ 1 τ 1 = τ 1 * 0 + P 1 d I m ( λ ) d τ 1 τ 1 = τ 1 * 0 = R 2
where
P 1 = A 1 + ( B 1 τ 1 * 0 B 2 ) cos ( ω 10 * τ 1 * 0 ) τ 1 * 0 B 1 ω 10 * sin ( ω 10 * τ 1 * 0 ) + ( C 1 C 2 τ 2 ) cos ( ω 10 * τ 2 ) C 1 ω 10 * τ 2 sin ( ω 10 * τ 2 ) D 1 ( τ 1 * 0 + τ 2 ) cos ( ω 10 * ( τ 1 * 0 + τ 2 ) ) , P 2 = 2 ω 10 * + ( B 1 τ 1 * 0 B 2 ) sin ( ω 10 * τ 1 * 0 ) + τ 1 * 0 B 1 ω 10 * cos ( ω 10 * τ 1 * 0 ) + ( C 1 τ 2 C 2 ) sin ( ω 10 * τ 2 ) + C 1 ω 10 * τ 2 cos ( ω 10 * τ 2 ) D 1 ( τ 1 * 0 + τ 2 ) sin ω 10 * ( τ 1 * 0 + τ 2 ) , R 1 = B 2 ω 10 * sin ( ω 10 * τ 1 * 0 ) + B 1 ω 10 * 2 cos ( ω 10 * τ 1 * 0 ) D 1 ω 10 * sin ( ω 10 * ( τ 1 * 0 + τ 2 ) ) , R 2 = B 1 ω 10 * 2 sin ( ω 10 * τ 1 * 0 ) B 2 ω 10 * cos ( ω 10 * τ 1 * 0 ) D 1 ω 10 * cos ( ω 10 * ( τ 1 * 0 + τ 2 ) ) .
By solving Equations (31) and (32), we get:
d R e ( λ ) d τ 1 τ 1 = τ 1 * 0 = P 1 R 1 P 2 R 2 P 1 2 + P 2 2 .
Therefore, Hopf bifurcation (around E * ) occurs at τ 1 * 0 if P 1 R 1 P 2 R 2 0 .
Theorem 17.
If equilibrium E * exists and is LAS in absence of both the time delays, then for τ 2 ( 0 , τ 2 0 ) if δ 4 < 0 , there exists τ 1 = τ 1 * 0 such that equilibrium E * is LAS when 0 < τ 1 < τ 1 * 0 and unstable when τ 1 > τ 1 * 0 , where
τ 1 * j = 1 ω 10 * c o s 1 ( η 1 ω 10 * 2 η 2 ) + ( η 3 ω 10 * 2 η 4 ) cos ( ω 10 * τ 2 ) + η 5 ω 10 * sin ( ω 10 * τ 2 ) B 2 2 + D 1 2 + B 1 2 ω 10 * 2 + 2 B 2 D 1 cos ( ω 10 * τ 2 ) 2 B 1 D 1 ω 10 * sin ( ω 10 * τ 2 ) + 2 π j ω 10 *
for j = 0 , 1 , 2 , . Further, system (14) will exhibit a Hopf bifurcation around E * for τ 1 = τ 1 * 0 , provided P 1 R 1 P 2 R 2 0 .
  • Case V
τ 1 ( 0 , τ 1 0 ) , τ 2 > 0 .
We can state the following theorem relying on the same analysis as in Case IV.
Theorem 18.
If equilibrium E * exists and is LAS in absence of both the time delays, then for τ 1 ( 0 , τ 1 0 ) if δ ¯ 4 < 0 , there exists τ 2 = τ 2 * 0 such that equilibrium E * is LAS when 0 < τ 2 < τ 2 * 0 and unstable when τ 2 > τ 2 * 0 , where
τ 2 * j = 1 ω 20 * c o s 1 ( η ¯ 1 ω 20 * 2 η ¯ 2 ) + ( η ¯ 3 ω 20 * 2 η ¯ 4 ) cos ( ω 20 * τ 2 ) + η ¯ 5 ω 20 * sin ( ω 20 * τ 2 ) C 2 2 + D 1 2 + C 1 2 ω 20 * 2 + 2 C 2 D 1 cos ( ω 20 * τ 2 ) 2 C 1 D 1 ω 20 * sin ( ω 20 * τ 2 ) + 2 π j ω 20 *
for j = 0 , 1 , 2 , . Further, system (14) will exhibit a Hopf bifurcation around equilibrium E * for τ 2 = τ 2 * 0 , provided P ¯ 1 R ¯ 1 P ¯ 2 R ¯ 2 0 where
η ¯ 1 = ( C 2 + A 1 C 1 ) , η ¯ 2 = ( D 1 B 2 + A 2 C 2 ) , η ¯ 3 = ( D 1 B 1 C 1 ) , η ¯ 4 = ( A 2 D 1 + B 2 C 2 ) , η ¯ 5 = ( A 1 D 1 B 1 C 2 + B 2 C 1 ) , P ¯ 1 = A 1 + ( C 1 τ 2 * 0 C 2 ) cos ( ω 20 * τ 2 * 0 ) τ 2 * 0 C 1 ω 20 * sin ( ω 20 * τ 2 * 0 ) + ( B 1 B 2 τ 1 ) cos ( ω 20 * τ 1 ) B 1 ω 20 * τ 1 sin ( ω 20 * τ 1 ) D 1 ( τ 2 * 0 + τ 1 ) cos ( ω 20 * ( τ 2 * 0 + τ 1 ) ) , P ¯ 2 = 2 ω 20 * + ( C 1 τ 2 * 0 C 2 ) sin ( ω 20 * τ 2 * 0 ) + τ 2 * 0 C 1 ω 20 * cos ( ω 20 * τ 2 * 0 ) + ( B 1 τ 1 B 2 ) sin ( ω 20 * τ 1 ) + B 1 ω 20 * τ 1 cos ( ω 20 * τ 1 ) D 1 ( τ 2 * 0 + τ 1 ) sin ω 20 * ( τ 2 * 0 + τ 1 ) , R ¯ 1 = C 2 ω 20 * sin ( ω 20 * τ 2 * 0 ) + C 1 ω 20 * 2 cos ( ω 20 * τ 2 * 0 ) D 1 ω 20 * sin ( ω 20 * ( τ 2 * 0 + τ 1 ) ) , R ¯ 2 = C 1 ω 20 * 2 sin ( ω 20 * τ 2 * 0 ) C 2 ω 20 * cos ( ω 20 * τ 2 * 0 ) D 1 ω 20 * cos ( ω 20 * ( τ 2 * 0 + τ 1 ) ) , δ ¯ 4 = A 2 2 + B 2 2 C 2 2 D 1 2 + 2 ( A 2 B 2 C 2 D 1 ) cos v τ 1 .

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: { b = 5.5 , d 1 = 5.6 , d 2 = 0.2 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.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 1 [ T C 2 ] = 5.5 , i.e., when d 1 = b , one eigenvalue is zero, and the stability of trivial equilibrium point, E 0 ( 0 , 0 ) , exchanges with that of the axial equilibrium point, E 1 , i.e., transcritical bifurcation occurs. If we further decrease d 1 , at d 1 [ T C 1 ] = 5.444 , the axial equilibrium point exchanges its stability with that of the positive interior equilibrium point E * , and finally at d 1 [ H 0 ] = 2.952 , 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: { b = 5.5 , d 1 = 0.1 , d 2 = 0.95 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.2 , h = 1 and γ = 0.8 } . Here, b > d 1 and the stability condition γ α b d 1 β 1 + α h b d 1 β < d 2 is satisfied. Therefore, E 1 b d 1 β , 0 E 1 ( 27 , 0 ) 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 b [ T C ] = 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 1 [ T C ] = 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 [ T C ] = 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 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 0.5 , k 2 = 0.01 , β = 0.6 , α = 2.5 , h = 1.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 1 [ H 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 α [ T C ] = 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.
Considering a parametric set as { b = 5.5 , d 1 = 0.1 , d 2 = 0.2 , k 1 = 2 , β = 0.2 , α = 1.2 h = 3.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.

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: { b = 8.7 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 1.5 , k 2 = 0.1 , β = 0.2 , α = 2 , h = 1 , γ = 0.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).
For the delay (Case II), when τ 2 = 0 and τ 1 > 0 , we have analytically found that delayed system, (14), exhibits a Hopf bifurcation at some critical value of τ 1 = τ 1 0 = 2.125 . 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 τ 1 = 1.5 < τ 1 0 , 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 E * ( x * , y * ) E * ( 2.5 , 1.564 ) for τ 1 = 2.5 > τ 1 0 . Therefore, the delay parameter τ 1 can make system (14) go from stable to unstable as we increase the value of τ 1 from zero onward while taking τ 2 = 0 .
For Case III: when τ 1 = 0 and τ 2 > 0 , we found that delayed system (14) undergoes a Hopf bifurcation for some critical value of τ 2 = τ 2 0 = 0.813 , as depicted in Figure 20. Again, we took τ 2 = 0.5 < τ 2 0 = 0.813 and τ 2 = 2.1 > τ 2 0 = 0.813 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 τ 1 when τ 2 = 0.2 ( 0 , τ 2 0 = 0.813 ) and Hopf bifurcation occurs at τ 1 * 0 = 1.51 . Additionally, we have drawn the time series of the delayed system, (14), with τ 1 = 1.2 < τ 1 * 0 (see Figure 24), τ 1 = 2 > τ 1 * 0 (see Figure 25) and fixed τ 2 = 0.2 .
It can be observed in Figure 24 and Figure 25 that when 0 < τ 2 < τ 2 0 , the system is stable for τ 1 < τ 1 * 0 and unstable otherwise. However, numerically it can be noticed that when τ 2 τ 2 0 , for any positive value of τ 1 , the system is always unstable and forms a limit cycle around the coexistence state E * ( 2.5 , 1.564 ) (see Figure 26). As a consequence, system (14) does not exhibit any switching behavior for any τ 2 τ 2 0 , regardless of the value of τ 1 .
Finally, in Case V: for τ 1 = 1 ( 0 , τ 1 0 ) = ( 0 , 2.125 ) and τ 2 > 0 , we have drawn the bifurcation diagram with τ 2 . It can be observed that system (14) exhibits a Hopf bifurcation at some critical point τ 2 * 0 = 0.456 (see Figure 27). We have drawn the time series for system (14) with τ 2 = 0.25 < τ 2 * 0 = 0.456 (see Figure 28) and τ 2 = 1 > τ 2 * 0 = 0.466 (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 τ 2 < τ 2 * 0 provided 0 < τ 1 < τ 1 0 , and unstable otherwise. However, it is evident from numerical simulations that an unstable limit cycle is generated when τ 1 τ 1 0 for any positive value of τ 2 (see Figure 30). This leads to the conclusion that system (14) has no switching behavior for any τ 1 τ 1 0 , irrespective of the values of τ 2 .

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 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:
γ i ( t ) = 0 and γ i ( t 1 ) γ i ( t 2 ) = μ i 2 δ i ( t 1 t 2 ) for i = 1 , 2 .
Here, μ i > 0 is the intensity or strength of γ i , and the Dirac delta function δ is defined as follows:
{ δ ( x ) = 0 , for x 0 , δ ( x ) d x = 1
where · is the ensemble average of the stochastic process under consideration. Thus, system (3) was modified as follows:
d x d t = b x 1 + k 1 y ( d 1 + γ 1 ( t ) ) x β ( 1 + k 2 y ) x 2 α x y 1 + α h x , d y d t = γ α x y 1 + α h x ( d 2 + γ 2 ( t ) ) y ,
i . e . , d x d t = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x μ 1 x d w 1 d t , d y d t = γ α x y 1 + α h x d 2 y μ 2 y d w 2 d t .
We took γ 1 = μ 1 d w 1 d t and γ 2 = μ 2 d w 2 d t , 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:
d x = b x 1 + k 1 y d 1 x β ( 1 + k 2 y ) x 2 α x y 1 + α h x d t μ 1 x d w 1 , d y = γ α x y 1 + α h x d 2 y d t μ 2 y d w 2 .
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 3.169 , 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 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.

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 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 [ t c ] = γ α b d 1 β 1 + α h b d 1 β , the system experiences transcritical bifurcation. The predator-free equilibrium point exchanges its stability with an interior equilibrium point at this threshold value, d 2 [ t c ] . The effect of the death rate of the predators, d 2 , has been numerically studied, as shown in Figure 7a. Whenever the value of death rate d 2 exceeds its threshold value, d 2 [ t c ] , 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 2 [ t c ] , the coexistence equilibrium point is stable up to some threshold value of d 2 [ H ] . The prey and predator populations oscillate periodically if d 2 diminishes further. Again, Theorem 12 states that at γ [ t c ] = d 2 1 + β α ( b d 1 ) , 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, d y * d k 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, d y * d k 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 [ T C ] < 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) oscillate around the deterministic equilibrium values of 2 and 3.169 , 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 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 Figure 2 and Figure 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 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.

Author Contributions

All the authors have participate equally in all the aspects of this paper: conceptualization, methodology, investigation, formal analysis, writing—original draft preparation, writing—review and editing. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that were utilized to support the outcomes of this investigation are contained in the article.

Acknowledgments

The authors are grateful to the anonymous referees and the learned editor, for their careful reading, valuable comments and helpful suggestions, which have helped them to improve the presentation of this work significantly.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Holmes, E.E.; Lewis, M.A.; Banks, J.E.; Veit, R.R. Partial Differential Equations in Ecology: Spatial Interactions and Population Dynamics. Ecology 1994, 75, 17–29. [Google Scholar] [CrossRef]
  2. Murray, J.D. Mathematical Biology II. Spatial Models and Biological Applications; Springer: New York, NY, USA, 2003. [Google Scholar]
  3. Murray, J.D. Mathematical Biology I. An Introduction; Springer: New York, NY, USA, 2002. [Google Scholar]
  4. Shigesada, N.; Kawasaki, K.; Teramoto, E. Spatial segregation of interacting species. J. Theor. Biol. 1979, 79, 83–99. [Google Scholar] [CrossRef]
  5. Wang, W.; Zhang, Y.; Liu, C. Analysis of a discrete-time predator–prey system with Allee effect. Ecol. Complex. 2011, 8, 81–85. [Google Scholar] [CrossRef]
  6. Yin, C.; Chen, Y.; Zhong, S. Fractional-order sliding mode based extremum seeking control of a class of nonlinear systems. Automatica 2014, 50, 3173–3181. [Google Scholar] [CrossRef]
  7. Perc, M.; Szolnoki, A.; Szabó, G. Cyclical interactions with alliance-specific heterogeneous invasion rates. Phys. Rev. E 2007, 75, 052102. [Google Scholar] [CrossRef] [Green Version]
  8. Zhang, T.; Zang, H. Delay-induced Turing instability in reaction-diffusion equations. Phys. Rev. E 2014, 90, 052908. [Google Scholar] [CrossRef]
  9. Gross, T.; Blasius, B. Adaptive coevolutionary networks: A review. J. R. Soc. Interface 2008, 5, 259–271. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Malthus, T.R. An Essay on the Principle of Population, as It Affects the Future Imporvement of Society, with Remarks on the Speculations of Mr. Godwin, M. Condorcet, and Other Writers; The Lawbook Exchange, Ltd.: Clark, NJ, USA, 2007. [Google Scholar]
  11. Verhulst, P. Notice sur la loi que la population suit dans son accroissement. Corresp. Math. Phys. 1838, 10, 113–126. [Google Scholar]
  12. Lotka, A.J. Elements of Physical Biology; Williams & Wilkins: Philadelphia, PA, USA, 1925. [Google Scholar]
  13. Mondal, S.; Samanta, G. Time-delayed predator-prey interaction with the benefit of antipredation response in presence of refuge. Z. Naturforschung A 2021, 76, 23–42. [Google Scholar] [CrossRef]
  14. Sahoo, D.; Samanta, G.P. Impact of fear effect in a two prey-one predator system with switching behaviour in predation. Differ. Equ. Dyn. Syst. 2021, 1–23. [Google Scholar] [CrossRef]
  15. Samanta, G.; Mondal, A.; Sahoo, D.; Dolai, P. A prey-predator system with herd behaviour of prey in a rapidly fluctuating environment. Math. Appl. Sci. Eng. 2020, 1, 16–26. [Google Scholar] [CrossRef]
  16. Creel, S.; Christianson, D. Relationships between direct predation and risk effects. Trends Ecol. Evol. 2008, 23, 194–201. [Google Scholar] [CrossRef] [PubMed]
  17. Lima, S.L. Predators and the breeding bird: Behavioral and reproductive flexibility under the risk of predation. Biol. Rev. 2009, 84, 485–513. [Google Scholar] [CrossRef]
  18. Wang, X.; Zanette, L.; Zou, X. Modelling the fear effect in predator–prey interactions. J. Math. Biol. 2016, 73, 1179–1204. [Google Scholar] [CrossRef] [PubMed]
  19. Cresswell, W. Predation in bird populations. J. Ornithol. 2011, 152, 251–263. [Google Scholar] [CrossRef]
  20. Pettorelli, N.; Coulson, T.; Durant, S.M.; Gaillard, J. Predation, individual variability and vertebrate population dynamics. Oecologia 2011, 167, 305–314. [Google Scholar] [CrossRef]
  21. Hua, F.; Sieving, K.E.; Fletcher, R.J., Jr.; Wright, C.A. Increased perception of predation risk to adults and offspring alters avian reproductive strategy and performance. Behav. Ecol. 2014, 25, 509–519. [Google Scholar] [CrossRef]
  22. Altendorf, K.B.; Laundré, J.W.; López González, C.A.; Brown, J.S. Assessing effects of predation risk on foraging behavior of mule deer. J. Mammal. 2001, 82, 430–439. [Google Scholar] [CrossRef]
  23. Zanette, L.Y.; White, A.F.; Allen, M.C.; Clinchy, M. Perceived predation risk reduces the number of offspring songbirds produce per year. Science 2011, 334, 1398–1401. [Google Scholar] [CrossRef]
  24. Creel, S.; Christianson, D.; Liley, S.; Winnie, J.A., Jr. Predation risk affects reproductive physiology and demography of elk. Science 2007, 315, 960. [Google Scholar] [CrossRef] [Green Version]
  25. Zhang, H.; Cai, Y.; Fu, S.; Wang, W. Impact of the fear effect in a prey-predator model incorporating a prey refuge. Appl. Math. Comput. 2019, 356, 328–337. [Google Scholar] [CrossRef]
  26. Wang, Y.; Zou, X. On a predator–prey system with digestion delay and anti-predation strategy. J. Nonlinear Sci. 2020, 30, 1579–1605. [Google Scholar] [CrossRef]
  27. Das, B.K.; Sahoo, D.; Samanta, G.P. Impact of fear in a delay-induced predator–prey system with intraspecific competition within predator species. Math. Comput. Simul. 2022, 191, 134–156. [Google Scholar] [CrossRef]
  28. Kuang, Y. Delay Differential Equations: With Applications in Population Dynamics; Academic Press: Cambridge, MA, USA, 1993. [Google Scholar]
  29. MacDonald, N.; MacDonald, N. Biological Delay Systems: Linear Stability Theory; Cambridge University Press: Cambridge, UK, 2008. [Google Scholar]
  30. Ruan, S. On nonlinear dynamics of predator-prey models with discrete delay. Math. Model. Nat. Phenom. 2009, 4, 140–188. [Google Scholar] [CrossRef]
  31. Gurney, W.S.C.; Nisbet, R.M. Fluctuation periodicity, generation separation, and the expression of larval competition. Theor. Popul. Biol. 1985, 28, 150–180. [Google Scholar] [CrossRef]
  32. Cooke, K.L.; van den Driessche, P.; Zou, X. Interaction of maturation delay and nonlinear birth in population and epidemic models. J. Math. Biol. 1999, 39, 332–352. [Google Scholar] [CrossRef]
  33. Fan, G.; Liu, J.; van den Driessche, P.; Wu, J.; Zhu, H. The impact of maturation delay of mosquitoes on the transmission of West Nile virus. Math. Biosci. 2010, 228, 119–126. [Google Scholar] [CrossRef]
  34. Holling, C.S. The components of predation as revealed by a study of small-mammal predation of the european pine sawfly. Can. Entomol. 1959, 91, 293–320. [Google Scholar] [CrossRef]
  35. Holling, C.S. Some Characteristics of Simple Types of Predation and Parasitism. Can. Entomol. 1959, 91, 385–398. [Google Scholar] [CrossRef]
  36. Holling, C.S. The Functional Response of Predators to Prey Density and its Role in Mimicry and Population Regulation. Mem. Entomol. Soc. Can. 1965, 97, 5–60. [Google Scholar] [CrossRef] [Green Version]
  37. Dutta, P.; Sahoo, D.; Mondal, S.; Samanta, G. Dynamical complexity of a delay-induced eco-epidemic model with Beddington–DeAngelis incidence rate. Math. Comput. Simul. 2022, 197, 45–90. [Google Scholar] [CrossRef]
  38. Mondal, S.; Samanta, G.P. Dynamics of a delayed predator–prey interaction incorporating nonlinear prey refuge under the influence of fear effect and additional food. J. Phys. A Math. Theor. 2020, 53, 295601. [Google Scholar] [CrossRef]
  39. Sarkar, K.; Khajanchi, S. Impact of fear effect on the growth of prey in a predator-prey interaction model. Ecol. Complex. 2020, 42, 100826. [Google Scholar] [CrossRef]
  40. Halder, S.; Bhattacharyya, J.; Pal, S. Predator-prey interactions under fear effect and multiple foraging strategies. Discret. Contin. Dyn. Syst. B 2022, 27, 3779. [Google Scholar] [CrossRef]
  41. Freedman, H.I.; Ruan, S. Uniform persistence in functional differential equations. J. Differ. Equ. 1995, 115, 173–192. [Google Scholar] [CrossRef] [Green Version]
  42. Wei, J.; Ruan, S. Stability and bifurcation in a neural network model with two delays. Phys. D Nonlinear Phenom. 1999, 130, 255–272. [Google Scholar] [CrossRef]
  43. Ruan, S.; Wei, J. On the zeros of transcendental functions with applications to stability of delay differential equations with two delays. Dyn. Contin. Discret. Impuls. Syst. Ser. A 2003, 10, 863–874. [Google Scholar]
  44. Mondal, S.; Maiti, A.; Samanta, G.P. Effects of fear and additional food in a delayed predator–prey model. Biophys. Rev. Lett. 2018, 13, 157–177. [Google Scholar] [CrossRef]
  45. Das, A.; Samanta, G.P. Modeling the fear effect on a stochastic prey–predator system with additional food for the predator. J. Phys. A Math. Theor. 2018, 51, 465601. [Google Scholar] [CrossRef]
  46. Mondal, S.; Samanta, G. A comparison study of predator–prey system in deterministic and stochastic environments influenced by fear and its carry-over effects. Eur. Phys. J. Plus 2022, 137, 70. [Google Scholar] [CrossRef]
Figure 1. Graphical representation of predator–prey nullclines for the parametric set: { b = 10 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 0.5 , k 2 = 0.01 , β = 0.6 , α = 2.5 , h = 1.2 , γ = 0.7 } .
Figure 1. Graphical representation of predator–prey nullclines for the parametric set: { b = 10 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 0.5 , k 2 = 0.01 , β = 0.6 , α = 2.5 , h = 1.2 , γ = 0.7 } .
Mathematics 10 03795 g001
Figure 2. Stability nature of E 0 ( 0 , 0 ) under the parametric set: { b = 5.5 , d 1 = 5.6 , d 2 = 0.2 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.2 , h = 1 , γ = 0.8 } .
Figure 2. Stability nature of E 0 ( 0 , 0 ) under the parametric set: { b = 5.5 , d 1 = 5.6 , d 2 = 0.2 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.2 , h = 1 , γ = 0.8 } .
Mathematics 10 03795 g002
Figure 3. Bifurcation diagram regarding parameter d 1 , and all other parameters are the same as in Figure 2.
Figure 3. Bifurcation diagram regarding parameter d 1 , and all other parameters are the same as in Figure 2.
Mathematics 10 03795 g003
Figure 4. Stability nature of the predator-free steady state, E 1 b d 1 β , 0 , under the parametric set: { b = 5.5 , d 1 = 0.1 , d 2 = 0.95 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.2 , h = 1 , γ = 0.8 } .
Figure 4. Stability nature of the predator-free steady state, E 1 b d 1 β , 0 , under the parametric set: { b = 5.5 , d 1 = 0.1 , d 2 = 0.95 , k 1 = 2 , k 2 = 1 , β = 0.2 , α = 1.2 , h = 1 , γ = 0.8 } .
Mathematics 10 03795 g004
Figure 5. Bifurcation diagram regarding parameters (a) b and (b) d 1 . All other parameters are the same as in Figure 4.
Figure 5. Bifurcation diagram regarding parameters (a) b and (b) d 1 . All other parameters are the same as in Figure 4.
Mathematics 10 03795 g005
Figure 6. Bifurcation diagram regarding parameter h, and all other parameters are the same as in Figure 4.
Figure 6. Bifurcation diagram regarding parameter h, and all other parameters are the same as in Figure 4.
Mathematics 10 03795 g006
Figure 7. Bifurcation diagrams regarding parameters (a) d 2 and (b) γ . All other parameters are the same as in Figure 4.
Figure 7. Bifurcation diagrams regarding parameters (a) d 2 and (b) γ . All other parameters are the same as in Figure 4.
Mathematics 10 03795 g007
Figure 8. Stability nature of the coexistence steady state E * ( x * , y * ) E * ( 2 , 3.169 ) under the parametric set: { b = 10 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 0.5 , k 2 = 0.01 , β = 0.6 , α = 2.5 , h = 1.2 , γ = 0.7 } . (a) Time series; (b) Phase portrait.
Figure 8. Stability nature of the coexistence steady state E * ( x * , y * ) E * ( 2 , 3.169 ) under the parametric set: { b = 10 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 0.5 , k 2 = 0.01 , β = 0.6 , α = 2.5 , h = 1.2 , γ = 0.7 } . (a) Time series; (b) Phase portrait.
Mathematics 10 03795 g008
Figure 9. Bifurcation diagram regarding parameter k 1 , the other parameters are the same as those in Figure 8. (a) Prey population; (b) Predator population.
Figure 9. Bifurcation diagram regarding parameter k 1 , the other parameters are the same as those in Figure 8. (a) Prey population; (b) Predator population.
Mathematics 10 03795 g009
Figure 10. Bifurcation diagram regarding parameter β when the other parameters are the same as those in Figure 8. (a) Prey population; (b) Predator population.
Figure 10. Bifurcation diagram regarding parameter β when the other parameters are the same as those in Figure 8. (a) Prey population; (b) Predator population.
Mathematics 10 03795 g010
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 11. Bifurcation diagram regarding parameter α when the other parameters are the same as those in Figure 8. (a) Prey population; (b) Predator population.
Mathematics 10 03795 g011
Figure 12. Nature of the growth of predator population with respect to k 2 when the remaining parameters are the same as those in Figure 8.
Figure 12. Nature of the growth of predator population with respect to k 2 when the remaining parameters are the same as those in Figure 8.
Mathematics 10 03795 g012
Figure 13. Bifurcation diagram of prey population regarding parameter b when the other parameters are the same as those in Figure 8. (a) k 1 = 0 , k 2 = 0 ; (b) k 1 = 0.5 , k 2 = 0 ; (c) k 1 = 0.5 , k 2 = 0.01 .
Figure 13. Bifurcation diagram of prey population regarding parameter b when the other parameters are the same as those in Figure 8. (a) k 1 = 0 , k 2 = 0 ; (b) k 1 = 0.5 , k 2 = 0 ; (c) k 1 = 0.5 , k 2 = 0.01 .
Mathematics 10 03795 g013
Figure 14. Bifurcation diagram regarding parameter k 2 when taking the parameter set: { b = 5.5 , d 1 = 0.1 , d 2 = 0.2 , k 1 = 2 , β = 0.2 , α = 1.2 h = 3.4 γ = 0.8 } . (a) Prey population; (b) Predator population.
Figure 14. Bifurcation diagram regarding parameter k 2 when taking the parameter set: { b = 5.5 , d 1 = 0.1 , d 2 = 0.2 , k 1 = 2 , β = 0.2 , α = 1.2 h = 3.4 γ = 0.8 } . (a) Prey population; (b) Predator population.
Mathematics 10 03795 g014
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 15. Two parametric bifurcation diagrams in the ( k 2 k 1 ) parametric plane with all other parameters the same as in Figure 14.
Mathematics 10 03795 g015
Figure 16. Stability nature of the coexistence equilibrium point E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = τ 2 = 0 and b = 8.7 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 1.5 , k 2 = 0.1 , β = 0.2 , α = 2 , h = 1 and γ = 0.6 .
Figure 16. Stability nature of the coexistence equilibrium point E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = τ 2 = 0 and b = 8.7 , d 1 = 1.5 , d 2 = 0.5 , k 1 = 1.5 , k 2 = 0.1 , β = 0.2 , α = 2 , h = 1 and γ = 0.6 .
Mathematics 10 03795 g016
Figure 17. Bifurcation diagram regarding delay parameter τ 1 when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Figure 17. Bifurcation diagram regarding delay parameter τ 1 when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g017
Figure 18. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1.5 < τ 1 0 = 2.125 when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Figure 18. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1.5 < τ 1 0 = 2.125 when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g018
Figure 19. Oscillatory nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 2.5 > τ 1 0 = 2.125 , when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Figure 19. Oscillatory nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 2.5 > τ 1 0 = 2.125 , when τ 2 = 0 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g019
Figure 20. Bifurcation diagram regarding delay parameter τ 2 when τ 1 = 0 and all other parameter values are the same as in Figure 16.
Figure 20. Bifurcation diagram regarding delay parameter τ 2 when τ 1 = 0 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g020
Figure 21. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 0 , τ 2 = 0.5 < τ 2 0 = 0.813 and all other parameter values are the same as in Figure 16.
Figure 21. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 0 , τ 2 = 0.5 < τ 2 0 = 0.813 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g021
Figure 22. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 0 , τ 2 = 2.1 > τ 2 0 = 0.813 and all other parameter values are the same as in Figure 16.
Figure 22. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 0 , τ 2 = 2.1 > τ 2 0 = 0.813 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g022
Figure 23. Bifurcation diagram regarding delay parameter τ 1 when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) and all other parameter values are the same as in Figure 16.
Figure 23. Bifurcation diagram regarding delay parameter τ 1 when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g023
Figure 24. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) , τ 1 = 1.2 < τ 1 * 0 = 1.51 and all other parameter values are the same as in Figure 16.
Figure 24. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) , τ 1 = 1.2 < τ 1 * 0 = 1.51 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g024
Figure 25. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) , τ 1 = 2 > τ 1 * 0 = 1.51 and all other parameter values are the same as in Figure 16.
Figure 25. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 2 = 0.2 ( 0 , τ 2 0 ) = ( 0 , 0.813 ) , τ 1 = 2 > τ 1 * 0 = 1.51 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g025
Figure 26. (a) Oscillatory nature with time t and (b) limit cycle around E * ( 2.5 , 1.564 ) when τ 2 τ 2 0 and τ 1 > 0 . All parameters are the same as in Figure 16.
Figure 26. (a) Oscillatory nature with time t and (b) limit cycle around E * ( 2.5 , 1.564 ) when τ 2 τ 2 0 and τ 1 > 0 . All parameters are the same as in Figure 16.
Mathematics 10 03795 g026
Figure 27. Bifurcation diagram regarding delay parameter τ 2 when τ 1 = 1 and all other parameter values are the same as in Figure 16.
Figure 27. Bifurcation diagram regarding delay parameter τ 2 when τ 1 = 1 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g027
Figure 28. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1 ( 0 , τ 1 0 ) = ( 0 , 2.125 ) , τ 2 = 0.25 < τ 2 * 0 = 0.456 and all other parameter values are the same as in Figure 16.
Figure 28. Stability nature of E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1 ( 0 , τ 1 0 ) = ( 0 , 2.125 ) , τ 2 = 0.25 < τ 2 * 0 = 0.456 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g028
Figure 29. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1 ( 0 , τ 1 0 ) = ( 0 , 2.125 ) , τ 2 = 0.5 > τ 2 * 0 = 0.456 and all other parameter values are the same as in Figure 16.
Figure 29. Oscillatory nature of the delayed system, (14), around E * ( x * , y * ) E * ( 2.5 , 1.564 ) when τ 1 = 1 ( 0 , τ 1 0 ) = ( 0 , 2.125 ) , τ 2 = 0.5 > τ 2 * 0 = 0.456 and all other parameter values are the same as in Figure 16.
Mathematics 10 03795 g029
Figure 30. (a) Oscillatory nature with time t and (b) limit cycle around E * ( 2.5 , 1.564 ) when τ 1 τ 1 0 and τ 2 > 0 . All parameters are same as in Figure 16.
Figure 30. (a) Oscillatory nature with time t and (b) limit cycle around E * ( 2.5 , 1.564 ) when τ 1 τ 1 0 and τ 2 > 0 . All parameters are same as in Figure 16.
Mathematics 10 03795 g030
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 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.
Mathematics 10 03795 g031
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 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.
Mathematics 10 03795 g032
Figure 33. Histograms of prey and predator populations corresponding to system (34).
Figure 33. Histograms of prey and predator populations corresponding to system (34).
Mathematics 10 03795 g033
Figure 34. Extinction of both species when μ 1 = μ 2 = 0.01 and all other parameter values are the same as in Figure 2.
Figure 34. Extinction of both species when μ 1 = μ 2 = 0.01 and all other parameter values are the same as in Figure 2.
Mathematics 10 03795 g034
Figure 35. Extinction of the predators when μ 1 = μ 2 = 0.01 and all other parameter values are the same as in Figure 4.
Figure 35. Extinction of the predators when μ 1 = μ 2 = 0.01 and all other parameter values are the same as in Figure 4.
Mathematics 10 03795 g035
Table 1. List of model parameters with units and their biological significance [37].
Table 1. List of model parameters with units and their biological significance [37].
ParameterDescriptionUnit
bintrinsic birth rate of the preyTime −1
k i , i = 1 , 2 level of fear due to prey’s anti-predation behaviorMass −1
d 1 natural mortality rate of the preyTime −1
β decay rate due to intra-specific competition of preyMass−1 Time−1
α attack rate by predatorMass −1Time−1
hprey handling timeTime
γ conversion efficiency of predatordimensionless
d 2 natural mortality rate of predatorTime −1
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Santra, N.; Mondal, S.; Samanta, G. Complex Dynamics of a Predator–Prey Interaction with Fear Effect in Deterministic and Fluctuating Environments. Mathematics 2022, 10, 3795. https://doi.org/10.3390/math10203795

AMA Style

Santra N, Mondal S, Samanta G. Complex Dynamics of a Predator–Prey Interaction with Fear Effect in Deterministic and Fluctuating Environments. Mathematics. 2022; 10(20):3795. https://doi.org/10.3390/math10203795

Chicago/Turabian Style

Santra, Nirapada, Sudeshna Mondal, and Guruprasad Samanta. 2022. "Complex Dynamics of a Predator–Prey Interaction with Fear Effect in Deterministic and Fluctuating Environments" Mathematics 10, no. 20: 3795. https://doi.org/10.3390/math10203795

APA Style

Santra, N., Mondal, S., & Samanta, G. (2022). Complex Dynamics of a Predator–Prey Interaction with Fear Effect in Deterministic and Fluctuating Environments. Mathematics, 10(20), 3795. https://doi.org/10.3390/math10203795

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop