Next Article in Journal
EEG-Based Person Identification and Authentication Using Deep Convolutional Neural Network
Next Article in Special Issue
Dynamics in an n-Species Lotka–Volterra Cooperative System with Delays
Previous Article in Journal
Functional Epistemology “Nullifies” Dyson’s Rebuttal of Perturbation Theory
Previous Article in Special Issue
Stability Analysis of Fractional-Order Predator-Prey System with Consuming Food Resource
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Dynamic Behavior of a Predator–Prey Model with Double Delays and Beddington–DeAngelis Functional Response

1
College of Science, Guilin University of Technology, Guilin 541004, China
2
School of Mathematical Sciences, Inner Mongolia University, Hohhot 010021, China
*
Author to whom correspondence should be addressed.
Axioms 2023, 12(1), 73; https://doi.org/10.3390/axioms12010073
Submission received: 25 October 2022 / Revised: 9 December 2022 / Accepted: 6 January 2023 / Published: 11 January 2023
(This article belongs to the Special Issue Impulsive, Delay and Fractional Order Systems)

Abstract

:
In the predator–prey system, predators can affect the prey population by direct killing and predation fear. In the present study, we consider a delayed predator–prey model with fear and Beddington–DeAngelis functional response. The model incorporates not only the fear of predator on prey with an intraspecific competition relationship, but also fear delay and pregnancy delay. Apart from the local stability analysis of the equilibrium points of the model, we find that time delay can change the stability of the system and cause Hopf bifurcation. Taking time delay as the bifurcation parameter, the critical values of delays in several cases are derived. In addition, we extend it to the random environment and study the stochastic ultimate boundedness of the stochastic process. Finally, our theoretical results are validated by numerical simulation.
MSC:
34C23; 60H10; 92D25

1. Introduction

Predator–prey interaction is an important relationship in population dynamics, which is a hot research topic in ecology and biomathematics. The study of predator–prey models can help us prevent and control large-scale species changes, protect species diversity, and optimize ecosystems.
An important component of the predator–prey relationship is the functional response of the predator. The functional response refers to the effect of the amount of predation by a predator in a certain period of time on the change of prey density. Numerous experiments [1,2] have shown that the functional response depends not only on the density of prey but also on the density of predators. Pal et al. [3] studied the predator–prey model using the Beddington–DeAngelis functional response, and they proved the asymptotic dynamic properties and explored the influence of fear on the stability of the system. Huang et al. [4] considered a delayed virus dynamics model with the Beddington–DeAngelis functional response and studied its stability analysis.
For a long time, predator–prey systems were studied only in terms of the effect of direct killing by predators [5,6]. However, experiments with song sparrows [7] and elk populations [8,9,10,11] have shown that fear of predators also affects prey populations, and sometimes the influence was even more serious than direct predation [12,13]. In addition, most studies have considered that fear reduces the reproductive rate of prey, whereas intraspecies competition is not affected [14,15,16]. In real life, fear may also have an effect on intraspecific competition for prey. Therefore, it is natural to include the fear of predator in our research model.
In fact, there exists time delay [17] in all kinds of biological processes, such as digestion of food, conversion of energy, gestation, maturation, inducible defense of prey groups, and so on. When prey senses a chemical signal or a sound signal, they need a certain amount of time to react, known as the fear delay. There is also a time lag between eating the prey and producing offspring, known as delayed gestation. However, few studies have investigated both fear delay and pregnancy delay [15,18,19,20]. Moreover, time delay has an important effect on the stability of the system [21,22,23]. Panday et al. [24] mainly studied the local stability of the system and the direction and stability of Hopf bifurcation. Furthermore, they found that when the time delay exceeded the critical value, the system developed Hopf bifurcation and became unstable. Kumar and Dubey [25] investigated a prey–predator system with prey refuge and gestation delay, and they proved the global asymptotic stability of the model and investigated the Hopf bifurcation behavior induced by the fear effect of prey.
Based on the above analysis, we propose a predator–prey model with double delays with fear and Beddington–DeAngelis functional response:
d x d t = k x 1 + f y ( t τ 1 ) α x 2 p x y a x + b y + c , d y d t = μ p x ( t τ 2 ) y ( t τ 2 ) a x ( t τ 2 ) + b y ( t τ 2 ) + c d y h y 2 ,
with the initial conditions
x ( θ ) = φ 1 ( θ ) , y ( θ ) = φ 2 ( θ ) , φ ( θ ) = ( φ 1 ( θ ) , φ 2 ( θ ) ) C ( [ τ , 0 ] , R + 2 ) ,
where τ 1 is the fear delay and τ 2 is the gestation delay, τ = max { τ 1 , τ 2 } , R + 2 = { ( x , y ) : x > 0 , y > 0 } , k denotes the intrinsic growth rate of prey species, f is the level of fear caused by predators, we denotes p by the per capita predator consumption rate, α is the decay rate of prey due to intraspecies competition, a is the time of the predator for each prey that is consumed, b measures the mutual interference between predators, c is the half-saturation constant of the predator population, μ is the conversion rate of prey biomass to predator, d is the natural mortality rate of predators, and h represented the death rate of predators from intraspecific competition. We define that | | φ | | = max { | φ ( θ ) | : θ [ τ , 0 ] } .
Because the parameters of model (1) are fixed and the influence of environmental noise is not considered, we incorporated random perturbation in the model (1) to investigate the impact of white noise on population dynamics. We establish the stochastic differential equation by perturbing the birth rate of prey population and natural death rate of predator population [26]. Then, we acquire the following stochastic model:
d x = k x 1 + f y ( t τ 1 ) α x 2 p x y a x + b y + c d t + σ 1 x 1 + f y ( t τ 1 ) d B 1 ( t ) , d y = μ p x ( t τ 2 ) y ( t τ 2 ) a x ( t τ 2 ) + b y ( t τ 2 ) + c d y h y 2 d t + σ 2 y d B 2 ( t ) ,
where B 1 ( t ) and B 2 ( t ) are standard and mutually independent Brownian motions defined on a complete probability space ( Ω , F , P ) with a filtration { F t } t 0 satisfying the usual conditions and σ i 2 ( i = 1 , 2 ) represents the intensity of the white noise.
The rest of this paper is organized as follows. In Section 2, we establish the positivity and boundedness of the solution of the delay model (1). In Section 3, we discuss the existence of equilibrium points, study the local stability of each equilibrium point and the Hopf bifurcation with different time delays. We consider the stochastic delay model (3) and prove the globally unique existence and stochastic ultimate boundedness of positive solutions in Section 4. We will perform some numerical simulations in Section 5 to verify our theoretical results. Finally, we draw conclusions based on our research findings and put forward some suggestions for future work in Section 6.

2. Positivity and Boundedness

To ensure that the model has a biological background, we study the positivity and boundedness of the delay Equation (1), that is, the solution is positive and invariant in the first quadrant and does not go beyond the given interval.

2.1. Positivity

Theorem 1. 
Let ( x ( t ) , y ( t ) ) be the solution of system (1) with the initial conditions (2). Then, ( x ( t ) , y ( t ) ) remains positive for any time t 0 .
Proof. 
System (1) is rewritten as a matrix
H ˙ = G ( H ) ,
where H = x y , G ( H ) = G 1 ( H ) G 2 ( H ) = k x 1 + f y ( t τ 1 ) α x 2 p x y a x + b y + c μ p x ( t τ 2 ) y ( t τ 2 ) a x ( t τ 2 ) + b y ( t τ 2 ) + c d y h y 2 with initial condition
H ( θ ) = ( φ 1 ( θ ) , φ 2 ( θ ) ) C ( [ τ , 0 ] , R + 2 ) , φ i ( 0 ) > 0 , i = 1 , 2 .
We can easily check the system (4) whenever taking H ( θ ) R + 2 such that x = y = 0 . Then, we obtain
G i ( H ) | h i = 0 , H R + 2 0
with h 1 ( t ) = x ( t ) , h 2 ( t ) = y ( t ) . According to Lemma 4 in [27], we derive that every solution of (4) with the initial condition (5) is positive, that is to say, any solution of the system (1) belongs to the region R + 2 and remains positive for any t 0 . □

2.2. Boundedness

Theorem 2. 
All solutions of system (1) that start in R + 2 are bounded.
Proof. 
Let W ( t ) = x ( t ) + y ( t + τ 2 ) μ . The time derivative of W ( t ) along the solution of (1) is
d W ( t ) d t = d x ( t ) d t + 1 μ d y ( t + τ 2 ) d t = k x 1 + f y ( t τ 1 ) α x 2 d μ y ( t + τ 2 ) h μ ( y ( t + τ 2 ) ) 2 k x α x 2 d μ y ( t + τ 2 ) .
Choose a constant σ such that σ < d . Then,
d W ( t ) d t + σ W ( t ) x ( k + σ α x ) d σ μ y ( t + τ 2 ) ( k + σ ) 2 4 α : = M .
Applying differential inequality theory results, we obtain
0 W ( t ) M σ ( 1 exp ( σ t ) ) + W ( x ( 0 ) , y ( 0 ) ) exp ( σ t ) .
As t , we have 0 < W ( t ) M σ . Therefore, all the solutions of system (1) are confined in the region
Z = { ( x , y ) R + 2 : 0 W ( t ) M σ } .

3. Stability Analysis

In this section, we mainly study the stability of equilibrium points and Hopf bifurcation around E * ( x * , y * ) of (1).

3.1. Equilibrium Points and Existence Criterion

The model (1) has three positive equilibrium points:
(1)
the trivial equilibrium point E 0 ( 0 , 0 ) ;
(2)
the free equilibrium point E 1 ( k α , 0 ) ;
(3)
the coexistence equilibrium point E * ( x * , y * ) satisfying the following equation:
k 1 + f y α x p y a x + b y + c = 0 , μ p x a x + b y + c d h y = 0 ,
By solving the above equation, we obtain x * = ( b y * + c ) ( d + h y * ) μ p a ( d + h y * ) . The interior equilibrium point exists only if μ p a ( d + h y * ) > 0 and y * is the equation A y 4 + B y 3 + C y 2 + D y + E = 0 , where
A = a 2 h 2 f p + b 2 α f h μ p , B = μ p b α ( 2 c f h + b d f + b h ) + a 2 h 2 p + 2 a h f p ( a d μ p ) , C = μ p b α ( b d + 2 c d f + 2 c h ) + h μ p ( a b k + c 2 α f ) + f p ( a d μ p ) 2 + 2 a h p ( a d μ p ) , D = c μ p ( 2 b d α + a h k + c α h + c d α f ) + μ p b k ( a d μ p ) + p ( a d μ p ) 2 , E = μ p c 2 d α + μ p c k ( a d μ p ) .
By the Descartes rule of sign and A > 0 , we find that the existence of at least one positive root when E < 0 is guaranteed, that is, k ( a d μ p ) + c d α < 0 . For more details, we refer the reader to [26].

3.2. Local Stability Analysis and Hopf Bifurcation of Equilibria

In this section, we perform the local stability on the dynamics of the model system (1) around the prescribed equilibrium points.
Using the transformations X = x x * , Y = y y * , the linearized form of system (1) is
d X ( t ) d t = J 0 X ( t ) + J 1 X ( t τ 1 ) + J 2 X ( t τ 2 ) ,
where X ( t ) = [ x ( t ) , y ( t ) ] T ,
J 0 = k 1 + f y * 2 α x * p y * ( b y * + c ) ( a x * + b y * + c ) 2 p x * ( a x * + c ) ( a x * + b y * + c ) 2 0 ( d + 2 h y * ) , J 1 = 0 k f x * ( 1 + f y * ) 2 0 0 , J 2 = 0 0 μ p y * ( b y * + c ) ( a x * + b y * + c ) 2 μ p x * ( a x * + c ) ( a x * + b y * + c ) 2 .
(1)
At E 0 :
J E 0 = k 0 0 d .
The Jacobian matrix has the eigenvalues k and d , which shows that E 0 is always unstable.
Therefore, we draw the following conclusion:
Theorem 3. 
The trivial equilibrium point E 0 is always unstable.
(2)
At E 1 :
J E 1 = k p k a k + c α f k 2 α e λ τ 1 0 d + μ p k a k + c α e λ τ 2 ,
and the characteristic equation is
( k λ ) ( d + μ p k a k + c α e λ τ 2 λ ) = 0 .
Hence, we can derive λ 1 = k , λ 2 = d + μ p k a k + c α e λ τ 2 .
Case I: τ 2 =0, λ 2 = d + μ p k a k + c α , E 1 is locally asymptotically stable if d > μ p k a k + c α .
Case II: τ 2 > 0 ,
d + μ p k a k + c α e λ τ 2 λ = 0 .
Assume that λ = i ξ ( ξ > 0 ) is a root of (8). Then, we have
d ( a k + c α ) + μ p k ( c o s ξ τ 2 i s i n ξ τ 2 ) i ξ ( a k + c α ) = 0 .
Then, we separate the real and imaginary parts as follows:
d ( a k + c α ) + μ p k c o s ξ τ 2 = 0 ,
μ p k s i n ξ τ 2 + ξ ( a k + c α ) = 0 .
By squaring and adding (9) and (10), we have
( a k + c α ) 2 ( d 2 + ξ 2 ) = ( μ p k ) 2 .
Therefore, we obtain
ξ 2 = ( μ p k a k + c α ) 2 d 2 .
If μ p k a k + c α d > 0 , the above Equation (8) has a positive root, which implies that E 1 is unstable;
If μ p k a k + c α d < 0 , Equation (8) contains one negative real root and imaginary roots with negative real parts, so E 1 is locally asymptotically stable for any τ 2 > 0 .
Therefore, we have the following conclusion:
Theorem 4. 
The free equilibrium point E 1 of the system (1) is locally asymptotically stable for any τ 2 0 if μ p k a k + c α d < 0 .
Remark 1. 
E 1 is locally asymptotically stable if a d + c d μ p > 0 , which contradicts the existence of E * . In other words, the instability of E 1 guarantees the existence of E * .
(3)
At E * :
J E * = a 11 b 12 e λ τ 1 + a 12 c 21 e λ τ 2 a 22 + c 22 e λ τ 2 ,
where a 11 = k 1 + f y * 2 α x * p y * ( b y * + c ) ( a x * + b y * + c ) 2 = α x * + a p x * y * ( a x * + b y * + c ) 2 , b 12 = f k x * ( 1 + f y * ) 2 , a 12 = p x * ( a x * + c ) ( a x * + b y * + c ) 2 , c 21 = μ p y * ( b y * + c ) ( a x * + b y * + c ) 2 , a 22 = ( d + 2 h y * ) , c 22 = μ p x * ( a x * + c ) ( a x * + b y * + c ) 2 .
Characteristic equation at E * is given by
λ 2 + A 1 λ + A 2 λ e λ τ 2 + A 3 e λ τ 2 + A 4 + A 5 e λ ( τ 1 + τ 2 ) = 0 ,
where A 1 = ( a 11 + a 22 ) , A 2 = c 22 , A 3 = a 11 c 22 a 12 c 21 , A 4 = a 11 a 22 , A 5 = b 12 c 21 .
Case I: τ 1 = 0 , τ 2 = 0 .
The characteristic Equation (11) becomes
λ 2 + ( A 1 + A 2 ) λ + A 3 + A 4 + A 5 = 0 .
According to the Routh–Hurwitz criteria, we can infer the J E * has a negative real part if A 1 + A 2 > 0 , A 3 + A 4 + A 5 > 0 , which implies a p y * α ( a x * + b y * + c ) 2 < 1 . Hence, E * is locally asymptotically stable if a p y * α ( a x * + b y * + c ) 2 < 1 .
Thus, we have the following theorem:
Theorem 5. 
In the absence of both τ 1 and τ 2 , the coexistence equilibrium point E * of the system (1) is locally asymptotically stable if a p y * α ( a x * + b y * + c ) 2 < 1 .
Case II: τ 1 > 0 , τ 2 = 0 . The characteristic Equation (11) becomes
λ 2 + ( A 1 + A 2 ) λ + A 3 + A 4 + A 5 e λ τ 1 = 0 .
Taking λ = i ξ * . Then, we obtain the real and imaginary parts, respectively, as follows:
( ξ * ) 2 + A 3 + A 4 + A 5 c o s ξ * τ 1 = 0 ,
( A 1 + A 2 ) ξ * A 5 s i n ξ * τ 1 = 0 .
Squaring (13) and (14) and adding, we obtain
( ξ * ) 4 + R 1 ( ξ * ) 2 + R 2 = 0 ,
where R 1 = ( A 1 + A 2 ) 2 2 ( A 3 + A 4 ) , R 2 = ( A 3 + A 4 ) 2 A 5 2 . By calculating, we have
( ξ * ) 2 = R 1 ± R 1 2 4 R 2 2 .
Let Z ( ξ * ) = ( ξ * ) 4 + R 1 ( ξ * ) 2 + R 2 , then Z ( 0 ) = R 2 < 0 , that is, ( A 3 + A 4 ) 2 < A 5 2 . Equation (16) has at least positive root ( ξ 0 * ) 2 . So, Equation (12) has a pair of imaginary roots ± i ξ 0 * . Substituting ( ξ 0 * ) 2 in (13) and (14), we obtain
τ 1 k = 1 ξ 0 * a r c t a n { ( A 1 + A 2 ) ξ 0 * ( ξ 0 * ) 2 + A 3 + A 4 } + 2 k π ξ 0 * , k = 0 , 1 , 2 ,
According to Butler’s lemma [28], E * remains locally asymptotically stable for 0 < τ 1 < τ 1 * ( = min k 0 τ 1 k ) and unstable if τ 1 > τ 1 * .
Now, we take the derivative of Equation (12) with respect to τ 1 as
d λ d τ 1 1 = 2 λ + A 1 + A 2 τ 1 A 5 e λ τ 1 λ A 5 e λ τ 1 = 2 λ + A 1 + A 2 λ A 5 e λ τ 1 τ 1 λ .
Furthermore, because e λ τ 1 = λ 2 + ( A 1 + A 2 ) λ + A 3 + A 4 A 5 , we can derive
d λ d τ 1 1 = 2 λ + A 1 + A 2 λ ( λ 2 + ( A 1 + A 2 ) λ + A 3 + A 4 ) τ 1 λ .
At τ 1 = τ 1 * , ξ * = ξ 0 * ,
d d τ 1 R e λ ( τ 1 ) 1 = ( ( A 1 + A 2 ) 2 2 ( A 3 + A 4 ) ) ξ 0 * + 2 ( ξ 0 * ) 4 ( A 1 + A 2 ) 2 ( ξ 0 * ) 4 + ( ξ 0 * ) 2 ( A 3 + A 4 ) 2 ( ξ 0 * ) 2 = A 5 2 ( A 3 + A 4 ) 2 + ( ξ 0 * ) 4 ( A 1 + A 2 ) 2 ( ξ 0 * ) 4 + ( ( ξ 0 * ) 2 ( A 3 + A 4 ) ) 2 ( ξ 0 * ) 2 > 0 .
Hence, the system (1) has Hopf bifurcation at τ 1 = τ 1 * .
Based on the above analysis, we have the following conclusion:
Theorem 6. 
In the absence of τ 2 and ( A 3 + A 4 ) 2 < A 5 2 , the coexistence equilibrium point E * is locally asymptotically stable for 0 < τ 1 < τ 1 * and unstable for τ 1 > τ 1 * . In addition, the system (1) will undergo a Hopf bifurcation at τ 1 = τ 1 * .
Case III: τ 1 = 0 , τ 2 > 0 . The characteristic Equation (11) becomes
λ 2 + A 1 λ + ( A 2 λ + A 3 + A 5 ) e λ τ 2 + A 4 = 0 .
Taking λ = i ξ ¯ * . Then, we obtain the real and imaginary parts, respectively, as follows:
( ξ ¯ * ) 2 + ( A 3 + A 5 ) c o s ξ ¯ * τ 2 + A 2 ξ ¯ * s i n ξ ¯ * τ 2 + A 4 = 0 ,
A 1 ξ ¯ * + A 2 ξ ¯ * c o s ξ ¯ * τ 2 ( A 3 + A 5 ) s i n ξ ¯ * τ 2 = 0 .
Squaring (19) and (20) and adding, we obtain
( ξ ¯ * ) 4 + R 3 ( ξ ¯ * ) 2 + R 4 = 0 .
Define R 3 = 2 A 4 + A 1 2 A 2 2 , R 4 = ( A 3 + A 5 ) 2 + A 4 2 . By calculating, we have
( ξ ¯ * ) 2 = R 3 ± R 3 2 4 R 4 2 .
Furthermore, because
R 3 = α x * + a p x * y * ( a x * + b y * + c ) 2 2 + ( h y * ) 2 + 2 ( d + h y * ) h y * + ( μ p x * ) 2 2 ( a x * + c ) b y * + ( b y * ) 2 ( a x * + b y * + c ) 4 > 0 .
When R 4 > 0 , Equation (21) has no positive roots and no real ξ ¯ * exists. So, E * is locally asymptotically stable for any τ 2 > 0 ;
When R 4 < 0 , Equation (22) has unique positive root ( ξ ¯ 0 * ) 2 . So, Equation (18) exists a pair of imaginary roots ± i ξ ¯ 0 * . Putting ( ξ ¯ 0 * ) 2 in (19) and (20), we obtain
τ 2 k = 1 ξ ¯ 0 * c o s 1 { ( ( ξ ¯ 0 * ) 2 A 4 ) A 2 ξ ¯ 0 * + A 1 ξ ¯ 0 * ( A 3 + A 5 ) A 2 ( ξ ¯ 0 * ) 2 + ( A 3 + A 5 ) 2 } + 2 k π ξ ¯ 0 * , k = 0 , 1 , 2 ,
According to Butler’s lemma, E * is locally asymptotically stable for 0 < τ 2 < τ 2 * ( = min k 0 τ 2 k ) and unstable if τ 2 > τ 2 * .
Now, take the derivative of Equation (18) with respect to τ 2
d λ d τ 2 1 = 2 λ + A 1 λ ( A 2 λ + A 3 + A 5 ) e λ τ 2 + A 2 λ ( A 2 λ + A 3 + A 5 ) τ 2 λ .
Furthermore, because e λ τ 2 = λ 2 + A 1 λ + A 4 A 2 λ + A 3 + A 5 , we can derive
d λ d τ 2 1 = 2 λ + A 1 λ ( λ 2 + A 1 λ + A 4 ) + A 2 λ ( A 2 λ + A 3 + A 5 ) τ 2 λ .
At τ 2 = τ 2 * , ξ ¯ * = ξ ¯ 0 * ,
d d τ 2 R e λ ( τ 2 ) 1 = A 1 2 + 2 ( ( ξ ¯ 0 * ) 2 A 4 ) ( A 1 ξ ¯ 0 * ) 2 + ( ( ξ ¯ 0 * ) 2 A 4 ) 2 A 2 2 ( A 2 ξ ¯ 0 * ) 2 + ( A 3 + A 5 ) 2 = A 1 2 A 2 2 2 A 4 + 2 ( ξ ¯ 0 * ) 2 ( A 2 ξ ¯ 0 * ) 2 + ( A 3 + A 5 ) 2 = R 3 + 2 ( ξ ¯ 0 * ) 2 A 2 2 ( ξ ¯ 0 * ) 2 + ( A 3 + A 5 ) 2 > 0 .
Hence, the system (1) undergoes Hopf bifurcation at τ 2 = τ 2 * .
According to the above analysis, we have the following theorem:
Theorem 7. 
In the absence of τ 1 , for system (1),
(1) If R 4 > 0 , the coexistence equilibrium point E * is locally asymptotically stable for any τ 2 > 0 ;
(2) If R 4 < 0 , the coexistence equilibrium point E * is locally asymptotically stable for 0 < τ 2 < τ 2 * and unstable if τ 2 > τ 2 * . In addition, the system (1) will undergo a Hopf bifurcation at τ 2 = τ 2 * .
Case IV: τ 1 is fixed in ( 0 , τ 1 * ) , τ 2 > 0 . Assume λ = i ξ ˜ , which is put in (11) and separate the real and imaginary parts as follows:
ξ ˜ 2 + A 2 ξ ˜ s i n ξ ˜ τ 2 + A 3 c o s ξ ˜ τ 2 + A 4 + A 5 ( c o s ξ ˜ τ 1 c o s ξ ˜ τ 2 s i n ξ ˜ τ 1 s i n ξ ˜ τ 2 ) = 0 ,
A 1 ξ ˜ + A 2 ξ ˜ c o s ξ ˜ τ 2 A 3 s i n ξ ˜ τ 2 A 5 ( s i n ξ ˜ τ 1 c o s ξ ˜ τ 2 + s i n ξ ˜ τ 2 c o s ξ ˜ τ 1 ) = 0 .
The above formulas can be arranged as:
( A 2 ξ ˜ A 5 s i n ξ ˜ τ 1 ) s i n ξ ˜ τ 2 + ( A 3 + A 5 c o s ξ ˜ τ 1 ) c o s ξ ˜ τ 2 = ξ ˜ 2 A 4 ,
( A 3 + A 5 c o s ξ ˜ τ 1 ) s i n ξ ˜ τ 2 + ( A 2 ξ ˜ A 5 s i n ξ ˜ τ 1 ) c o s ξ ˜ τ 2 = A 1 ξ ˜ .
Squaring (24) and (25) and adding to eliminate τ 2 , we have
ξ ˜ 4 + R 3 ξ ˜ 2 + R 5 ξ ˜ + R 6 = 0 ,
where R 5 = 2 A 2 A 5 s i n ξ ˜ τ 1 , R 6 = A 5 2 A 3 2 + A 4 2 2 A 3 A 5 c o s ξ ˜ τ 1 = 0 .
Define Z ( ξ ˜ ) = ξ ˜ 4 + R 3 ξ ˜ 2 + R 5 ξ ˜ + R 6 = 0 . Then, Z ( 0 ) = A 5 2 A 3 2 + A 4 2 2 A 3 A 5 = ( A 3 + A 5 ) 2 + A 4 2 , and Z ( ) = . Suppose Z ( 0 ) < 0 , that is to say, A 4 2 < ( A 3 + A 5 ) 2 . Then, Equation (26) has only one positive root. Hence, there exist the roots ± i ξ ˜ * in the characteristic Equation (11). From (24) and (25), we obtain
τ ˜ 2 k = 1 ξ ˜ * a r c s i n F 1 F 3 + F 2 F 4 F 1 2 + F 2 2 + 2 k π ξ ˜ * , k = 0 , 1 , 2 ,
where F 1 = A 2 ξ ˜ A 5 s i n ξ ˜ τ 1 , F 2 = A 3 + A 5 c o s ξ ˜ τ 1 , F 3 = ξ ˜ 2 A 4 , F 4 = A 1 ξ ˜ .
According to Butler’s lemma, E * is locally asymptotically stable for 0 < τ 2 < τ ˜ 2 * ( = min k 0 τ ˜ 2 k ) and unstable if τ ˜ 2 * .
Now, take the derivative of Equation (11) with respect to τ 2
d λ d τ 2 1 = 2 λ + A 1 + A 2 e λ τ 2 A 5 τ 1 e λ ( τ 1 + τ 2 ) λ ( A 2 λ + A 3 + A 5 e λ τ 1 ) e λ τ 2 τ 2 λ .
Furthermore, because e λ τ 2 = λ 2 + A 1 λ + A 4 A 2 λ + A 3 + A 5 e λ τ 2 , we can derive
d λ d τ 2 1 = 2 λ + A 1 λ ( λ 2 + A 1 λ + A 4 ) + A 2 A 5 τ 1 e λ τ 1 λ ( A 2 λ + A 3 + A 5 e λ τ 1 ) τ 2 λ .
At τ 2 = τ ˜ 2 * , ξ ˜ = ξ ˜ * ,
d d τ 2 R e λ ( τ 2 ) 1 = A 1 2 + 2 F 3 F 4 2 + F 3 2 + F 1 F 5 + F 6 F 2 ξ ˜ * ( F 1 2 + F 2 2 ) 2 ,
where F 5 = A 2 A 5 τ 1 c o s ξ ˜ * τ 1 , F 6 = A 5 τ 1 s i n ξ ˜ * τ 1 .
So, d d τ 2 R e λ ( τ 2 ) 1 > 0 leads to F 6 F 2 F 1 F 5 > 0 . Hence, there exists Hopf bifurcation at τ 2 = τ ˜ 2 * in the system (1).
Based on the above analysis, we have the following theorem:
Theorem 8. 
Assume that τ 1 is fixed in ( 0 , τ 1 * ) and A 4 2 < ( A 3 + A 5 ) 2 . Then, the coexistence equilibrium point E * is locally asymptotically stable for 0 < τ 2 < τ ˜ 2 * and unstable for τ 2 > τ ˜ 2 * . In addition, the system (1) will undergo a Hopf bifurcation at τ 2 = τ ˜ 2 * provided F 6 F 2 F 1 F 5 > 0 , where F i are all defined in the proof.

4. Stochastic Delay Model Analysis

4.1. Existence and Uniqueness of Positive Solution

In this subsection, we prove the unique existence of a global positive solution by means of a random comparison theorem.
Theorem 9. 
For any given initial value ( x ( 0 ) , y ( 0 ) ) R + 2 , there is a unique solution ( x ( t ) , y ( t ) ) to system (3) on t 0 . Furthermore, the solution will remain in R + 2 with probability 1, that is to say, ( x ( t ) , y ( t ) ) R + 2 almost surely.
Proof. 
Taking m = ln x , n = ln y , applying Itô’s formula we have
d m = k 1 + f e n ( t τ 1 ) α e m p e n a e m + b e n + c σ 1 2 2 ( 1 + f e n ( t τ 1 ) ) 2 d t + σ 1 1 + f e n ( t τ 1 ) d B 1 ( t ) , d n = μ p e m ( t τ 2 ) e n ( t τ 2 ) e n ( a e m ( t τ 2 ) + b e n ( t τ 2 ) + c ) d h e n σ 2 2 2 + σ 2 d B 2 ( t ) ,
where m ( 0 ) = ln x ( 0 ) and n ( 0 ) = ln y ( 0 ) . We notice that the coefficients of the above equation satisfy the local Lipschitz condition, so it possesses a unique local solution ( m ( t ) , n ( t ) ) on t [ 0 , τ e ) , where τ e is the explosion time. Hence, x = e m ( t ) , y = e n ( t ) is the unique positive local solution to (3) with initial value x 0 > 0 , y 0 > 0 . To show that the solution is global, we only need to prove τ e = .
According to the first equation of (3), we have
d x x ( t ) ( k α x ( t ) ) d t + σ 1 x ( t ) d B 1 ( t ) .
Let Υ ( t ) be the unique solution of the equation
d Υ ( t ) = Υ ( t ) ( k α Υ ( t ) ) d t + σ 1 Υ ( t ) d B 1 ( t ) , Υ ( 0 ) = x ( 0 ) ,
then,
Υ ( t ) = e ( k σ 1 2 2 ) t + σ 1 B 1 ( t ) 1 x ( 0 ) + α 0 t e ( k σ 1 2 2 ) s + σ 1 B 1 ( s ) d s .
Hence, by the comparison theorem [29] for stochastic equations, we obtain
x ( t ) Υ ( t ) .
According to the second equation of (3), we have
d y ( t ) y ( t ) ( d h y ( t ) ) d t + σ 2 y ( t ) d B 2 ( t ) .
Let Ψ ( t ) is the unique solution of the equation
d Ψ ( t ) = Ψ ( t ) ( d h Ψ ( t ) ) d t + σ 2 Ψ ( t ) d B 2 ( t ) , Ψ ( 0 ) = y ( 0 ) ,
then,
Ψ ( t ) = e ( d σ 2 2 2 ) t + σ 2 B 2 ( t ) 1 y ( 0 ) + h 0 t e ( d σ 2 2 2 ) s + σ 2 B 2 ( s ) d s .
Therefore,
y ( t ) Ψ ( t ) .
On the other hand,
d y ( t ) μ p Υ ( t τ 2 ) b d y ( t ) d t + σ 2 y ( t ) d B 2 ( t ) .
Let Φ ( t ) be the unique solution of the equation
d Φ ( t ) = μ p Φ ( t τ 2 ) b d y ( t ) d t + σ 2 y ( t ) d B 2 ( t ) , Φ ( 0 ) = y ( 0 ) .
Let N ( t ) = 1 Φ ( t ) with N ( 0 ) = 1 Φ ( 0 ) , by Itô’s formula
d N ( t ) = μ p Υ ( t τ 2 ) b N 2 ( t ) + d N ( t ) + σ 2 2 N ( t ) d t σ 2 N ( t ) d B 2 ( t ) = d + σ 2 2 μ p Υ ( t τ 2 ) b N ( t ) N ( t ) d t σ 2 N ( t ) d B 2 ( t ) ,
then,
N ( t ) = e ( d + σ 2 2 2 ) t σ 2 B 2 ( t ) Φ ( 0 ) + μ p b 0 t Υ ( s τ 2 ) e ( d + σ 2 2 2 ) s σ 2 B 2 ( s ) d s .
Hence, Φ ( t ) = 1 N ( t ) = Φ ( 0 ) e ( d + σ 2 2 2 ) t + σ 2 B 2 ( t ) + μ p b e ( d + σ 2 2 2 ) t + σ 2 B 2 ( t ) 0 t Υ ( s τ 2 ) e ( d + σ 2 2 2 ) s σ 2 B 2 ( s ) d s .
Similarly, we can derive
y ( t ) Φ ( t ) .
Again, from the first equation of (3), we have
d x ( t ) x ( t ) k 1 + f Φ ( t τ 1 ) α x p b d t + σ 1 x ( t ) 1 + f Φ ( t τ 1 ) d B 1 ( t ) .
We can safely infer
x ( t ) e p b t + 0 t k 1 + f Φ ( s τ 1 ) σ 1 2 2 ( 1 + f Φ ( s τ 1 ) ) 2 d s + 0 t σ 1 1 + f Φ ( s τ 1 ) d B 1 ( s ) 1 x ( 0 ) + α 0 t e p b s + 0 s k 1 + f Φ ( ρ τ 1 ) σ 1 2 2 ( 1 + f Φ ( ρ τ 1 ) ) 2 d ρ + 0 s σ 1 1 + f Φ ( ρ τ 1 ) d B 1 ( ρ ) d s : = Θ ( t ) ,
then, we have
x ( t ) Θ ( t ) .
Therefore, we have
Θ ( t ) x ( t ) Υ ( t ) , Ψ ( t ) y ( t ) Φ ( t ) .
Hence, τ e = , that is to say, the solution of (3) is global. This completes the proof. □

4.2. Stochastic Ultimate Boundedness

Definition 1. 
Consider a stochastic differential equation d X ( t ) = F ( X ( t ) ) d t + S ( X ( t ) ) d B ( t ) . Its solution is said to be stochastic ultimate bounded if for any ϵ ( 0 , 1 ) , there exists a positive Q = Q ( ϵ ) such that for any initial value ( x ( 0 ) , y ( 0 ) ) R + 2 , the solution satisfies
lim sup t P { | ( x ( t ) , y ( t ) ) | > Q } < ϵ .
Theorem 10. 
For any r ( 0 , 1 ) there exists a positive Q = Q ( r ) such that the solution ( x ( t ) , y ( t ) ) of the model (3) satisfies
lim sup t E | ( x ( t ) , y ( t ) ) | r Q ,
then we can derive the solution of the model (3) is stochastic and ultimately bounded.
Proof. 
Let V 1 ( x , y ) = x r + y r , from Itô’s formula we have
d V 1 ( x , y ) = L V 1 ( x , y ) d t + r x r σ 1 1 + f y ( t τ 1 ) d B 1 ( t ) + r σ 2 y r d B 2 ( t ) ,
where
L V 1 ( x , y ) = r x r k 1 + f y ( t τ 1 ) α x p y a x + b y + c + 1 2 r ( r 1 ) x r σ 1 2 ( 1 + f y ( t τ 1 ) ) 2 + r y r 1 μ p x ( t τ 2 ) y ( t τ 2 ) a x ( t τ 2 ) + b y ( t τ 2 ) + c d y h y 2 + 1 2 r ( r 1 ) σ 2 2 y r k r x r 1 + f y ( t τ 1 ) α r x r + 1 1 2 r ( 1 r ) x r σ 1 2 ( 1 + f y ( t τ 1 ) ) 2 + x r + r μ p x ( t τ 2 ) y ( t τ 2 ) a x ( t τ 2 ) + b y ( t τ 2 ) + c y r 1 h r y r + 1 d r y r 1 2 r ( 1 r ) σ 2 2 y r + y r V 1 ( x , y ) H 1 V 1 ( x , y ) ,
where H 1 is a positive constant. Hence, we can obtain
d V 1 ( x , y ) ( H 1 V 1 ( x , y ) ) d t + r x r σ 1 1 + f y ( t τ 1 ) d B 1 ( t ) + r σ 2 y r d B 2 ( t ) .
Define V 2 ( x , y ) = e t V 1 ( x , y ) , using Itô’s formula, we obtain
d V 2 ( x , y ) e t H 1 d t + r x r σ 1 1 + f y ( t τ 1 ) d B 1 ( t ) + r σ 2 y r d B 2 ( t ) .
Integrating both sides of the above inequality and taking the expectation, then
e t E V 1 ( x , y ) ( e t 1 ) H 1 + V 1 ( x ( 0 ) , y ( 0 ) ) .
So,
E V 1 ( x , y ) ( 1 e t ) H 1 + e t V 1 ( x ( 0 ) , y ( 0 ) ) .
Hence, we can obtain
lim sup t E V 1 ( x , y ) H 1 .
Because
| ( x ( t ) , y ( t ) ) | r = ( x 2 ( t ) + y 2 ( t ) ) r 2 2 r 2 max { x r , y r } 2 r 2 V 1 ( x , y ) ,
we can derive
lim sup t E | ( x ( t ) , y ( t ) ) | r 2 r 2 lim sup t E V 1 ( x , y ) 2 r 2 H 1 : = Q .
For any ϵ > 0 , let H 1 = ( H ϵ ) 1 r . According to the Markov inequality, we have
P { | ( x ( t ) , y ( t ) ) | > H 1 } < E | ( x ( t ) , y ( t ) ) | r H 1 r < ϵ .
The proof is complete. □

5. Numerical Simulations

In this section, as a general method in stochastic differential equations, we validate our mathematical findings by performing numerical simulations using Milstein’s high-order method [30] and MATLAB 2019a, that is, as proven in the previous sections. We select the following parameters:
f = 1.5 , p = 4 , a = 4.8 , b = 5 , c = 2.1 , μ = 1.2 , d = 0.1 , h = 0.01 , k = 0.9 , α = 0.6 ,
and take the initial data x ( 0 ) = 0.5 , y ( 0 ) = 0.5 . For the parameters (32), the system (1) has three equilibrium points, which are the trivial equilibrium point E 0 ( 0 , 0 ) , the free equilibrium point E 1 ( 1.5 , 0 ) , and the coexistence equilibrium point E * ( x * , y * ) = ( 0.1219 , 0.5702 ) . By Theorem (3), E 0 is always unstable. Furthermore, the existence of E * confirms that E 2 is unstable for any τ 1 0 and τ 2 0 .
Case I: Assume τ 1 = 0 and τ 2 = 0 and keep other parameters the same as (32). By simple calculation, a p y * α ( a x * + b y * + c ) 2 1 = 0.4047 < 0 . Then, the condition of Theorem 5 is satisfied, which implies that E * is locally asymptotically stable, as shown in Figure 1.
Case II: Assume τ 2 = 0 and slowly increase the value of τ 1 . Leave the other parameters remain (32) unchanged. By virtue of Theorem 6, we can obtain the critical value of τ 1 = τ 1 * = 4.5923 . When τ 1 = 3.5 < τ 1 * , the system (1) is locally asymptotically stable (as shown in Figure 2) and unstable for τ 1 = 5 > τ 1 * . Furthermore, the occurrence of oscillation behavior and limit cycle is illustrated by Figure 3. Furthermore, the system (1) undergoes a Hopf bifurcation at τ 1 = τ 1 * as shown in Figure 4.
Case III: Assume τ 1 = 0 and gradually increase the value of τ 2 , keeping other parameters remaining in (32) unchanged. After numerical calculation R 4 = 0.0015 < 0 , according to Theorem 7, we can obtain the critical value of τ 2 = τ 2 * = 2.3941 . When τ 2 = 1.5 < τ 2 * , the system (1) is locally asymptotically stable as shown in Figure 5 and unstable for τ 2 = 2.5 > τ 2 * , which is illustrated by Figure 6. Furthermore, the system (1) experiences a Hopf bifurcation at τ 2 = τ 2 * (as shown in Figure 7).
Case IV: We take a fixed τ 1 = 3.5 ( 0 , τ 1 * ) = ( 0 , 4.5923 ) and vary the parameter τ 2 while keeping other parameters the same as in (32). From Theorem 8, we can see that the critical value of τ 2 = τ ˜ 2 * = 0.4830 . When τ 2 = 0.2 < τ ˜ 2 * , the system (1) is locally asymptotically stable as shown in Figure 8 and unstable for τ 2 = 0.5 > τ ˜ 2 * , which is illustrated by Figure 9. Figure 10 reveals that a Hopf bifurcation occurs at τ 2 = τ ˜ 2 * .
For stochastic delay systems (3), according to Theorem 10, it is random and ultimately bounded. Figure 11 confirms our results.

6. Conclusions

In this paper, we mainly investigate a predator–prey model with Beddington–DeAngelis functional response and fear by predator. At the same time, because predator–prey interactions do not occur immediately, we introduced the fear delay and the pregnancy delay to make the model more natural.
First, we prove the positivity and boundedness of system (1). Then, we discuss the existence criterion and local stability of equilibrium point, and find that the existence of E * ensures that E 1 is unstable. Finally, the existence of a Hopf bifurcation with the fear delay τ 1 and the pregnancy delay τ 2 as the bifurcation parameter is studied, and the critical values of the bifurcation parameters are derived in several possible cases. We find the relationship between local stability and critical bifurcation value of the system (1). When the delay is less than the critical bifurcation value, both prey and predator oscillate periodically for a finite time and then reach equilibrium. When the delay exceeds the critical bifurcation value, Hopf bifurcation occurs in system (1), and periodic oscillation and limit cycle are generated. At this point, the system (1) switches from a stable state to an unstable state. Numerical simulations confirm our theoretical findings. Furthermore, for stochastic delay system (3), we study the unique existence of the global positive solution and explore the stochastic ultimate boundedness.
In future studies, we can generalize the model (1) to a multi-population model and introduce prey refuge to investigate its impact on the stability and persistence of population dynamics. In the meantime, it may be possible to study prey refuge as a Hopf bifurcation parameter. We leave these interesting questions for further study.

Author Contributions

Writing—original draft preparation, M.C.; writing—review and editing, Y.S., R.X. and J.Z.; visualization, M.C.; supervision, Y.S. and R.X.; funding acquisition: Y.S. and J.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work is supported by Science and Technology Program of Inner Mongolia Autonomous Region (2022YFHH0063, 2022YFHH0017) and National Natural Science Foundation of China (11861027, 12161062).

Data Availability Statement

Not applicable.

Acknowledgments

The authors are very grateful to the anonymous referees and the editors for their careful reading and valuable comments, which have helped improve the presentation of this work significantly.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. 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]
  2. Mondal, B.; Sarkar, S.; Ghosh, U. Complex dynamics of a generalist predator–prey model with hunting cooperation in predator. Eur. Phys. J. Plus 2022, 137, 43. [Google Scholar] [CrossRef]
  3. Pal, S.; Majhi, S.; Mandal, S.; Pal, N. Role of fear in a predator–prey model with Beddington–DeAngelis functional response. Z. FüR Nat. A 2019, 74, 581–595. [Google Scholar] [CrossRef]
  4. Huang, G.; Ma, W.; Takeuchi, Y. Global analysis for delay virus dynamics model with Beddington–DeAngelis functional response. Appl. Math. Lett. 2011, 24, 1199–1203. [Google Scholar] [CrossRef] [Green Version]
  5. Alam, S. Risk of disease-selective predation in an infected prey-predator system. J. Biol. Syst. 2009, 17, 111–124. [Google Scholar] [CrossRef]
  6. Roy, J.; Alam, S. Dynamics of an autonomous food chain model and existence of global attractor of the associated non-autonomous system. Int. J. Biomath. 2019, 12, 1950082. [Google Scholar] [CrossRef]
  7. 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]
  8. 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] [PubMed] [Green Version]
  9. Eggers, S.; Griesser, M.; Nystr, M.; Ekman, J. Predation risk induces changes in nest-site selection and clutch size in the Siberian jay. PRoceedings R. Soc. B Biol. Sci. 2006, 273, 701–706. [Google Scholar] [CrossRef] [Green Version]
  10. Orrock, J.L.; Fletcher, R.J., Jr. An island-wide predator manipulation reveals immediate and long-lasting matching of risk by prey. Proc. R. Soc. B Biol. Sci. 2014, 281, 20140391. [Google Scholar] [CrossRef]
  11. Sheriff, M.J.; Krebs, C.J.; Boonstra, R. The sensitive hare: Sublethal effects of predator stress on reproduction in snowshoe hares. J. Anim. Ecol. 2009, 78, 1249–1258. [Google Scholar] [CrossRef] [PubMed]
  12. Creel, S.; Christianson, D. Relationships between direct predation and risk effects. Trends Ecol. Evol. 2008, 23, 194–201. [Google Scholar] [CrossRef] [PubMed]
  13. 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] [PubMed]
  14. 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]
  15. Shao, Y. Global stability of a delayed predator–prey system with fear and Holling-type II functional response in deterministic and stochastic environments. Math. Comput. Simul. 2022, 200, 65–77. [Google Scholar] [CrossRef]
  16. Xia, Y.; Yuan, S. Survival analysis of a stochastic predator–prey model with prey refuge and fear effect. J. Biol. Dyn. 2020, 14, 871–892. [Google Scholar] [CrossRef]
  17. Liu, S.; Beretta, E.; Breda, D. Predator–prey model of Beddington–DeAngelis type with maturation and gestation delays. Nonlinear Anal. Real World Appl. 2010, 11, 4072–4091. [Google Scholar] [CrossRef]
  18. Bhattacharyyai, A.; Pal, A.K. Complex dynamics of delay induced two-prey one-predator model with Beddington-DeAngelis response function. Nonlinear Stud. 2021, 28. [Google Scholar]
  19. Li, J.; Liu, X.; Wei, C. The impact of fear factor and self-defence on the dynamics of predator-prey model with digestion delay. Math. Biosci. Eng. 2021, 18, 5478–5504. [Google Scholar] [CrossRef]
  20. Djilali, S.; Cattani, C.; Guin, L.N. Delayed predator–prey model with prey social behavior. Eur. Phys. J. Plus 2021, 136, 940. [Google Scholar] [CrossRef]
  21. Shu, H.; Hu, X.; Wang, L.; Watmough, J. Delay induced stability switch, multitype bistability and chaos in an intraguild predation model. J. Math. Biol. 2015, 71, 1269–1298. [Google Scholar] [CrossRef]
  22. Yuan, J.; Zhao, L.; Huang, C.; Xiao, M. Stability and bifurcation analysis of a fractional predator–prey model involving two nonidentical delays. Math. Comput. Simul. 2021, 181, 562–580. [Google Scholar] [CrossRef]
  23. Kong, W.; Shao, Y. The long time behavior of equilibrium status of a predator-prey system with delayed fear in deterministic and stochastic scenarios. J. Math. 2022, 2022, 3214358. [Google Scholar] [CrossRef]
  24. Pay, P.; Samanta, S.; Pal, N.; Chattopadhyay, J. Delay induced multiple stability switch and chaos in a predator–prey model with fear effect. Math. Comput. Simul. 2020, 172, 134–158. [Google Scholar]
  25. Kumar, A.; Dubey, B. Modeling the effect of fear in a prey–predator system with prey refuge and gestation delay. Int. J. Bifurc. Chaos 2019, 29, 1950195. [Google Scholar] [CrossRef]
  26. Roy, J.; Alam, S. Fear factor in a prey–predator system in deterministic and stochastic environment. Phys. A Stat. Mech. Its Appl. 2020, 541, 123359. [Google Scholar] [CrossRef]
  27. Yang, X.; Chen, L.; Chen, J. Permanence and positive periodic solution for the single-species nonautonomous delay diffusive models. Comput. Math. Appl. 1996, 32, 109–116. [Google Scholar] [CrossRef] [Green Version]
  28. Freedman, H.L.; Sree Hari Rao, V. The trade-off between mutual interference and time lags in predator-prey systems. Bull. Math. Biol. 1983, 45, 991–1004. [Google Scholar] [CrossRef]
  29. Peng, S.; Zhu, X. Necessary and sufficient condition for comparison theorem of 1-dimensional stochastic differential equations. Stoch. Process. Their Appl. 2006, 116, 370–380. [Google Scholar] [CrossRef] [Green Version]
  30. Higham, D.J. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev. 2001, 43, 525–546. [Google Scholar] [CrossRef]
Figure 1. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 0 and τ 2 = 0 . Other parameters are the same as in (32).
Figure 1. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 0 and τ 2 = 0 . Other parameters are the same as in (32).
Axioms 12 00073 g001
Figure 2. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0 . Other parameters are the same as in (32).
Figure 2. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0 . Other parameters are the same as in (32).
Axioms 12 00073 g002
Figure 3. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 5 and τ 2 = 0 . Other parameters are same as in (32).
Figure 3. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 5 and τ 2 = 0 . Other parameters are same as in (32).
Axioms 12 00073 g003
Figure 4. The Hopf bifurcation diagram of (1) for τ 1 and τ 2 = 0 .
Figure 4. The Hopf bifurcation diagram of (1) for τ 1 and τ 2 = 0 .
Axioms 12 00073 g004
Figure 5. (a) Time series and (b) phase portrait of E * for system (1) when τ 2 = 1.5 and τ 1 = 0 . Other parameters are the same as in (32).
Figure 5. (a) Time series and (b) phase portrait of E * for system (1) when τ 2 = 1.5 and τ 1 = 0 . Other parameters are the same as in (32).
Axioms 12 00073 g005
Figure 6. (a) Time series and (b) phase portrait of E * for system (1) when τ 2 = 2.5 and τ 1 = 0 . Other parameters are the same as in (32).
Figure 6. (a) Time series and (b) phase portrait of E * for system (1) when τ 2 = 2.5 and τ 1 = 0 . Other parameters are the same as in (32).
Axioms 12 00073 g006
Figure 7. The Hopf bifurcation diagram of (1) for τ 2 and τ 1 = 0 .
Figure 7. The Hopf bifurcation diagram of (1) for τ 2 and τ 1 = 0 .
Axioms 12 00073 g007
Figure 8. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0.2 . Other parameters are the same as in (32).
Figure 8. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0.2 . Other parameters are the same as in (32).
Axioms 12 00073 g008
Figure 9. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0.5 . Other parameters are the same as in (32).
Figure 9. (a) Time series and (b) phase portrait of E * for system (1) when τ 1 = 3.5 and τ 2 = 0.5 . Other parameters are the same as in (32).
Axioms 12 00073 g009
Figure 10. The Hopf bifurcation diagram of (1) for τ 2 and τ 1 = 3.5 .
Figure 10. The Hopf bifurcation diagram of (1) for τ 2 and τ 1 = 3.5 .
Axioms 12 00073 g010
Figure 11. Dynamical behaviors of (3). The parameters are taken as (32) and σ 1 = 0.03 , σ 2 = 0.03 .
Figure 11. Dynamical behaviors of (3). The parameters are taken as (32) and σ 1 = 0.03 , σ 2 = 0.03 .
Axioms 12 00073 g011
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Cui, M.; Shao, Y.; Xue, R.; Zhao, J. Dynamic Behavior of a Predator–Prey Model with Double Delays and Beddington–DeAngelis Functional Response. Axioms 2023, 12, 73. https://doi.org/10.3390/axioms12010073

AMA Style

Cui M, Shao Y, Xue R, Zhao J. Dynamic Behavior of a Predator–Prey Model with Double Delays and Beddington–DeAngelis Functional Response. Axioms. 2023; 12(1):73. https://doi.org/10.3390/axioms12010073

Chicago/Turabian Style

Cui, Minjuan, Yuanfu Shao, Renxiu Xue, and Jinxing Zhao. 2023. "Dynamic Behavior of a Predator–Prey Model with Double Delays and Beddington–DeAngelis Functional Response" Axioms 12, no. 1: 73. https://doi.org/10.3390/axioms12010073

APA Style

Cui, M., Shao, Y., Xue, R., & Zhao, J. (2023). Dynamic Behavior of a Predator–Prey Model with Double Delays and Beddington–DeAngelis Functional Response. Axioms, 12(1), 73. https://doi.org/10.3390/axioms12010073

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