Next Article in Journal
X-ray Self-Emission Imaging of Hydrodynamic Laser-Induced Astrophysical Phenomena
Next Article in Special Issue
Symmetries of Non-Linear ODEs: Lambda Extensions of the Ising Correlations
Previous Article in Journal
Be Careful What You Wish for: Cost Function Sensitivity in Predictive Simulations for Assistive Device Design
Previous Article in Special Issue
A Computational Scheme for the Numerical Results of Time-Fractional Degasperis–Procesi and Camassa–Holm Models
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Qualitative Analysis in a Beddington–DeAngelis Type Predator–Prey Model with Two Time Delays

1
School of Mathematical Sciences, Jiangsu University, Zhenjiang 212013, China
2
School of Medical Informatics and Engineering, Xuzhou Medical University, 209 Tongshan Road, Xuzhou 221004, China
3
Department of Basic Science, Obour High Institute for Engineering and Technology, Cairo 11828, Egypt
*
Author to whom correspondence should be addressed.
Symmetry 2022, 14(12), 2535; https://doi.org/10.3390/sym14122535
Submission received: 30 October 2022 / Revised: 21 November 2022 / Accepted: 21 November 2022 / Published: 30 November 2022

Abstract

:
In this paper, we investigate a delayed predator–prey model with a prey refuge where the predator population eats the prey according to the Beddington–DeAngelis type functional response. Firstly, we consider the existence of equilibrium points. By analyzing the corresponding characteristic equations, the local stability of the trivial equilibrium, the predator–extinction balance, and the coexistence equilibrium of the system are discussed, and the existence of Hopf bifurcations concerning both delays at the coexistence equilibrium are established. Then, in accordance with the standard form method and center manifold theorem, the explicit formulas which determine the direction of Hopf bifurcation and stability of bifurcating period solutions are derived. Finally, representative numerical simulations are performed to validate the theoretical analysis.

1. Introduction

The dynamic relationship between predator and prey has always existed in our ecosystem due to its ubiquity and significance. The analysis of this relationship has turned into one of the research hotspots in biomathematics. Since Lotka [1] and Volterra [2] proposed the most concise model for predator–prey interaction, the classical two-dimensional Lotka–Volterra system has then been described by
d x d t = a 0 x b 0 x y , d y d t = b 0 x y c 0 y ,
where x and y denote the densities of prey and predator population, and the parameters a 0 , b 0 , and c 0 are three positive parameters representing the interaction of two populations. In view of this model, many scholars studied the relevant dynamics by considering unique structures, such as time delay term, functional response function term, prey refuge, etc. [3,4,5,6,7,8,9,10].
Functional response is an essential factor in the predator–prey model, which reflects the reaction of each predator to prey density. Different applicable response functions can induce various dynamic behaviors, which can explain the complexity of the ecosystem. Based on Holling’s experiments [11], three diverse functional response functions were introduced concerning different species: Holling I, Holling II, and Holling III. These functional response functions are only correlated with the density of prey but not with the density of predator, which can be known as the prey-dependent functional response. In [12], Liu et al. studied a predator–prey model with Holling type II functional response α x 1 1 + w x 1 , where parameter α (units: 1/time) and w (units: 1/prey) are positive constants which represent the influence of capture rate and handling time. Under this circumstance, predators do not interfere with each other’s activities. In reality, there will be competition or cooperation among the predators. The ideal functional response function should be related not only to the density of prey but also to the density of the predator.
Beddington [13] and DeAngelis [14] established a functional response which is defined by Beddington–DeAngelis functional response, expressed as q x a + b x + c y . Through comparison, we find that it is akin to Holling type II functional response, which can be described as its generalization. Several scholars have studied the predator–prey model with the Beddington–DeAngelis functional response function and obtained the corresponding conclusions, such as the existence of limit cycles, bifurcation problems, and so on [15,16,17,18,19]. In [15], Hwang reported a predator–prey system with the Beddington–DeAngelis functional response, and the uniqueness of limit cycles was proved. Tripathi et al. [17] explored a Beddington–DeAngelis type one-predator two-prey competitive approach with help and presented relevant local stability, uniform persistence, and permanence.
Furthermore, considering another essential factor from an ecological point of view, if all prey can be captured, then the prey population may be on the verge of extinction. Due to the continuous evolution of organisms, researchers have observed that, to avoid predation better, prey have developed many defensive abilities, such as color protection, concealment, and so on, that is, the situation where the prey refuge is produced [4,20,21,22,23,24]. Tang et al. [20] investigated the global analysis of a Holling type II predator–prey model with a prey refuge. In [22], Peng et al. established a predator–prey system with Holling type III functional response by introducing a prey refuge and selective harvesting.
In addition, due to the lag phenomenon of population interaction, this lag exists in the stage of each population. Therefore, it is necessary to consider lag in studying population dynamics, which can more reasonably reflect the interaction between populations. Some scholars have gradually considered time delay in the population dynamics model, including feedback delay [25,26], gestation delay [27,28], etc. Moreover, the diverse types of delays have been presented [29,30,31], in which they discussed the local stability, Hopf bifurcation, and the bifurcation properties.
Based on the discussions above and inspired by the work [12], we first establish a predator–prey model with Beddington–DeAngelis functional response and incorporate a prey refuge. Then we consider two time delays: one is feedback delay about the prey population, and the other is gestation delay about the predator population, which is in the following form:
x ˙ = x ( t ) r ρ 1 x t τ 1 α 1 ( 1 m ) x ( t ) y ( t ) a 1 + b 1 ( 1 m ) x ( t ) + c 1 y ( t ) , y ˙ = α 2 ( 1 m ) x t τ 2 y t τ 2 a 1 + b 1 ( 1 m ) x t τ 2 + c 1 y t τ 2 d 1 y ( t ) ρ 2 y 2 ( t ) .
where x ( t ) and y ( t ) can be understood as the population densities of prey and predator at time t, respectively. All the parameters r, α 1 , α 2 , a 1 , b 1 , c 1 , ρ 1 , ρ 2 , d 1 , m are positive constants. We also regard the assumptions as below:
(i)
For prey: r denotes an intrinsic growth rate, ρ 1 means the intraspecific competition rate, m [ 0 , 1 ) denotes refuge parameter, and ( 1 m ) x ( t ) number of prey are available to predator, τ 1 denotes the feedback time delay of the prey.
(ii)
For predator: α 1 ( 1 m ) shows the capturing rate, α 2 α 1 denotes the conversion rate of the predator, d 1 means the death rate, ρ 2 denotes the intraspecific competition rate, the predator consumes the prey with the Beddington–DeAngelis functional response and prey refuge α 1 ( 1 m ) x ( t ) y ( t ) a 1 + b 1 ( 1 m ) x ( t ) + c 1 y ( t ) . τ 2 is the time delay due to the gestation of the predator.
The initial conditions for the system (1) is given by
x ( 0 ) > 0 , y ( 0 ) > 0 .
In the light of the fundamental theory of functional differential equations [32,33], the system (1) has a unique solution ( x ( t ) , y ( t ) ) contenting initial conditions (2). It can easily be shown that all solutions of system (1) with initial conditions (2) are defined on [ 0 , + ) and keep positive for whole τ 1 > 0 and τ 2 > 0 .
For simplicity, by the re-scaling t r t , x ρ 1 r , y y , β 1 = α 1 r , β 2 = α 2 ρ 1 , a = a 1 , b = b 1 r ρ 1 , c = c 1 , d = d 1 r , e = ρ 2 r , then system (1) in a non-dimensional form is:
x ˙ = x ( t ) 1 x t τ 1 β 1 ( 1 m ) x ( t ) y ( t ) a + b ( 1 m ) x ( t ) + c y ( t ) y ˙ = β 2 ( 1 m ) x t τ 2 y t τ 2 a + b ( 1 m ) x t τ 2 + c y t τ 2 d y ( t ) e y 2 ( t ) .
The framework of this paper is arranged below. In Section 2, we study the existence of each possible equilibrium point and their local stability analysis of equilibria and Hopf bifurcation. In Section 3, we obtain the nature of Hopf bifurcation by utilizing the standard form method and center manifold theorem. In Section 4, some numerical simulations are conducted to check the theoretical analysis. Finally, we give a concise summary of this paper.

2. Stability Analysis and Hopf Bifurcation

In this section, the existence of equilibria and their local stability, as well as the occurrence of Hopf bifurcation will be explored.
The system (3) has three non-negative equilibrium points.
(i)
It is not difficult to verify that the system invariably exists a trivial equilibrium E 0 ( 0 , 0 ) .
(ii)
The predator–extinction equilibrium the E 1 ( 1 , 0 ) biological explanation for this boundary equilibrium is that, without predators, the number of prey reaches its carrying capacity.
(iii)
The positive equilibrium is defined by E 2 ( x * , y * ) ; note that E 2 satisfies the following equations:
1 x * β 1 ( 1 m ) y * a + b ( 1 m ) x * + c y * = 0 , β 2 ( 1 m ) x * a + b ( 1 m ) x * + c y * d e y * = 0 .
With the help of [16], if and only if the condition:
(H1) ( 1 m ) β 2 b d > a d , holds, then system (3) possesses a unique positive equilibrium E 2 ( x * , y * ) .
In what follows, we analyse the local stability about the equilibria E 0 , E 1 , E 2 , respectively.
To begin with, we discuss the local stability about the trivial equilibrium E 0 .
Theorem 1. 
The equilibrium E 0 is unstable.
Proof. 
The characteristic equation of system (3) at E 0 is
( λ 1 ) ( λ + d ) = 0 .
Obviously, Equation (5) has a positive root λ 1 = 1 , indicating that E 0 is unstable equilibrium. This completes the proof. □
Then, we investigate the stability concerning the predator-extinction equilibrium E 1 .
Theorem 2. 
For system (3),
(i) 
If (H1) ( 1 m ) ( β 2 b d ) > a d holds, then the equilibrium E 1 is unstable for total τ 2 0 .
(ii) 
If (H2) ( 1 m ) ( β 2 b d ) < a d holds, then the equilibrium E 1 is locally asymptotically stable when τ 2 = 0 , and for whole τ 2 > 0 , the E 1 is unstable.
Proof. 
The characteristic equation of system (3) at E 1 turns into
λ + e λ τ 1 λ + d β 2 ( 1 m ) e λ τ 2 a + b ( 1 m ) = 0 .
Evidently, λ 1 = e λ τ 1 always keep negative, which indicates that the other root of Equation (6) is decided by
λ + d β 2 ( 1 m ) e λ τ 2 a + b ( 1 m ) = 0 .
Because the above equation exists time delay τ 2 , we provide the following two cases.
For τ 2 = 0 , we obtain λ = d + β 2 ( 1 m ) a + b ( 1 m ) = ( 1 m ) β 2 b d a d , if (H2) holds, according to the Routh–Hurwitz criterion, E 1 is locally asymptotically stable. Furthermore, if (H1) is satisfied, E 1 is unstable.
When τ 2 0 , we consider i ω 0 ω 0 > 0 be a root of Equation (7), and by a calculation, it can obtain that
β 2 ( 1 m ) a + b ( 1 m ) sin ω 0 τ 2 = ω 0 , β 2 ( 1 m ) a + b ( 1 m ) cos ω 0 τ 2 = d ,
which implies
ω 0 2 + d 2 β 2 ( 1 m ) a + b ( 1 m ) 2 = 0 .
In above equation, the coefficient of ω 0 is zero, then in the light of the Routh–Hurwitz criterion, when τ 2 > 0 , the equilibrium E 1 is always unstable. The proof is complete. □
Further, we discuss the stability of a unique positive equilibrium E 2 .
Consider x ^ ( t ) = x ( t ) x * , y ^ ( t ) = y ( t ) y * , and denote x ^ ( t ) , y ^ ( t ) by x ( t ) , y ( t ) , respectively. Taylor expansion is utilized to spread the system (3) at the equilibrium E 2 ; thus, we obtain
x ˙ ( t ) = a 11 x ( t ) + b 11 x t τ 1 + a 12 y ( t ) + i + j + k 2 f 1 ( i j k ) x i ( t ) y j ( t ) x k t τ 1 y ˙ ( t ) = a 21 x t τ 2 + a 22 y ( t ) + b 22 y t τ 2 + i + j + k 2 f 2 ( i j k ) x i t τ 2 y j ( t ) y k t τ 2
where
a 11 = 1 x * β 1 ( 1 m ) y * a + c y * a + b ( 1 m ) x * + c y * 2 , b 11 = x * , a 12 = β 1 ( 1 m ) x * a + b ( 1 m ) x * a + b ( 1 m ) x * + c y * 2 , a 21 = β 2 ( 1 m ) y * a + c y * a + b ( 1 m ) x * + c y * 2 , a 22 = d 2 e y * , b 22 = β 2 ( 1 m ) x * a + b ( 1 m ) x * a + b ( 1 m ) x * + c y * 2 , f 1 ( i j k ) = 1 i ! j ! k ! i + j + k f 1 x i ( t ) y j ( t ) x k t τ 1 x * , y * , f 2 ( i j k ) = 1 i ! j ! k ! i + j + k f 2 x i t τ 2 y j ( t ) y k t τ 2 x * , y * , f 1 = x ( t ) 1 x t τ 1 β 1 ( 1 m ) x ( t ) y ( t ) a + b ( 1 m ) x ( t ) + c y ( t ) , f 2 = β 2 ( 1 m ) x t τ 2 y t τ 2 a + b ( 1 m ) x t τ 2 + c y t τ 2 d y ( t ) e y 2 ( t ) .
Then, we linearize the system (9) as follows
x ˙ ( t ) = a 11 x ( t ) + b 11 x t τ 1 + a 12 y ( t ) , y ˙ ( t ) = a 21 x t τ 2 + a 22 y ( t ) + b 22 y t τ 2 .
The corresponding characteristic equation of system (10) becomes
λ 2 + A λ + B + ( C λ + D ) e λ τ 1 + ( E λ + F ) e λ τ 2 + G e λ τ 1 + τ 2 = 0 ,
where A = a 11 + a 22 , B = a 11 a 22 , C = b 11 , D = a 22 b 11 , E = b 22 , F = a 11 b 22 a 12 a 21 , G = b 11 b 22 .
To analyze the root distribution of the transcendental Equation (11) in the following, we give the result of Ruan and Wei [34].
Lemma 1. 
Regard to the transcendental equation
p λ , e λ τ 1 , , e λ τ m = λ n + p 1 ( 0 ) λ n 1 + + p n 1 ( 0 ) λ + p n ( 0 ) + p 1 ( 1 ) λ n 1 + + p n 1 ( 1 ) λ + p n ( 1 ) e λ τ 1 + + p 1 ( m ) λ n 1 + + p n 1 ( m ) λ + p n ( m ) e λ τ m = 0 ,
when ( τ 1 , τ 2 , τ 3 , , τ m ) change, the sum of orders of the zeros of p λ , e λ τ 1 , , e λ τ m in the open right half plane can vary, and only a zero emerges on or passes through the imaginary axis.
Due to the existence about τ 1 and τ 2 , we analyze six different situations as follows.
Situation 1: τ 1 = τ 2 = 0 . Equation (11) turns to
λ 2 + A 1 λ + B 1 = 0 ,
where A 1 = A + C + E , B 1 = B + D + F + G .
Then total roots of Equation (12) have negative real parts if and only if (H3) A 1 > 0 , B 1 > 0 .
Hence, the positive equilibrium E 2 is locally asymptotically stable when (H3) is satisfied.
Situation 2: τ 1 = 0 , τ 2 > 0 ; then Equation (11) reduces to
λ 2 + A 2 λ + B 2 + C 2 λ + D 2 e λ τ 2 = 0 ,
where A 2 = A + C , B 2 = B + D , C 2 = E , D 2 = F + G .
Let λ = i ω 2 ( ω 2 > 0 ) be a solution of Equation (13), which implies that
D 2 sin ω 2 τ 2 C 2 ω 2 cos ω 2 τ 2 = A 2 ω 2 , D 2 cos ω 2 τ 2 + C 2 ω 2 sin ω 2 τ 2 = ω 2 2 B 2 ,
which yields
ω 2 4 + p 21 ω 2 2 + p 22 = 0 ,
where p 21 = A 2 2 2 B 2 C 2 2 , p 22 = B 2 2 D 2 2 .
Making ω 2 2 = v 1 , then Equation (15) can be re-written by
v 1 2 + p 21 v 1 + p 22 = 0 ,
We can observe that if (H4) p 21 > 0 , p 22 > 0 is satisfied, then Equation (15) has no positive root, indicating that total roots of Equation (13) have negative real parts. Consequently, E 2 is locally asymptotically stable when τ 2 0 under the conditions (H3) and (H4).
If Conditions (H3) and (H5) p 21 > 0 , p 22 < 0 are satisfied, indicating that Equation (15) exists a positive root ω 2 2 , by calculation, we obtain
τ 2 ( n ) = 1 ω 2 arccos D 2 A 2 C 2 ω 2 2 B 2 D 2 D 2 2 + C 2 2 ω 2 2 + 2 n π ω 2 , n = 0 , 1 , 2 ,
On the basis of the Hopf bifurcation theorem [35,36,37], we require to prove the transversality condition.
Lemma 2. 
If Condition (H6) f ω 20 2 0 holds, then d ( Re λ ) d τ 2 λ = i ω 20 0 .
Proof. 
Differentiating both sides of Equation (13) with regard to τ 2 , we obtain
d λ d τ 2 = 2 λ + A 2 λ λ 2 + A 2 λ + B 2 + C 2 λ C 2 λ + D 2 τ 2 λ ,
which leads to
Re d λ d τ 2 λ = i ω 20 1 = Re 2 λ + A 2 λ λ 2 + A 2 λ + B 2 + Re C 2 λ C 2 λ + D 2 = 2 ω 20 2 + A 2 2 2 B 2 ω 20 2 B 2 2 + A 2 2 ω 20 2 + C 2 2 C 2 2 ω 20 2 + D 2 2 .
We obtain from Equation (14) that
ω 20 2 B 2 2 + A 2 2 ω 20 2 = C 2 2 ω 20 2 + D 2 2 .
Noting that d ( Re λ ) d τ 2 λ = i ω 20 and Re ( d λ d τ 2 ) 1 λ = i ω 20 have the same sign, then
sign d ( Re λ ) d τ 2 λ = i ω 20 = sign Re d λ d τ 2 1 λ = i ω 20 = 2 ω 20 2 + A 2 2 2 B 2 C 2 2 C 2 2 ω 20 2 + D 2 2 = f ω 20 2 C 2 2 ω 20 2 + D 2 2 .
Therefore, d ( Re λ ) d τ 2 λ = i ω 20 0 if (H6) f ω 20 2 0 holds. The proof of this lemma is now finished.
Through the above research, the corresponding results are concluded as below.
Theorem 3. 
For system (3), when τ 1 = 0 , suppose that (H1) and (H3) are satisfied,
(i) 
If (H4) holds, then the equilibrium E 2 ( x * , y * ) is asymptotically stable for total τ 2 0 .
(ii) 
If (H5) and (H6) hold, then the equilibrium is E 2 asymptotically stable for whole τ 2 [ 0 , τ 20 ) and unstable for [ τ 2 > τ 20 ] . Additionally, the system occurs a Hopf bifurcation at the positive equilibrium E 2 when τ 2 = τ 20 .
Situation 3: τ 1 > 0 , τ 2 = 0 .
Theorem 4. 
For system (3), the positive equilibrium is E 2 ( x * , y * ) asymptotically stable for total τ 2 [ 0 , τ 20 ) and unstable for τ 1 > τ 10 . Additionally, the system (3) goes through a Hopf bifurcation at the equilibrium E 2 when τ 1 = τ 10 , where τ 10 means the minimum critical value of time delay τ 1 when Hopf bifurcation occurs.
Proof. 
Because we adopt a same method which is applied in Situation 2, here we omit the process. □
Situation 4: τ 1 = τ 2 = τ , then Equation (11) can be written by
λ 2 + A 4 λ + B 4 + C 4 λ + D 4 e λ τ + E 4 e 2 λ τ = 0 ,
where A 4 = A , B 4 = B , C 4 = C + E , D 4 = D + F , E 4 = G .
Multiplying both sides of Equation (17) of e λ τ , we obtain
C 4 λ + D 4 + λ 2 + A 4 λ + B 4 e λ τ + E 4 e λ τ = 0 .
Let i ω ( ω > 0 ) be a root of Equation (18); thus, we have
A 41 sin ( ω τ ) + A 42 cos ( ω τ ) = A 45 , A 43 cos ( ω τ ) + A 44 sin ( ω τ ) = A 46 ,
where A 41 = ω 2 + B 4 E 4 , A 42 = A 4 ω , A 43 = ω 2 + B 4 + E 4 , A 44 = A 4 ω , A 45 = C 4 ω , A 46 = D 4 .
This results in
sin ( ω τ ) = B 41 ω 3 + B 42 ω ω 4 + B 45 ω 2 + B 46 , cos ( ω τ ) = B 43 ω 2 + B 44 ω 4 + B 45 ω 2 + B 46 ,
where B 41 = C 4 , B 42 = A 4 D 4 B 4 C 4 C 4 E 4 , B 43 = D 4 A 4 C 4 , B 44 = D 4 E 4 B 4 D 4 , B 45 = 2 B 4 + A 4 2 , B 46 = B 4 E 4 2 .
Then we obtain
ω 8 + p 41 ω 6 + p 42 ω 4 + p 43 ω 2 + p 44 = 0 ,
with
p 41 = B 41 2 + 2 B 45 , p 42 = B 45 2 B 43 2 + 2 B 46 2 B 41 B 42 p 43 = 2 B 45 B 46 B 42 2 2 B 43 B 44 , p 44 = B 46 2 B 44 2 .
Consider ω 2 = v 2 , then, from Equation (21), we obtain
v 2 4 + p 41 v 2 3 + p 42 v 2 2 + p 43 v 2 + p 44 = 0 .
In general, we suppose Equation (22) has four positive roots, which are described by v 1 , v 2 , v 3 , and v 4 . Then, Equation (21) also exists four positive roots ω k = v 2 k , k = 1 , 2 , 3 , 4 . The corresponding critical value of time delay τ 2 k ( j ) is expressed by
τ k ( n ) = 1 ω k arccos B 43 ω k 2 + B 44 ω k 4 + B 45 ω k 2 + B 46 + 2 n π ω k , k = 1 , 2 , 3 , 4 ; n = 0 , 1 , 2 , .
Let λ ( τ ) = ς ( τ ) + i ω ( τ ) be a root of Equation (18) comprising ς τ k j = 0 , ω τ k j = ω k , and define τ 0 = min k 1 , 2 , 3 , 4 τ k ( n ) , ω 0 = ω k .
Lemma 3. 
If Condition (H7) Q 41 Q 43 + Q 42 Q 44 0 holds, then d ( Re λ ) d τ λ = i ω 0 0 .
Proof. 
Differentiating both sides of Equation (18) about τ , and by simplification, we have
d λ d τ 1 = C 4 + 2 λ + A 4 e λ τ λ C 4 λ + D 4 + 2 λ E 4 e λ τ τ λ ,
It follows that
Re d λ d τ λ = i ω 0 1 = Re C 4 + 2 λ + A 4 e λ τ λ C 4 λ + D 4 + 2 λ E 4 e λ τ = Re Q 41 + Q 42 i Q 43 + Q 44 i = Q 41 Q 43 + Q 42 Q 44 Q 43 2 + Q 44 2 ,
where
Q 41 = C 4 + A 4 cos ω 0 τ 0 2 ω 0 sin ω 0 τ 0 , Q 42 = A 4 sin ω 0 τ 0 + 2 ω 0 cos ω 0 τ 0 , Q 43 = 2 C 4 E 4 ω 0 sin ω 0 τ 0 2 D 4 E 4 ω 0 2 cos ω 0 τ 0 , Q 44 = 2 C 4 E 4 ω 0 3 cos ω 0 τ 0 + 2 D 4 E 4 ω 0 sin ω 0 τ 0 .
Evidently, if (H7) Q 41 Q 43 + Q 42 Q 44 0 is satisfied, we can obtain d ( Re λ ) d τ λ = i ω 0 0 . It is now obvious that the lemma holds.
To sum up the above discussion, we can draw the following conclusions.
Theorem 5. 
For system (3), when τ 1 = τ 2 = τ , assume that (H1) and (H7) hold, the positive equilibrium E 2 ( x * , y * ) is asymptotically stable for whole τ [ 0 , τ 0 ) and unstable for τ > τ 0 . In addition, a Hopf bifurcation occurs at the positive equilibrium E 2 when τ = τ 0 .
Situation 5: τ 1 [ 0 , τ 10 ) , τ 2 > 0 and τ 1 τ 2 .
Let τ 1 be set to a value selected from its stable interval, and assume that Equation (11) includes only τ 1 as a parameter. Let i ω 2 ( ω 2 > 0 ) be the root of Equation (11), and we obtain
A 51 cos ω 2 τ 2 + A 52 sin ω 2 τ 2 = A 53 , A 51 sin ω 2 τ 2 A 52 cos ω 2 τ 2 = A 54 ,
where
A 51 = E ω 2 G sin ω 2 τ 1 , A 52 = F G cos ω 2 τ 1 , A 53 = A ω 2 C ω 2 cos ω 2 τ 1 + D sin ω 2 τ 1 , A 54 = ω 2 2 B C ω 2 sin ω 2 τ 1 D cos ω 2 τ 1 .
From Equation (24), we have
ω 2 4 + p 51 ω 2 2 + p 52 + p 53 ω 2 2 + p 54 cos ω 2 τ 1 + p 55 ω 2 3 + p 56 ω 2 sin ω 2 τ 1 = 0 ,
with
p 51 = A 2 + C 2 2 B E 2 , p 52 = B 2 + D 2 F 2 G 2 , p 53 = 2 D + 2 A C , p 54 = 2 B D 2 F G , p 55 = 2 C , p 56 = 2 A D + 2 B C + 2 E G .
We make the following assumption. (H8) Equation (25) has at least a finite positive root.
Assume that (H8) holds, we write the positive roots of (25) down as ω 2 ( 1 ) , , ω 2 ( k ) . For every ω 2 ( i ) ( i = 1 , , k ) , the corresponding critical value of time delay τ 2 i ( j ) ( j = 0 , 1 , 2 , ) is shown by
τ 2 i ( j ) = 1 ω 2 i arccos A 51 A 53 A 52 A 54 A 51 2 + A 52 2 + 2 n π ω 2 i , i = 1 , , k ; j = 0 , 1 , 2
Let τ 20 = min τ 2 i ( j ) i = 1 , k ; j = 0 , 1 , 2 , ω 20 is the corresponding root of Equation (25): with τ 20 .
Next, we verify the transversality condition.
Lemma 4. 
If Condition (H9) Q 51 Q 53 + Q 52 Q 54 0 holds, then d ( Re λ ) d τ 2 λ = i ω 20 0 .
Proof. 
Differentiating both sides of Equation (11) about τ 2 , we can derive that
d λ d τ 2 1 = 2 λ + A + C e λ τ 1 τ 1 ( C λ + D ) e λ τ 1 + E e λ τ 2 G τ 1 e λ τ 1 + τ 2 λ λ 2 + A λ + B + ( C λ + D ) e λ τ 1 τ 2 λ .
This, hence, leads to
Re d λ d τ 2 i = i ω 20 1 = Re Q 51 + Q 52 i Q 53 + Q 54 i = Q 51 Q 53 + Q 52 Q 54 Q 53 2 + Q 54 2 ,
where
Q 51 = A + C τ 1 D cos ω 20 τ 1 C ω 20 τ 1 sin ω 20 τ 1 + cos ω 20 τ 20 E G τ 1 cos ω 20 τ 1   + sin ω 20 τ 20 E + G τ 1 sin ω 20 τ 1 , Q 52 = 2 ω 20 + τ 1 D C sin ω 20 τ 1 C ω 20 τ 1 cos ω 20 τ 1 + cos ω 20 τ 20 G τ 1 sin ω 20 τ 1   + sin ω 20 τ 20 G τ 1 cos ω 20 τ 1 , Q 53 = A ω 20 2 + C ω 20 2 cos ω 20 τ 1 D ω 20 sin ω 20 τ 1 , Q 54 = ω 20 3 B ω 20 C ω 20 2 sin ω 20 τ 1 D ω 20 cos ω 20 τ 1 .
Clearly, if (H9) Q 51 Q 53 + Q 52 Q 54 0 is satisfied, we can obtain d ( Re λ ) d τ 2 λ = i ω 20 0 . This completes the proof. □
Here we summarize the results as follows.
Theorem 6. 
For system (3), τ 2 > 0 , τ 1 [ 0 , τ 10 ) and τ 1 τ 2 . Assume that the conditions (H8) and (H9) are satisfied, then the positive equilibrium E 2 ( x * , y * ) is asymptotically stable for whole τ 2 [ 0 , τ 20 ) and unstable for τ 2 > τ 20 . Moreover, a Hopf bifurcation happens at the equilibrium E 2 when τ 2 = τ 20 .
Situation 6: τ 1 > 0 , τ 2 [ 0 , τ 20 ) , and τ 1 τ 2 .
We regard Equation (11) with τ 2 in its stable interval, and τ 1 is taken as the parameter. Like Situation 5, we directly give the following theorem.
Theorem 7. 
For system (3), τ 1 > 0 , τ 2 0 , τ 20 , and τ 1 τ 2 . The positive equilibrium E 2 x * , y * is asymptotically stable for all τ 1 0 , τ 10 and unstable for τ 1 > τ 10 . In addition, a Hopf bifurcation at the positive equilibrium E 2 occurs when τ 1 = τ 10 , where τ 10 denotes the minimum critical value of time delay τ 1 for the appearance of Hopf bifurcation when τ 2 0 , τ 20 .

3. Direction and Stability of Hopf Bifurcation

In the preceding section, we have discussed that system (3) goes through Hopf bifurcation for different combinations of τ 1 and τ 2 . In this section, we will utilize the central manifold theory and the normal form theorem of Hassard et al. [35] to derive the quality analysis of Hopf bifurcation at τ 2 = τ 20 , τ 1 * 0 , τ 10 .
Without loss of generality, we suppose that τ 1 * < τ 20 . Let τ 2 = τ 20 + μ , μ R , t = s τ 2 , x s τ 2 = x ^ ( s ) , y s τ 2 = y ^ ( s ) , and records as x = x ˜ , y = y ˜ , t = s , then system (3) is rewritten by a functional differential equation (FDE) in C ˜ = C ˜ [ 1 , 0 ] , R 2 :
U ( t ) = L μ U t + F μ , U t ,
where
U ( t ) = ( x ( t ) , y ( t ) ) T C ˜ , U t ( θ ) = ( x ( t + θ ) , y ( t + θ ) ) T C ˜ , and   L μ : C ˜ R 2 , F : R × C ˜ R 2
are described as
L μ ( ϕ ) = τ 20 + μ A ϕ ( 0 ) + B ϕ τ 1 * τ 20 + C ϕ ( 1 ) .
and
f ( ξ , ϕ ) = τ 20 + μ f 1 , f 2 T ,
where
ϕ ( θ ) = ϕ 1 ( θ ) , ϕ 2 ( θ ) T C ˜ , A = a 11 a 12 0 a 22 , B = b 11 0 0 0 , C = 0 0 a 21 b 22 , f 1 = k 1 ϕ 1 2 ( 0 ) + k 2 ϕ 1 ( 0 ) ϕ 1 τ 1 * τ 20 + k 3 ϕ 1 ( 0 ) ϕ 2 ( 0 ) + k 4 ϕ 1 2 ( 0 ) ϕ 2 ( 0 ) + k 5 ϕ 2 2 ( 0 ) + k 6 ϕ 2 2 ( 0 ) ϕ 1 ( 0 ) + k 7 ϕ 1 3 ( 0 ) + k 8 ϕ 2 3 ( 0 ) + , f 2 = g 1 ϕ 1 2 ( 1 ) + g 2 ϕ 1 ( 1 ) ϕ 2 ( 1 ) + g 3 ϕ 2 2 ( 0 ) + g 4 ϕ 1 2 ( 1 ) ϕ 2 ( 1 ) + g 5 ϕ 2 2 ( 1 ) + g 6 ϕ 2 2 ( 1 ) ϕ 1 ( 1 ) + g 7 ϕ 1 3 ( 1 ) + g 8 ϕ 2 3 ( 1 ) + , k 1 = β 1 b ( 1 m ) 2 y * a + c y * a + b ( 1 m ) x * + c y * 3 , k 2 = 1 , k 3 = β 1 ( 1 m ) a 2 + a b ( 1 m ) x * + a c y * + 2 b c ( 1 m ) x * y * a + b ( 1 m ) x * + c y * 3 , k 4 = β 1 b ( 1 m ) 2 a 2 + a b ( 1 m ) x * c 2 y * 2 + 2 b c ( 1 m ) x * y * a + b ( 1 m ) x * + c y * 4 , k 5 = β 1 c ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 3 , k 6 = β 1 c ( 1 m ) a + 2 b ( 1 m ) x * a + b ( 1 m ) x * + c y * 3 b ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 4 , k 7 = β 1 b 2 ( 1 m ) 3 a y * + c y * 2 a + b ( 1 m ) x * + c y * 4 , k 8 = β 1 c 2 ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 4 , g 1 = β 2 b ( 1 m ) 2 y * a + c y * a + b ( 1 m ) x * + c y * 3 , g 2 = β 2 ( 1 m ) a 2 + a b ( 1 m ) x * + a c y * + 2 b c ( 1 m ) x * y * a + b ( 1 m ) x * + c y * 3 , g 3 = e , g 4 = β 2 b ( 1 m ) 2 a 2 + a b ( 1 m ) x * c 2 y * 2 + 2 b c ( 1 m ) x * y * a + b ( 1 m ) x * + c y * 4 , g 5 = β 2 c ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 3 , g 6 = β 2 c ( 1 m ) a + 2 b ( 1 m ) x * a + b ( 1 m ) x * + c y * 3 b ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 4 , g 7 = β 2 b 2 ( 1 m ) 3 a y * + c y * 2 a + b ( 1 m ) x * + c y * 4 , g 8 = β 2 c 2 ( 1 m ) a x * + b ( 1 m ) x * 2 a + b ( 1 m ) x * + c y * 4 .
Hence, in light of the Riesz representation theorem, there exists a 2 × 2 matrix function η ( θ , μ ) of bounded variation when θ [ 1 , 0 ] , which implies
L μ ϕ = 1 0 d η ( θ , μ ) ϕ ( θ ) , for ϕ C ˜ .
As a matter of fact, we can choose
η ( θ , μ ) = τ 20 + μ A + B + C , θ = 0 , τ 20 + μ B + C , θ τ 1 * τ 20 , 0 , τ 20 + μ C , θ 1 , τ 1 * τ 20 , 0 , θ = 1 .
For ϕ C ˜ [ 1 , 0 ] , R 2 , define
A ˜ ( μ ) ϕ = d ϕ ( θ ) d θ , 1 θ < 0 , 1 0 d η ( θ , μ ) ϕ ( θ ) , θ = 0 .
and
R ˜ μ ( ϕ ) = 0 , 1 θ < 0 , f ( μ , ϕ ) , θ = 0 .
Then system (27) can take the form
U t = A ˜ ( μ ) U t + R ˜ ( μ ) U t .
For ψ C [ 1 , 0 ] , R 2 * , where R 2 * is the 2-dimensional space of row vectors, the adjoint operator A * of A ˜ ( 0 ) is further defined:
A * ψ ( s ) = d ψ ( s ) d s , s ( 0 , 1 ] , 1 0 d η T ( t , 0 ) ψ ( t ) , s = 0 .
For ϕ C ˜ [ 1 , 0 ] , R 2 and ψ C [ 1 , 0 ] , R 2 * , the bilinear form is then defined:
ψ ( s ) , ϕ ( θ ) = ψ ¯ ( 0 ) ϕ ( 0 ) 1 0 ζ = 0 θ ψ ¯ ( ζ θ ) d η ( θ ) ϕ ( ζ ) d ζ ,
where η ( θ ) = η ( θ , 0 ) , A ˜ = A ˜ ( 0 ) and A * are adjoint operators.
In accordance with the previous section, we observe that ± i ω 20 τ 20 are eigenvalues of A ˜ ( 0 ) . Hence, they are still the eigenvalues of A * . Supposing that q ( θ ) = ( 1 , δ ) T e i ω 20 τ 20 θ is the eigenvector of A ˜ ( 0 ) corresponding to i ω 20 τ 20 and q * ( s ) = M ˜ 1 , δ * e i ω 20 τ 20 s is the eigenvector of A * corresponding to i ω 20 τ 20 , after some algebra, we obtain
δ = i ω 20 a 11 b 11 e i ω 20 τ 1 * a 12 , δ * = i ω 20 + a 11 + b 11 e i ω 20 τ 1 * a 21 e i ω 20 τ 20
By a further calculation, we obtain from Equation (36)
q * ( s ) , q ( θ ) = M ˜ 1 , δ ¯ * ( 1 , δ ) T 1 0 ζ = 0 θ q ¯ * ( ζ θ ) d η ( θ ) q ( ζ ) d ζ = M ˜ 1 + δ δ ¯ * 1 0 1 , δ ¯ * θ e i ω 22 τ 20 θ d η ( θ ) ( 1 , δ ) T = M ˜ 1 + δ δ ¯ * + b 11 τ 1 * e i ω 20 τ 1 * + τ 20 e i ω 20 τ 20 a 21 δ ¯ * + b 22 δ δ ¯ * .
This implies that
M ˜ = 1 + δ δ ¯ * + b 11 τ 1 * e i ω 20 τ 1 * + τ 20 e i ω 20 τ 20 a 21 δ ¯ * + b 22 δ δ ¯ * 1 ,
such that q * ( s ) , q ( θ ) = 1 , q * ( s ) , q ¯ ( θ ) = 0 . According to the computing method in Hassard et al. [35] and the calculation process in [38,39,40], we derive the coefficients to determine the relevant bifurcation properties.
g 20 = 2 M ˜ τ 20 k 1 + k 2 e i ω 20 τ 1 * + k 3 δ + k 5 δ 2 + δ ¯ * g 1 e 2 i ω 20 σ 20 + g 2 δ e 2 i ω 20 τ 20 + g 3 δ 2 + g 5 δ 2 e 2 · ω 20 τ 20 , g 11 = M ˜ τ 20 2 k 1 + k 2 e i ω 20 τ 1 * + e i ω 20 τ 1 * + k 3 ( δ + δ ¯ ) + 2 k 5 δ δ ¯ + δ ¯ * 2 g 1 + g 2 ( δ + δ ¯ ) + 2 g 3 δ δ ¯ + 2 g 5 δ δ ¯ , g 02 = 2 M ˜ τ 20 k 1 + k 2 e i ω 20 τ 1 * + k 3 δ ¯ + k 5 δ ¯ 2 + δ ¯ * g 1 e 2 i ω 20 τ 20 + g 2 δ ¯ 2 i ω 20 τ 20 + g 3 δ ¯ 2 + g 5 δ ¯ 2 e 2 i ω 20 τ 20 , g 21 = 2 M ˜ τ 20 k 1 W 20 ( 1 ) ( 0 ) + 2 W 11 ( 1 ) ( 0 ) + k 2 W 11 ( 1 ) τ 1 * τ 20 + 1 2 W 20 ( 1 ) τ 1 * τ 20 + 1 2 W 20 ( 1 ) ( 0 ) e i ω 20 τ 1 * + W 11 ( 1 ) ( 0 ) e i ω 20 τ 1 * + k 3 W 11 ( 2 ) ( 0 ) + 1 2 W 20 ( 2 ) ( 0 ) + 1 2 W 20 ( 1 ) ( 0 ) δ ¯ + W 11 ( 1 ) ( 0 ) δ + k 4 ( 2 δ + δ ¯ ) + k 5 2 δ W 11 ( 2 ) ( 0 ) + δ ¯ W 20 ( 2 ) ( 0 ) + k 6 δ 2 + 2 δ δ ¯ + 3 k 7 + 3 k 8 δ 2 δ ¯ + δ ¯ * g 1 W 20 ( 1 ) ( 1 ) e i ω 20 τ 20 + 2 W 11 ( 1 ) ( 1 ) e i ω 20 τ 20 + g 2 W 11 ( 2 ) ( 1 ) e i ω 20 τ 20 + 1 2 W 20 ( 2 ) ( 1 ) e i ω 20 τ 20 + 1 2 δ ¯ W 20 ( 1 ) ( 1 ) e i ω 20 τ 20 + δ W 11 ( 1 ) ( 1 ) e i ω 20 τ 20 + g 3 2 δ W 11 ( 2 ) ( 0 ) + δ ¯ W 20 ( 2 ) ( 0 ) + g 4 2 δ e i ω 20 τ 20 + δ ¯ e i ω 20 τ 20 + g 5 2 δ W 11 ( 2 ) ( 1 ) e i ω 20 τ 20 + δ ¯ W 20 ( 2 ) ( 1 ) e i ω 20 τ 20 + g 6 δ 2 e i ω 22 τ 20 + 2 δ δ ¯ e i ω 20 τ 20 + 3 g 7 e i ω 20 τ 20 + 3 g 8 δ 2 δ ¯ e i ω 20 τ 20 .
and
W 20 ( θ ) = i g 20 ω 20 τ 20 q ( 0 ) e i ω 20 τ 20 θ + i g ¯ 02 3 ω 20 τ 20 q ¯ ( 0 ) e i ω 20 τ 20 θ + M 1 e 2 i ω 20 τ 20 θ , W 11 ( θ ) = i g 11 ω 20 τ 20 q ( 0 ) e i ω 20 τ 20 θ + i g ¯ 11 ω 20 τ 20 q ¯ ( 0 ) e i ω 20 τ 20 θ + M 2 .
where M 1 =   ( M 1 ( 1 ) , M 1 ( 2 ) ) T R 2 and M 2 =   ( M 2 ( 1 ) , M 2 ( 2 ) ) T R 2 are also constant vectors: and to determine the values of M 1 and M 2 , we give
2 i ω 20 a 11 b 11 e 2 i ω 20 τ 1 * a 12 a 21 e 2 i ω 20 τ 20 2 i ω 20 a 22 b 22 e 2 i ω 20 τ 20 M 1 = 2 H 1 H 2 , a 11 b 11 a 12 a 21 a 22 b 22 M 2 = H 3 H 4 ,
with
H 1 = k 1 + k 2 e i ω 20 τ 1 * + k 3 δ + k 5 δ 2 H 2 = g 1 e 2 i ω 20 τ 20 + g 2 ε 2 i ω 20 τ 20 + g 3 δ 2 + g 5 δ 2 e 2 i ω 20 τ 20 H 3 = 2 k 1 + k 2 e i ω 20 τ 1 * + e i ω 20 τ 1 * + k 3 ( δ + δ ¯ ) + 2 k 5 δ δ ¯ H 4 = 2 g 1 + g 2 ( δ + δ ¯ ) + 2 g 3 δ δ ¯ + 2 g 5 δ δ ¯
Consequently, the coefficients g 21 can be solved and we can obtain the values:
c ˜ 1 ( 0 ) = i 2 ω 20 τ 20 g 20 g 11 2 g 11 2 g 02 2 3 + g 21 2 , μ ˜ 2 = Re c ˜ 1 ( 0 ) Re λ τ 20 , β ˜ 2 = 2 Re c ˜ 1 ( 0 ) , T ˜ 2 = Im c ˜ 1 ( 0 ) + μ ˜ 2 Im λ τ 20 ω 20 τ 20 ,
which decides the nature of Hopf bifurcation at the critical value τ 2 = τ 20 . Therefore, we have the following results.
Theorem 8. 
For system (3),
(1) 
The sign of μ ˜ 2 determines the direction of Hopf bifurcation: if μ ˜ 2 > 0 , then the Hopf bifurcation is supercritical; if μ ˜ 2 < 0 , then the Hopf bifurcation is subcritical.
(2) 
The sign of β ˜ 2 decides the stability of bifurcating periodic solutions: if β ˜ 2 < 0 , the bifurcating periodic solutions are stable; if β ˜ 2 > 0 , the bifurcating periodic solutions are unstable.
(3) 
The sign of T ˜ 2 determines the period of bifurcating periodic solutions: if T ˜ 2 > 0 , the period increases; if T ˜ 2 < 0 , the period decreases.

4. Numerical Simulations

This section gives numerical examples and simulations to verify the above analytical findings using Matlab software.
We take the system Matlab software to give numerical examples and simulations to verify the above analytical finding s 1 = 2 , c 1 = 1 , α 2 = 2 , d 1 = 0.25 , ρ 2 = 0.05 , m = 0.5 , and then we obtain β 1 = 8 3 , a = 0.1 , b = 2 , c = 1 , β 2 = 2 , m = 0.5 , d = 0.25 , e = 0.01 , whereby the system (3) can be expressed by
x ˙ = x ( t ) 1 x t τ 1 8 3 ( 1 0.5 ) x ( t ) y ( t ) 0.1 + 2 · ( 1 0.5 ) x ( t ) + 1 · y ( t ) , y ˙ = 2 · ( 1 0.5 ) x t τ 2 y t τ 2 0.1 + 2 · ( 1 0.5 ) x t τ 2 + 1 · y t τ 2 0.25 y ( t ) 0.01 y 2 ( t ) .
System (40) exists the equilibrium E 0 ( 0 , 0 ) and E 1 ( 1 , 0 ) , and both are unstable from the above-mentioned Theorems 1 and 2. When the condition (H1) holds, the unique positive equilibrium E 2 ( 0.1872 , 0.4485 ) is obtained.
(i)
For τ 1 = τ 2 = 0 , by calculation, the condition (H3) holds, consequently the equilibrium E 2 is locally asymptotically stable.
(ii)
For τ 1 = 0 , τ 2 > 0 , we can derive that the conditions (H1), (H3), (H5) are satisfied, the corresponding simulation results are ω 20 = 0.1986 , τ 20 = 2.5714 . According to Theorem 3, we know that the positive equilibrium E 2 is asymptotically stable when τ 2 [ 0 , τ 20 ) , and for τ 2 > τ 20 , E 2 is unstable, which are shown in Figure 1 and Figure 2. In the same way, we obtain ω 10 = 0.2980 , τ 10 = 2.8701 (see Figure 3 and Figure 4).
(iii)
For τ 1 = τ 2 = τ 0 , we further obtain ω 0 = 0.2707 , τ 0 = 1.3689 . From Theorem 5, the positive equilibrium E 2 is asymptotically stable as the time delay τ increases from 0 to τ 0 . When the time delay τ goes through the critical value τ 0 , the positive equilibrium E 2 will lose stability, and a Hopf bifurcation appears (see Figure 5 and Figure 6).
(iv)
For τ 2 > 0 , τ 1 = 1 0 , τ 10 , we have ω 20 = 0.2296 , τ 20 = 1.8564 . From Theorem 6, E 2 is asymptotically stable when τ 2 0 , τ 20 and unstable when τ 2 > τ 20 . By the formulas obtained in the above section, we can obtain c ˜ 1 ( 0 ) = 3.7126 + 0.2034 I , β ˜ 2 = 7.4252 , μ ˜ 2 = 159.3391 , T ˜ 2 = 6.7752 . By Theorem 8, the nature of Hopf bifurcation is supercritical; the bifurcating periodic solutions are stable, and its period increases, which can be depicted in Figure 7 and Figure 8. Numerical simulations illustrate our theoretical analysis. Due to the bifurcation of periodic solutions being stable, that is, the number of predators and prey implies stability and coexistence. Similarly, for τ 1 > 0 , τ 2 = 0.5 [ 0 , τ 20 ) , we can have ω 10 = 0.2842 , τ 10 = 2.3894 (see Figure 9 and Figure 10).

5. Conclusions

This paper proposes and analyzes a predator–prey model with two time delays and a prey refuge. As mentioned earlier, the Beddington–DeAngelis type functional response is superior to other types. Then we describe the interaction between prey and predator as a Beddington–DeAngelis type functional response. First, we study the existence of the trivial equilibrium, the predator–extinction equilibrium, the coexistence equilibrium, and the relevant local stability of each feasible equilibrium of the system (3). We found that the equilibrium E 0 is unstable, and if (H2) ( 1 m ) ( β 2 b d ) < a d holds, then the equilibrium E 1 is locally asymptotically stable when τ 2 = 0 ; for whole τ 2 > 0 , the E 1 is unstable. Afterward, we discuss the local stability of a unique positive equilibrium E 2 . According to the existing two time delays and discussing six different cases, we know that the positive equilibrium will lose its original stability, and a Hopf bifurcation will occur. A family of periodic solutions bifurcates from E 2 when the time delay passes through some critical values, which means that the existence of time delay can cause the number of predator and prey populations to change periodically. Therefore, the time delay has a great influence on the dynamics of system (3). For τ 2 > 0 , τ 1 = 1 0 , τ 10 , we obtain ω 20 = 0.2296 , τ 20 = 1.8564 . From Theorem 6, E 2 is asymptotically stable when τ 2 0 , τ 20 and unstable when τ 2 > τ 20 . Moreover, the properties of Hopf bifurcation are analyzed by normal form theory and center manifold theorem. In the end, numerical simulations verify our theoretical analysis; for instance, we obtain c ˜ 1 ( 0 ) = 3.7126 + 0.2034 I , β ˜ 2 = 7.4252 , μ ˜ 2 = 159.3391 , T ˜ 2 = 6.7752 . Based on Theorem 8, the nature of Hopf bifurcation is supercritical, the bifurcating periodic solutions are stable, and its period increases. From a biological point of view, a stable bifurcation periodic solution implies the coexistence of both populations in an oscillatory mode. In addition, the population can be divided into two stages: immature and mature. As for the question of what the system’s dynamic behavior will be if we are concerned with the combined effects of stage structure for predator and prey, we leave it for our future research.

Author Contributions

Methodology, Z.Z. and Y.C.; formal analysis, M.P.; investigation, R.L. and M.M.A.K.; resources, Z.Z.; data curation, M.P. and R.L.; supervision, Y.C. and M.M.A.K.; project administration, M.M.A.K.; funding acquisition, Z.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by the National Natural Science Foundation of China (Grant Nos. 12102148 and 11872189), and the Natural Science Research of Jiangsu Higher Education Institutions of China (Grant No. 21KJB110010).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments

The authors are thankful to the referees for the careful reading and valuable suggestions that improve the quality and description of this manuscript.

Conflicts of Interest

The authors declare that there are no conflict of interest regarding the publication of this paper.

References

  1. Lotka, A.J. Elements of Physical Biology; Williams and Wilkins: Baltimore, MD, USA, 1925. [Google Scholar]
  2. Volterra, V. Variazioni e Fluttuazioni del Numero D’individui in Specie Animali Conviventi; Società Anonima Tipografica “Leonardo da Vinci”: Città di Castello, Italy, 1925. [Google Scholar]
  3. Berryman, A.A. The orgins and evolution of predator-prey theory. Ecology 1992, 73, 1530–1535. [Google Scholar] [CrossRef] [Green Version]
  4. Waezizadeh, T.; Mehrpooya, A. A stochastic model for dynamics of two populations and its stability. In Proceedings of the 47th Annual Iranian Mathematics Conference (AIMC47), Karaj, Iran, 28–31 August 2016; Volume 10681071. [Google Scholar]
  5. Gokila, C.; Sambath, M.; Balachran, K.; Ma, Y.K. Analysis of stochastic predator-prey model with disease in the prey and Holling type II functional response. Adv. Math. Phys. 2020, 2020, 3632091. [Google Scholar] [CrossRef]
  6. Song, Y.; Xiao, W.; Qi, X.Y. Stability and Hopf bifurcation of a predator-prey model with stage structure and time delay for the prey. Nonlinear Dyn. 2016, 83, 1409–1418. [Google Scholar] [CrossRef]
  7. Peng, M.; Zhang, Z.D.; Wang, X.D. Hybrid control of Hopf bifurcation in a Lotka-Volterra predator-prey model with two delays. Adv. Differ. Equ. 2017, 387, 1–20. [Google Scholar] [CrossRef] [Green Version]
  8. Peng, M.; Zhang, Z.D.; Lim, C.W.; Wang, X.D. Hopf bifurcation and hybrid control of a delayed ecoepidemiological model with nonlinear incidence rate and Holling type II functional response. Math. Probl. Eng. 2018, 3, 1–12. [Google Scholar] [CrossRef]
  9. Chen, J.L.; Chen, Y.M.; Zhu, Z.L.; Chen, F.D. Stability and bifurcation of a discrete predator-prey system with Allee effect and other food resource for the predators. J. Appl. Math. Comput. 2022. [Google Scholar] [CrossRef]
  10. Yao, P.; Wang, Z.C.; Wang, L.S. Stability Analysis of a Ratio-Dependent Predator-Prey Model. J. Math. 2022, 2022, 4605267. [Google Scholar] [CrossRef]
  11. Holling, C.S. The functional response of predator to prey density and its role in mimicry and population regulation. Mem. Entomol. Soc. Can. 1959, 91, 385–398. [Google Scholar] [CrossRef]
  12. Liu, B.; Chen, L.S.; Zhang, Y.J. The dynamics of a prey-dependent consumption model concerning impulsive control strategy. Appl. Math. Comput. 2005, 169, 305–320. [Google Scholar] [CrossRef]
  13. Beddington, J.R. Mutual interference between parasites or predators and its effects on searching effiency. J. Anim. Ecol. 1975, 44, 331–340. [Google Scholar] [CrossRef]
  14. De-Angelis, D.L.; Goldstein, R.A.; O’Neill, R.V. A model for tropic interaction. Ecology 1975, 56, 881–892. [Google Scholar] [CrossRef]
  15. Hwang, T.W. Uniqueness of limit cycles of the predator-prey system with Beddington-Deangelis functional response. J. Math. Anal. Appl. 2004, 290, 113–122. [Google Scholar] [CrossRef] [Green Version]
  16. Li, H.Y.; Takeuchi, Y. Dynamics of the density dependent predator-prey system with Beddington-DeAngelis functional response. J. Math. Anal. Appl. 2011, 374, 644–654. [Google Scholar] [CrossRef] [Green Version]
  17. Tripathi, J.P.; Jana, D.; Tiwari, V. A Beddington-DeAngelis type one-predator two-prey competitive system with help. Nonlinear Dyn. 2018, 94, 553–573. [Google Scholar] [CrossRef]
  18. Pelen, N.N. On the Dynamics of Impulsive Predator-Prey Systems with Beddington- Deangelis-Type Functional Response. Ukr. Math. J. 2021, 73, 610–634. [Google Scholar] [CrossRef]
  19. Feng, X.Z.; Li, C.T.; Sun, H.; Wang, Y.Z. Global Bifurcation Structure of a Predator-Prey System with a Spatial Degeneracy and B-D Functional Response. Complexity 2021, 2021, 9970255. [Google Scholar] [CrossRef]
  20. Tang, G.Y.; Tang, S.Y.; Cheke, R.A. Global analysis of a Holling type II predator-prey model with a constant prey refuge. Nonlinear Dyn. 2014, 76, 635–647. [Google Scholar] [CrossRef]
  21. Khajanchi, S.; Banerjee, S. Role of constant prey refuge on stage structure predator-prey model with ratio dependent functional response. Appl. Math. Comput. 2017, 314, 193–198. [Google Scholar] [CrossRef]
  22. Peng, M.; Zhang, Z.D.; Wang, X.D.; Liu, X.Y. Hopf bifurcation analysis for a delayed predator-prey system with a prey refuge and selective harvesting. J. Appl. Anal. Comput. 2018, 8, 982–997. [Google Scholar]
  23. Zhang, H.S.; Cai, Y.L.; Fu, S.M.; Wang, W.M. Impact of the fear effect in a prey-predator model incorporating a prey refuge. Appl. Math. Comput. 2019, 356, 328–337. [Google Scholar] [CrossRef]
  24. Qi, H.K.; Meng, X.Z. Threshold behavior of a stochastic predator-prey system with prey refuge and fear effect. Appl. Math. Lett. 2021, 113, 106846. [Google Scholar] [CrossRef]
  25. Yang, R.Z.; Liu, M.; Zhang, C.R. A delayed-diffusive predator-prey model with a ratio-dependent functional response. Commun. Nonlinear Sci. Numer. Simul. 2017, 53, 94–110. [Google Scholar] [CrossRef]
  26. Li, L.Y.; Mei, Y.Y.; Cao, J.D. Hopf bifurcation analysis and stability for a ratio-dependent predator-prey diffusive system with time delay. Int. J. Bifurc. Chaos 2020, 30, 2050037. [Google Scholar] [CrossRef]
  27. Ma, Z.H.; Wang, S.F. A delay-induced predator-prey model with Holling type functional response and habitat complexity. Nonlinear Dyn. 2018, 93, 1519–1544. [Google Scholar] [CrossRef]
  28. Xiao, Z.W.; Xie, Z.D.; Xue, Y.L. Stability and bifurcation in a Holling type II predator-prey model with Allee effect and time delay. Adv. Differ. Equ. 2018, 288, 1–21. [Google Scholar] [CrossRef]
  29. Zheng, T.; Zhang, L.; Luo, Y.T.; Zhou, X.R.; Li, H.L.; Teng, Z.D. Stability and Hopf Bifurcation of a Stage-Structured Cannibalism Model with Two Delays. Int. J. Bifurc. Chaos 2021, 31, 2150242. [Google Scholar] [CrossRef]
  30. Wang, X.D.; Peng, M.; Liu, X.Y. Stability and Hopf bifurcation analysis of a ratio-dependent predator-prey model with two time delays and Holling type III functional response. Appl. Math. Comput. 2015, 268, 496–508. [Google Scholar] [CrossRef]
  31. Du, Y.F.; Niu, B.; Wei, J.J. Two delays induce Hopf bifurcation and double Hopf bifurcation in a diffusive Leslie-Gower predator-prey system. Chaos 2019, 29, 013101. [Google Scholar] [CrossRef] [Green Version]
  32. Brikhoff, G.; Rota, G.C. Ordinary Differential Equations; John Wiley and Sons: New York, NY, USA, 1982. [Google Scholar]
  33. Hale, J.K. Theory of Functional Differential Equations; Springer: New York, NY, USA, 1977. [Google Scholar]
  34. Ruan, S.G.; Wei, J.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]
  35. Hassard, B.D.; Kazarinoff, N.D.; Wan, Y.H. Theory and Application of Hopf Bifurcation; Cambridge University Press: Cambridge, UK, 1981. [Google Scholar]
  36. Kajiwara, T.; Sasaki, T.; Takeuchi, Y. Construction of lyapunov functionals for delay differential equations in virology and epidemiology. Nonlinear Anal. Real World Appl. 2012, 13, 1802–1826. [Google Scholar] [CrossRef]
  37. Manna, K.; Chakrabarty, S.P. Global stability of one and two discrete delay models for chronic hepatitis b infection with HBV DNA-containing capsids. Comput. Appl. Math. 2017, 36, 525–536. [Google Scholar] [CrossRef]
  38. Chen, X.X.; Wang, X.D. Qualitative analysis and control for predator-prey delays system. Chaos Solitons Fractals 2019, 123, 361–372. [Google Scholar] [CrossRef]
  39. Peng, M.; Zhang, Z.D.; Qu, Z.F.; Bi, Q.S. Qualitative analysis in a delayed Van del Pol oscillator. Phys. A Stat. Mech. Its Appl. 2019, 544, 12348. [Google Scholar]
  40. Zhu, L.H.; Wang, X.W.; Zhang, Z.D.; Shen, S.L. Global stability and bifurcation analysis of a rumor propagation model with two discrete delays in social networks. Int. J. Bifurc. Chaos 2020, 30, 2050175. [Google Scholar] [CrossRef]
Figure 1. Behavior and phase portrait of the system (40) when τ 2 = 2.0 < τ 20 = 2.5714 , E 2 is asymptotically stable.
Figure 1. Behavior and phase portrait of the system (40) when τ 2 = 2.0 < τ 20 = 2.5714 , E 2 is asymptotically stable.
Symmetry 14 02535 g001
Figure 2. Behavior and phase portrait of the system (40) when τ 2 = 2.7 > τ 20 = 2.5714 , E 2 is unstable.
Figure 2. Behavior and phase portrait of the system (40) when τ 2 = 2.7 > τ 20 = 2.5714 , E 2 is unstable.
Symmetry 14 02535 g002
Figure 3. Behavior and phase portrait of the system (40) when τ 1 = 2.7 < τ 10 = 2.8701 , E 2 is asymptotically stable.
Figure 3. Behavior and phase portrait of the system (40) when τ 1 = 2.7 < τ 10 = 2.8701 , E 2 is asymptotically stable.
Symmetry 14 02535 g003
Figure 4. Behavior and phase portrait of the system (40) when τ 1 = 2.9 > τ 10 = 2.8701 , E 2 is unstable.
Figure 4. Behavior and phase portrait of the system (40) when τ 1 = 2.9 > τ 10 = 2.8701 , E 2 is unstable.
Symmetry 14 02535 g004
Figure 5. Behavior and phase portrait of the system (40) when τ = 1.0 < τ 0 = 1.3689 , E 2 is asymptotically stable.
Figure 5. Behavior and phase portrait of the system (40) when τ = 1.0 < τ 0 = 1.3689 , E 2 is asymptotically stable.
Symmetry 14 02535 g005
Figure 6. Behavior and phase portrait of the system (40) when τ = 1.5 > τ 0 = 1.3689 , E 2 is unstable.
Figure 6. Behavior and phase portrait of the system (40) when τ = 1.5 > τ 0 = 1.3689 , E 2 is unstable.
Symmetry 14 02535 g006
Figure 7. E 2 is asymptotically stable for τ 2 = 1.1 < τ 20 = 1.8564 and τ 1 = 1.0 .
Figure 7. E 2 is asymptotically stable for τ 2 = 1.1 < τ 20 = 1.8564 and τ 1 = 1.0 .
Symmetry 14 02535 g007
Figure 8. E 2 is unstable for τ 2 = 1.9 > τ 20 = 1.8564 and τ 1 = 1.0 .
Figure 8. E 2 is unstable for τ 2 = 1.9 > τ 20 = 1.8564 and τ 1 = 1.0 .
Symmetry 14 02535 g008
Figure 9. E 2 is asymptotically stable for τ 1 = 1.5 < τ 10 = 2.3894 and τ 2 = 0.5 .
Figure 9. E 2 is asymptotically stable for τ 1 = 1.5 < τ 10 = 2.3894 and τ 2 = 0.5 .
Symmetry 14 02535 g009
Figure 10. E 2 is unstable for τ 1 = 2.4 > τ 10 = 2.3894 and τ 2 = 0.5 .
Figure 10. E 2 is unstable for τ 1 = 2.4 > τ 10 = 2.3894 and τ 2 = 0.5 .
Symmetry 14 02535 g010
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Peng, M.; Lin, R.; Chen, Y.; Zhang, Z.; Khater, M.M.A. Qualitative Analysis in a Beddington–DeAngelis Type Predator–Prey Model with Two Time Delays. Symmetry 2022, 14, 2535. https://doi.org/10.3390/sym14122535

AMA Style

Peng M, Lin R, Chen Y, Zhang Z, Khater MMA. Qualitative Analysis in a Beddington–DeAngelis Type Predator–Prey Model with Two Time Delays. Symmetry. 2022; 14(12):2535. https://doi.org/10.3390/sym14122535

Chicago/Turabian Style

Peng, Miao, Rui Lin, Yue Chen, Zhengdi Zhang, and Mostafa M. A. Khater. 2022. "Qualitative Analysis in a Beddington–DeAngelis Type Predator–Prey Model with Two Time Delays" Symmetry 14, no. 12: 2535. https://doi.org/10.3390/sym14122535

APA Style

Peng, M., Lin, R., Chen, Y., Zhang, Z., & Khater, M. M. A. (2022). Qualitative Analysis in a Beddington–DeAngelis Type Predator–Prey Model with Two Time Delays. Symmetry, 14(12), 2535. https://doi.org/10.3390/sym14122535

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