Next Article in Journal
Extended Gevrey Regularity via Weight Matrices
Previous Article in Journal
A Reliable Technique for Solving Fractional Partial Differential Equation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Stability and Hopf Bifurcation Analysis of a Stage-Structured Predator–Prey Model with Delay

School of Mathematics and Statistics, Xinyang Normal University, Xinyang 464000, China
Axioms 2022, 11(10), 575; https://doi.org/10.3390/axioms11100575
Submission received: 1 September 2022 / Revised: 26 September 2022 / Accepted: 17 October 2022 / Published: 20 October 2022

Abstract

:
In this work, a Lotka–Volterra type predator–prey system with time delay and stage structure for the predators is proposed and analyzed. By using the permanence theory for infinite dimensional system, we get that the system is permanent if some conditions are satisfied. The local and global stability of the positive equilibrium is presented. The existence of Hopf bifurcation around the positive equilibrium is observed. Further, by using the normal form theory and center manifold approach, we derive the explicit formulas determining the stability of bifurcating periodic solutions and the direction of Hopf bifurcation. Numerical simulations are carried out by Matlab software to explain the theoretical results. We find that combined time delay and stage structure can affect the dynamical behavior of the system.

1. Introduction

Differential equations are a powerful tool for characterizing natural phenomena [1,2]. The predator–prey model is a very classic model, which plays a key role in population ecology. Many predator–prey models have been investigated by some researchers [3,4,5,6,7,8,9]. From [10], we know that if the following classical autonomous Lotka–Volterra type predator prey model,
d x ( t ) d t = x ( t ) [ k x ( t ) α y ( t ) ] , d y ( t ) d t = y ( t ) [ r + β x ( t ) y ( t ) ] ,
exists in a positive equilibrium ( x 1 , y 1 ) , it must be globally asymptotically stable. Time lag is pervasive in nature. The stability issues for the Lotka-Volterra system with different types of time delays have been extensively studied. In [11], by using Lyapunov functional, He examined the global attraction for a kind of delayed n-species Lotka–Volterra-type system. In [12], Gopalsamy et al. examined the global stability of a delay nonautonomous n-species competition system. In [13], He investigated the global asymptotic stability of a nonautonomous Lotka–Volterra system with “pure-delay type”. In [14], Wang et al. proved that delays are harmless for the two-dimensional delayed Lotka–Volterra system. As a special case of Lotka–Volterra-type systems with delays, Chen et al. proposed a model of two species’ growth delays as a reasonable generalization of the Lotka–Volterra model, which takes the form [10]:
d x ( t ) d t = x ( t ) [ k x ( t ) α y ( t τ ) ] , d y ( t ) d t = y ( t ) [ r + β x ( t τ ) y ( t ) ] .
System (2) is one of the simplest predator–prey models with a delay. Its stability and Hopf bifurcations, both local and global, have been widely investigated. For example, Wang et al. [14] found that system (2) was uniformly persistent irrespective of the size of the delays. He [11] showed that the positive equilibrium is locally and globally asymptotically stable.
In the real world, however, many consumer species may go though multiple life phases as they progress from birth to death. In [15,16], the authors studied the delayed stage structure models. In those models, a constant time lag represented the time from birth to maturity. References [17,18] have examined the stage structure of species when the transformation rate of the mature population is proportional to the existing immature population. Motivated by the works of Chen [10], He et al. [11,12,13], Cui et al. [15] and Song et al. [16], we built a predator–prey model based on system (2), which includes a time delay due to negative feedback of prey and the stage structure for the predators. This paper’s purpose is to explore the combined effects of both delay and stage structure on the predator–prey system’s dynamical behavior.

2. The Model

We consider a delayed predator–prey system with a stage structure among predator populations of the following form:
d x ( t ) d t = x ( t ) [ k x ( t ) α y ( t τ ) ] , d y ( t ) d t = y ( t ) [ r + β x ( t τ ) y ( t ) ] + σ y i ( t ) , d y i ( t ) d t = c y ( t ) ( σ + γ ) y i ( t ) ,
where x ( t ) expresses prey density at time t, and y ( t ) and y i ( t ) represent densities of the mature and the immature predator species at time t. In model (3), α and β represent the capturing rates of the predators; k is the intrinsic rate of increase for the prey; r represents the mature predator’s death rate; σ represents the conversion rate; γ is the immature predator’s death rate; c denotes the birth rate of the immature predators; τ is a constant delay. All the parameters (i.e., k, α , r, σ , β , c, γ , and τ ) are positive constants.
The initial conditions for system (3) have the following form:
x ( ς ) = φ 1 ( ς ) 0 , y ( ς ) = φ 2 ( ς ) 0 , y i ( ς ) = φ 2 ( ς ) 0 , ς [ τ , 0 ] , φ 1 ( 0 ) > 0 , φ 2 ( 0 ) > 0 .
We suppose that z ( t ) = ( x ( t ) , y ( t ) , y i ( t ) ) is the solution of system (3) with the initial conditions (4). Obviously, under the initial conditions given in (4), the solution z ( t ) of system (3) exists in the interval [ 0 , + ) . Further, it remains positive for all t 0 . In fact, from the 1st equation of (3), we obtain
x ( t ) = x ( 0 ) exp ( 0 t ( k x ( s ) α y ( s τ ) ) d s ) > 0
for t 0 . Set Ξ ( t ) = r + β x ( t τ ) y ( t ) , we can rewrite the last two equations of (3) as
d y ( t ) d t = Ξ ( t ) y ( t ) + σ y i ( t ) , d y i ( t ) d t = c y ( t ) ( σ + γ ) y i ( t ) .
Obviously, there is a unique solution z ( t ) of system (3) in a maximal interval [ 0 , d ) [19]. We can prove that the interval is [ 0 , + ) . Since (5) is a quasimonotone system, ( 0 , 0 ) is a subsolution and w ( e a t , e a t ) with w > y ( 0 ) , y i ( 0 ) and a > σ , σ is a supersolution [19]. This shows that y ( t ) and y i ( t ) are bounded in [ 0 , d ) and hence exist for all t > 0 . Suppose y i ( b ) = 0 . We can obtain d y i ( t ) d t | t = b = 0 and y ( b ) = 0 . This shows that ( y ( t ) , y i ( t ) ) is a solution of (5) at t = b and hence it is zero in ( 0 , + ) . This is a contradiction. Hence, y ( t ) > 0 , y i ( t ) > 0 for t > 0 . Thereby, x ( t ) > 0 for t > 0 .
It is straightforward to see that in system (3) there exist four equilibrium points
E 0 ( 0 , 0 , 0 ) , E 1 ( k , 0 , 0 ) , E 2 ( 0 , y 0 , y i 0 ) , E * ( x * , y * , y i * ) ,
where
y 0 = r ( σ + γ ) + σ c σ + γ , y i 0 = c ( r ( σ + γ ) + σ c ) ( σ + γ ) 2 , x * = ( k + α r ) ( σ + γ ) σ c α ( 1 + α β ) ( σ + γ ) , y * = ( r + k β ) ( σ + γ ) + σ c ( 1 + α β ) ( σ + γ ) , y i * = c ( ( r + β k ) ( σ + γ ) + σ c ) ( 1 + α β ) ( σ + γ ) 2 .
If the conditions
( k + α r ) ( σ + γ ) c σ α > 0
and
r ( σ + γ ) + c σ > 0
hold, all the equilibria are nonnegative.

3. Permanence of System (3)

We first introduce the definition of permanence.
Definition 1
([20]). System (3) is permanent if there exist positive constants M 1 and M 2 such that any positive solution ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) with initial conditions (4) satisfies
M 1 lim inf t + x ( t ) lim sup t + x ( t ) M 2 ,
M 1 lim inf t + y ( t ) lim sup t + y ( t ) M 2 ,
M 1 lim inf t + y i ( t ) lim sup t + y i ( t ) M 2 .
Theorem 1.
Any positive solution ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) with initial condition (4) satisfies 0 < x ( t ) , y ( t ) , y i ( t ) M for t 0 .
Proof. 
Suppose ( x ( t ) , y ( t ) , y i ( t ) ) is a positive solution of system (3) with initial condition (4). According to the first equation of system (3), it follows from the positivity of the solution that
x ˙ ( t ) x ( t ) ( k x ( t ) ) .
The solution of the auxiliary equation
u ˙ ( t ) = u ( t ) ( k u ( t ) ) , u ( 0 ) = x ( 0 ) > 0
has the following properties: there exist ε > 0 and T > 0 such that u ( t ) < k + ε for t T . Hence, by comparing the theorems for ordinary differential equations, we have x ( t ) u ( t ) < k + ε for t T .
Denote M 1 = k + ε . Then,
x ( t ) M 1 .
It follows from (8) and the second equation of system (3) that for t T ,
y ˙ ( t ) y ( r + β M 1 y ) + σ y i .
We define
V 1 ( t ) = ( σ + γ ) y + σ y i .
Along the last two equations of system (3), we calculate the upper-right derivative of V 1 :
D + V 1 ( t ) = ( σ + γ ) y ˙ + σ y i ˙ ( σ + γ ) y ( r + β M 1 y ) + σ c y [ ( σ + γ ) ( r + β M 1 ) + σ c ] 2 4 ( σ + γ ) .
Then, there exist two positive constants M 2 , M 3 , such that
0 < y ( t ) M 2 , 0 < y i ( t ) M 3 .
Denote M = min { M 1 , M 2 , M 3 } . Then, 0 < x ( t ) , y ( t ) , y i ( t ) M for t 0 .
Next, we discuss the permanence of system (3).
To prove that system (3) is persistent, we will use Theorem 4.1 in [20]. □
Theorem 2.
If condition (6) and
r + k β > 0
hold, then system (3) is permanent.
Proof. 
We begin by showing that sets X 1 , X 2 , X 3 and X 4 repel the positive solution of system (3) uniformly. Let us define
C 1 = { ( φ 1 , φ 2 , φ 3 ) C ( [ τ , 0 ] , R + 3 ) : φ 2 ( θ ) = 0 } , C 2 = { ( φ 1 , φ 2 , φ 3 ) C ( [ τ , 0 ] , R + 3 ) : φ 1 ( θ ) = 0 } , X 1 = { ( x , y , y i ) R 3 x > 0 , y = 0 , y i > 0 } , X 2 = { ( x , y , y i ) R 3 x = 0 , y > 0 , y i > 0 } , X 3 = { ( x , y , y i ) R 3 x = 0 , y > 0 , y i = 0 } , X 4 = { ( x , y , y i ) R 3 x = 0 , y = 0 , y i > 0 } , C 0 = C 1 C 2 , C 0 = Int C ( [ τ , 0 ] , R + 3 ) , X 0 = X 1 X 2 X 3 X 4 , X 0 = { ( x , y , y i ) R 3 x > 0 , y > 0 , y i > 0 } .
This choice meets the conditions in Theorem 4.1 in [20]. It suffices to show that, for any solution z ( t ) = ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) initiating from C 0 , there exists an ε 0 > 0 such that lim inf t + d ( z ( t ) , X 0 ) ε 0 . To this end, we shall verify that all the conditions of Theorem 4.1 in [20] are satisfied. It is easy to see that C 0 and C 0 are positively invariant. Obviously, conditions ( i ) and ( i i ) of Theorem 4.1 in [20] are satisfied. In the following, we shall only need to validate conditions ( i i i ) and ( i v ) .
There are three constant solutions, E 0 ( 0 , 0 , 0 ) , E 1 ( k , 0 , 0 ) , E 2 ( 0 , y 0 , y i 0 ) , in X 0 .
In the set X 3 or X 4 , system (3) becomes
y ˙ = y ( r y )   or   y ˙ i = ( σ + γ ) y i .
Clearly, lim t + y ( t ) = 0 , lim t + y i ( t ) 0 . Hence, for any solution ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) initiating from C 0 , we obtain ( x ( t ) , y ( t ) , y i ( t ) ) E 0 ( 0 , 0 , 0 ) as t + .
In the set X 1 , system (3) becomes
x ˙ ( t ) = x ( t ) ( k x ( t ) ) , y i ˙ ( t ) = ( σ + γ ) y i ( t ) .
Obviously, E 1 ( k , 0 ) is globally asymptotically stable. Therefore, any solution ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) initiating from C 1 is such that ( x ( t ) , y ( t ) , y i ( t ) ) E 1 ( k , 0 , 0 ) as t + .
In the set X 2 , system (3) becomes
y ˙ ( t ) = y ( t ) ( r y ( t ) ) + σ y i ( t ) , y i ˙ ( t ) = c y ( t ) ( σ + γ ) y i ( t ) .
Clearly, E 2 ( y 0 , y i 0 ) is globally asymptotically stable. Hence, for any solution ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) which initiates from C 2 , we have ( x ( t ) , y ( t ) , y i ( t ) ) E 2 ( 0 , y 0 , y i 0 ) as t + .
Obviously, E 0 , E 1 , E 2 are isolated invariant, and { E 0 , E 1 , E 2 } is isolated and is an acyclic covering.
It is obvious that W s ( E 0 ) X 0 = . Next, we will show that W s ( E 1 ) X 0 = , W s ( E 2 ) X 0 = .
Assume W s ( E q ) X 0 . Then, in system (3) there exists a positive solution ( x ¯ ( t ) , y ¯ ( t ) , y ¯ i ( t ) ) such that ( x ¯ ( t ) , y ¯ ( t ) , y ¯ i ( t ) )   E 1 ( k , 0 , 0 ) as t + . Hence, we have lim t + x ( t ) = k , that is, for ε > 0 , there exists t 0 > 0 such that x ¯ ( t ) > k ε . From system (3), for t > t 0 + τ ,
y ˙ ( t ) y ( t ) ( r + β ( k ε ) y ( t ) ) .
From (10), we can find that y ¯ ( t ) r + k β . This is a contradiction. Therefore, we have W s ( E 1 ) X 0 = .
Assume W s ( E 2 ) X 0 . Then, there is a positive solution ( x ¯ ( t ) , y ¯ ( t ) , y ¯ i ( t ) ) of (3) that ( x ¯ ( t ) , y ¯ ( t ) , y ¯ i ( t ) ) E 2 ( 0 , y 0 , y i 0 ) as t + . Hence, y y 0 , y y i 0 as t + . From system (3), for t > t 0 + τ , we have
x ˙ ( t ) x ( t ) ( k x ( t ) α ( y 0 + ε ) ) .
From (6), we find that x ¯ ( t ) k α y 0 as t + . This is a contradiction. Therefore, we have W s ( E 2 ) X 0 = .
At this time, we can conclude that X 0 repels the positive solutions ( x ( t ) , y ( t ) , y i ( t ) ) of system (3) which initiates from C 0 uniformly.
Thus, there exists an ε 0 > 0 such that
lim inf t + x ( t ) ε 0 , lim inf t + y ( t ) ε 0 .
From the 3rd equation of system (3), we obtain
y ˙ i ( t ) c ε 0 ( γ + σ ) y i ( t ) .
Then y i c ε 0 γ + σ as t + . Denote ε 1 = c ε 0 γ + σ , then y i ε 1 . □

4. Stability of Equilibria

In this section, we will discuss the sufficient conditions for the stability of all the equilibrium points for system (3).
Firstly, we analyze the local stability of the equilibria E 0 , E 1 , E 2 .
The characteristic equation of equilibrium E 0 is
( λ k ) ( λ 2 + ( σ + γ + r ) λ + r ( σ + γ ) σ c ) = 0 .
Obviously, the above equation always has a positive eigenvalue λ = k . Hence, E 0 ( 0 , 0 , 0 ) is unstable.
The characteristic equation of equilibrium E 1 is
( λ + k ) [ λ 2 + ( r β k + σ + γ ) λ + ( r β k ) ( σ + γ ) σ c ] = 0 .
Let
f ( λ ) = λ 2 + ( r β k + σ + γ ) λ + ( r β k ) ( σ + γ ) σ c .
From (7), we can find ( r β k ) ( σ + γ ) σ c < 0 . Then, f ( + ) = + , f ( 0 ) < 0 . So (11) has a positive root, i.e., the equilibrium E 1 ( k , 0 , 0 ) is unstable.
The characteristic equation of equilibrium E 2 is
( λ k + α y 0 ) [ λ 2 + ( r + 2 y 0 + σ + γ ) λ + ( r + 2 y 0 ) ( σ + γ ) σ c ] = 0 .
Denote
f ( λ ) = λ 2 + ( r + 2 y 0 + σ + γ ) λ + ( r + 2 y 0 ) ( σ + γ ) σ c .
Then, f ( + ) = + , f ( 0 ) < 0 , so (13) has a positive root. Then the equilibrium E 2 is unstable.
In the following, we shall discuss the local and global properties of positive equilibrium E * .
Theorem 3.
The interior equilibrium E * ( x * , y * , y i * ) of system (3) is locally asymptotically stable if conditions (6) and (7) and the condition
α β < 1
hold.
Proof. 
Linearizing system (3) at the equilibrium E * ( x * , y * , y i * ) leads to
x ˙ ( t ) = x * x ( t ) α x * y ( t τ ) , y ˙ ( t ) = β y * x ( t τ ) + ( r + β x * 2 y * ) y ( t ) + σ y i ( t ) , y i ˙ ( t ) = c y ( t ) ( σ + γ ) y i ( t ) .
We obtain the characteristic equation of the form
Γ ( λ , τ ) = λ 3 + p λ 2 + q λ + s + e 2 λ τ ( l λ + n ) = 0 ,
where
p = σ + γ + c σ σ + γ + x * + y * > 0 , q = ( x * + y * ) ( σ + γ ) + x * ( c σ σ + γ + y * ) > 0 , s = x * y * ( σ + γ ) > 0 , l = x * y * α β > 0 , n = x * y * α β ( σ + γ ) > 0 .
When τ = 0 , we can easily check that all the characteristic roots have negative real parts. We will show that all the characteristic roots have negative real parts for all τ > 0 , which implies that E * ( x * , y * , y i * ) is locally asymptotically stable for (3). Obviously, the characteristic Equation (15) has no positive real parts roots. Now, we suppose that there exists a characteristic root of (15) on the imaginary axis of the complex plane for some τ 0 > 0 . Let λ = i ω be such a characteristic root. It is straightforward to see that ω 0 . Substituting ( λ , τ ) = ( i ω , τ 0 ) into (15) and separating the real and imaginary parts, we obtain
p ω 2 + s + n cos ( 2 ω τ 0 ) + l ω sin ( 2 ω τ 0 ) = 0 ,
and
ω 3 + q ω + ω l cos ( 2 ω τ 0 ) n sin ( 2 ω τ 0 ) = 0 .
Furthermore,
( l 2 ω 2 + n 2 ) sin ( 2 ω τ 0 ) = p l ω 3 l s ω n ω 3 + q ω n , ( l 2 ω 2 + n 2 ) cos ( 2 ω τ 0 ) = l ω 4 q l ω 2 + p ω 2 n n s .
Squaring and adding the two equations yields
( l 2 ω 2 + n 2 ) 2 = ( p l ω 3 l s ω n ω 3 + q ω n ) 2 + ( l ω 4 q l ω 2 + p ω 2 n n s ) 2 .
Let Ω = ω 2 and
f ( Ω ) = Ω 3 + ( p 2 2 q ) Ω 2 + ( q 2 2 p s l 2 ) Ω + s 2 n 2 = 0 .
Then, f ( Ω ) must have a positive zero Ω = ω 2 because (17) and ω 0 . What is more, the coefficients of Ω 2 , Ω and f ( 0 ) in (18) need to be positive. In fact, the coefficients of Ω 2 and Ω are expressed in the following ways:
p 2 2 q = ( σ + γ ) 2 + ( c σ σ + γ ) 2 + x * 2 + y * 2 + 2 c σ + 2 c σ σ + γ y * > 0 , q 2 2 p s l 2 = x * 2 ( σ + γ ) 2 + y * 2 ( σ + γ ) 2 + ( c σ x * σ + γ ) 2 + 2 x * 2 c σ + 2 c σ x * 2 y * σ + γ + ( 1 α 2 β 2 ) x * 2 y * 2 > 0 .
Furthermore,
f ( 0 ) = s 2 n 2 = ( σ + γ ) 2 x * 2 y * 2 ( 1 α 2 β 2 ) > 0 .
Hence, f ( Ω ) > 0 . Therefore, f ( Ω ) = 0 has no positive roots, which is a contradiction. We complete the proof. □
Theorem 4.
Suppose that (6) and (7) and
α + β < 2
hold. Then the positive equilibrium E * ( x * , y * , y i * ) of system (3) is globally asymptotically stable.
Proof. 
Let ( x ( t ) , y ( t ) , y i ( t ) ) be any solution of system (3) satisfying initial condition (4). Define
V ( t ) = λ 1 ( x x * x * ln x x * ) + λ 2 ( y y * y * ln y y * ) + λ 3 ( y i y i * y i * ln y i y i * ) + λ 2 β 2 t τ t ( x * x ( s ) ) 2 d s + λ 1 α 2 t τ t ( y * y ( s ) ) 2 d s ,
where λ 1 , λ 2 , λ 3 are suitable positive constants to be determined in the subsequent steps. It is easy to see that V is a positive definite function in the region R + 3 except at E * ( x * , y * , y i * ) where it vanishes. Further,
lim x 0 , y 0 , y i 0 V ( x , y , y i ) = lim x , y , y i V ( x , y , y i ) = .
Calculating the upper right derivative of V along the solutions of system (3), we have:
D + V ( t ) | ( 3 ) = λ 1 ( x x * ) ( k x α y ( t τ ) ) + λ 2 ( y y * ) [ ( r + β x ( t τ ) y ) + σ y i y ] + λ 3 ( y i y i * ) ( c y y i ( σ + γ ) ) + λ 2 β 2 [ ( x x * ) 2 ( x ( t τ ) x * ) 2 ] + λ 1 α 2 [ ( y * y ) 2 ( y ( t τ ) y * ) 2 ] = λ 1 ( x x * ) 2 + λ 1 ( x x * ) ( α y * α y ( t τ ) ) λ 2 ( y y * ) 2 + λ 2 ( y y * ) ( β x ( t τ ) β x * c σ σ + γ + σ y i y ) + λ 3 ( y i y i * ) ( c y y i ( σ + γ ) ) + λ 2 β 2 [ ( x x * ) 2 ( x ( t τ ) x * ) 2 ] + λ 1 α 2 [ ( y * y ) 2 ( y ( t τ ) y * ) 2 ] λ 1 ( x x * ) 2 + λ 1 α 2 ( x x * ) 2 λ 2 ( y y * ) 2 + λ 2 β 2 ( y y * ) 2 + λ 2 β 2 ( x x * ) 2 + λ 1 α 2 ( y y * ) 2 + λ 2 ( y y * ) ( σ + γ ) y [ ( σ + γ ) σ y i c σ y ] + λ 3 ( y i y i * ) y i ( c y ( σ + γ ) y ) = ( λ 1 + λ 1 α 2 + λ 1 β 2 ) ( x x * ) 2 + ( λ 2 + λ 2 α 2 + λ 2 β 2 ) ( y y * ) 2 [ λ 2 σ y * y i λ 3 ( σ + γ ) y y i * ( σ + γ ) y y i + λ 3 ( σ + γ ) y y i λ 2 σ y y i ( σ + γ ) y y i ] ( c y ( σ + γ ) y i ) .
Let λ 2 σ = λ 3 ( σ + γ ) . Then,
D + V ( t ) | ( 3 ) ( λ 1 + λ 1 α 2 + λ 2 β 2 ) ( x x * ) 2 + ( λ 2 + λ 1 α 2 + λ 2 β 2 ) ( y y * ) 2 λ 2 σ y * y i λ 3 ( σ + γ ) y y i * ( σ + γ ) y y i ( c y ( σ + γ ) y i ) = ( λ 1 + λ 1 α 2 + λ 2 β 2 ) ( x x * ) 2 + ( λ 2 + λ 1 α 2 + λ 2 β 2 ) ( y y * ) 2 λ 2 σ y * y i ( σ + γ ) 2 y y i ( c y ( σ + γ ) y i ) 2 .
Let λ 1 = λ 2 . If (19) holds, then
D + V ( t ) | ( 3 ) < 0 .
Therefore, by using a Lyapunov–Lasalle type theorem [21], we have x x * , y y * , y i y i * as t . □

5. Existence of Hopf Bifurcation

In this section, we shall find the conditions under which a Hopf bifurcation may occur around the positive equilibrium point E * when delay τ passes through some critical values.
From (18) we can see that if 1 α β < 0 and q 2 2 p s l 2 > 0 , then f ( 0 ) = s 2 n 2 < 0 , f ( + ) = + and f ( 0 ) > 0 . Hence, Equation (15) has a unique positive root. Then, we have the following lemma.
Lemma 1.
Let 1 α β < 0 and q 2 2 p s l 2 > 0 hold. Then Equation (15) with τ = τ 0 has a unique pair of ω 0 , τ 0 , ω 0 τ 0 < 2 π such that Γ ( i ω 0 , τ 0 ) = 0 , where ω 0 and τ 0 are given by Equations (16).
Let us suppose τ = τ k when ω = ω 0 . From (16) we have
τ k = 1 2 ω 0 arcsin p l ω 0 3 l s ω 0 n ω 0 3 + q ω 0 n l 2 ω 0 2 + n 2 + 2 k π ω 0 , k = 0 , 1 , 2 , .
Theorem 5.
Assume 1 α β < 0 and q 2 2 p s l 2 > 0 . If τ [ 0 , τ 0 ) , then the equilibrium E * of system (3) is asymptotically stable; if τ > τ 0 , then E * is unstable. τ = τ k ( k = 0 , 1 , 2 , ) are Hopf bifurcation values for system (3), where τ k is defined by (21), that is, a family of periodic solutions bifurcates from E * as τ passes through the critical value τ k .
Proof. 
Let λ ( τ ) = a ( τ ) + i ω ( τ ) be a root of Equation (15). Substituting λ = a + i ω in (15) and separating real and imaginary parts, we get the transcendental equations
a 3 3 a ω 2 + p a 2 p ω 2 + q a + s + e 2 τ a ( l a + n ) cos ( 2 τ ω ) + e 2 τ a l ω sin ( 2 τ ω ) = 0
and
3 a 2 ω ω 3 + 2 p a ω + q ω e 2 τ a ( l a + n ) sin ( 2 τ ω ) + e 2 τ a l ω cos ( 2 τ ω ) = 0 .
Obviously, E * is stable when τ = 0 . Hence, E * remains stable for τ < τ 0 ( τ 0 is the smallest value for which where is a solution to (15) with real part zero). We now have to show that d a d τ τ = τ 0 > 0 when ω = ω 0 .
Now differentiating the two sides of Equation (22) with respect to τ , at a = 0 , τ = τ 0 , ω = ω 0 , we can obtain:
[ 3 ω 0 2 + q 2 n τ 0 cos ( 2 τ 0 ω 0 ) + l cos ( 2 τ 0 ω 0 ) 2 l τ 0 ω 0 sin ( 2 τ 0 ω 0 ) ] d a d τ + [ 2 p ω 0 2 n τ 0 sin ( 2 τ 0 ω 0 ) + 2 l τ 0 ω 0 cos ( 2 τ 0 ω 0 ) + l sin ( 2 τ 0 ω 0 ) ] d ω d τ + [ 2 n ω 0 sin ( 2 τ 0 ω 0 ) + 2 l ω 0 2 cos ( 2 τ 0 ω 0 ) ] = 0
and
[ 2 p ω 0 + 2 n τ 0 sin ( 2 τ 0 ω 0 ) 2 l τ 0 ω 0 cos ( 2 τ 0 ω 0 ) l sin ( 2 τ 0 ω 0 ) ] d a d τ + [ 3 ω 0 2 + q 2 n τ 0 cos 2 τ 0 ω 0 + l cos ( 2 τ 0 ω 0 ) 2 l τ 0 ω 0 sin ( 2 τ 0 ω 0 ) ] d ω d τ + [ 2 n ω 0 cos ( 2 τ 0 ω 0 ) + 2 l ω 0 2 sin ( 2 τ 0 ω 0 ) ] = 0 .
Hence,
[ ( 3 ω 0 2 + q 2 n τ 0 cos ( 2 τ 0 ω 0 ) 2 l τ 0 ω 0 sin ( 2 τ 0 ω 0 ) + l cos ( 2 τ 0 ω 0 ) ) 2 + ( 2 p ω 0 + 2 n τ 0 sin ( 2 τ 0 ω 0 ) l sin ( 2 τ 0 ω 0 ) 2 τ 0 l ω 0 cos ( 2 τ 0 ω 0 ) ) 2 ] d a d τ | a = 0 , ω = ω 0 , τ = τ 0 = 1 l 2 ω 0 2 + n 2 [ 6 n 2 ω 0 2 + 6 l 2 ω 0 8 + 2 n 2 ω 0 2 ( q 2 2 p s l 2 ) + 2 l 2 ω 0 4 ( q 2 2 p s l 2 ) + 4 n 2 ω 0 4 ( p 2 2 q ) + 4 l 2 ω 0 6 ( p 2 2 q ) ] .
Since q 2 2 p s l 2 > 0 , then d a d τ > 0 at a = 0 , τ = τ 0 , ω = ω 0 , the transversality condition holds. Hence, Hopf bifurcation occurs at ω = ω 0 , τ = τ 0 . □

6. Direction and Stability of the Hopf Bifurcation

In this section, we shall discuss the direction and stability of the Hopf bifurcation via the method introduced in [22].
Let u 1 ( t ) = x ( t ) x * , u 2 ( t ) = y ( t ) y * , u 3 ( t ) = y i ( t ) y i * , x i ( t ) = u i ( τ t ) , (i = 1, 2, 3), τ = τ k + μ . System (3) is transformed into functional differential equation (FDE) in C = C ( [ 1 , 0 ] , R 3 ) given by
d u d t = L μ ( u t ) + f ( μ , u t ) ,
where u ( t ) = ( u 1 ( t ) , u 2 ( t ) , u 3 ( t ) ) R 3 and
L μ ( ξ ) = ( τ k + μ ) x * 0 0 0 r + β x * 2 y * σ 0 c σ γ ξ 1 ( 0 ) ξ 2 ( 0 ) ξ 3 ( 0 )
+ ( τ k + μ ) 0 α x * 0 β y * 0 0 0 0 0 ξ 1 ( 1 ) ξ 2 ( 1 ) ξ 3 ( 1 ) ,
f ( μ , ξ ) = ( τ k + μ ) ξ 1 2 ( 0 ) α ξ 1 ( 0 ) ξ 2 ( 1 ) β ξ 1 ( 1 ) ξ 2 ( 0 ) ξ 2 2 ( 0 ) 0 .
By the Riesz representation theorem, there exists a 3 × 3 matrix function η ( θ , μ ) of bounded variation for θ [ 1 , 0 ] , such that
L μ ξ = 1 0 d η ( θ , 0 ) ξ ( θ )
for ξ C .
In practice, one can choose:
η ( θ , μ ) = ( τ k + μ ) x * 0 0 0 r + β x * 2 y * σ 0 c σ γ ξ 1 ( 0 ) ξ 2 ( 0 ) ξ 3 ( 0 ) δ ( θ )
( τ k + μ ) 0 α x * 0 β y * 0 0 0 0 0 δ ( θ + 1 ) ,
where δ ( θ ) denotes the Dirac delta function. For ξ C ( [ 1 , 0 ] , R 3 ) , we define
A ( μ ) ξ = d ξ ( θ ) d θ , θ [ 1 , 0 ) , 1 0 d η ( μ , s ) ξ ( s ) , θ = 0 ,
and
R ( μ ) ξ = 0 , θ [ 1 , 0 ) , f ( μ , ξ ) , θ = 0 .
Then, system (24) is equivalent to
u ˙ t = A ( μ ) u t + R ( μ ) u t ,
where u t ( θ ) = u ( t + θ ) , θ [ 1 , 0 ] .
For ς C 1 ( [ 1 , 0 ] , ( R 3 ) * ) , we define the adjoint operator A * of A as
A * ς ( s ) = d ς ( s ) d s , s ( 0 , 1 ] , 1 0 d η ( t , 0 ) ς ( t ) , s = 0 ,
and a bilinear inner product given by
< ς ( s ) , ξ ( θ ) > = ς ¯ ( 0 ) ξ ( 0 ) 1 0 ζ θ θ ς ¯ ( ζ θ ) d η ( θ ) ξ ( ζ ) d ζ ,
where η ( θ ) = η ( θ , 0 ) . Clearly, A ( 0 ) and A * are a pair of adjoint operators. From the discussions in Section 5, we know that ± i ω 0 τ k are eigenvalues of A ( 0 ) . Thus, they are also eigenvalues of A * . In the next, we shall calculate the eigenvector q ( θ ) of A ( 0 ) and eigenvector q * ( s ) of A * corresponding to + i ω 0 τ k and i ω 0 τ k , respectively.
Let q ( θ ) = ( 1 , υ , ν ) e i ω 0 τ k θ be the eigenvector of A ( 0 ) corresponding to + i ω 0 τ k , then A ( 0 ) q ( θ ) = i ω 0 τ k q ( θ ) . From the definition of A ( 0 ) and (25), (27) and (28), we can obtain that
τ k i ω 0 + x * α x * e i ω 0 τ k 0 β y * e i ω 0 τ k i ω 0 + r β x * + 2 y * σ 0 c i ω 0 + σ + γ . q ( 0 ) = 0 0 0 .
From above, it is easy to obtain q ( 0 ) = ( 1 , υ , ν ) , where υ = α x * e i ω 0 τ k i ω 0 + x * , ν = c υ i ω 0 + σ + γ . Similarly, assuming that q * ( s ) = D ( 1 , υ * , ν * ) e i ω 0 τ k s is the eigenvector of A * corresponding to i ω 0 τ k , we can get υ * = β y * e i ω 0 τ k i ω 0 + x * , ν * = σ υ * i ω 0 + σ + γ by the definition of A * and (25)–(27). In order to ensure < q * ( s ) , q ( θ ) > = 1 , we need to determine the value of D. From (31), we have
< q * ( s ) , q ( θ ) > = D ¯ ( 1 , υ * ¯ , ν * ¯ ) ( 1 , υ , ν ) 1 0 ζ = 0 θ D ¯ ( 1 , υ * ¯ , ν * ¯ ) e i ω 0 τ k ( ζ θ ) d η ( θ ) ( 1 , υ , ν ) e i ω 0 τ k ζ d ζ = D ¯ { 1 + υ υ * ¯ + ν ν * ¯ 1 0 ( 1 , υ * ¯ , ν * ¯ ) θ e i ω 0 τ k θ d η ( θ ) ( 1 , υ , ν ) } = D ¯ { 1 + υ υ * ¯ + ν ν * ¯ + ( α x * υ β y * υ * ¯ ) τ k e i ω 0 τ k } .
Thus, we can select D as
D = 1 1 + υ υ * ¯ + ν ν * ¯ + ( α x * υ β y * υ * ¯ ) τ k e i ω 0 τ k .
Using the same notations as in [22], we first calculate the coordinates to describe the center manifold C 0 at μ = 0 . Let u t be the solution of (30) when μ = 0 . Define
z ( t ) = < q * , u t > , W ( w , θ ) = u t ( θ ) 2 R e { z ( t ) q ( θ ) } .
On the center manifold C 0 we have
W ( t , θ ) = W ( z ( t ) , z ¯ ( t ) , θ ) ,
where
W ( z , z ¯ , θ ) = W 20 ( θ ) z 2 2 + W 11 ( θ ) z z ¯ + W 02 ( θ ) z ¯ 2 2 + W 30 ( θ ) z 3 6 + ,
z and z ¯ are local coordinates for center manifold C 0 in the direction of q * and q * ¯ . Note that W is real if x t is real. We only consider real solutions. For solution x t C 0 of (30), since μ = 0 , we have
z ˙ ( t ) = i ω 0 τ k z + q * ¯ ( 0 ) f ( 0 , W ( z , z ¯ , θ ) + 2 R e { z q ( θ ) } ) = i ω 0 τ k z + q * ¯ ( 0 ) f 0 ( z , z ¯ ) .
We rewrite this equation as
z ˙ ( t ) = i ω 0 τ k z ( t ) + g ( z , z ¯ ) ,
where
g ( z , z ¯ ) = q * ¯ ( 0 ) f 0 ( z , z ¯ ) = g 20 z 2 2 + g 11 z z ¯ + g 02 z ¯ 2 2 + g 21 z 2 z ¯ 2 + .
It follows from (32) and (33) that:
x t ( θ ) = W ( t , θ ) 2 R e { z ( t ) q ( t ) } = W 20 ( θ ) z 2 2 + W 11 ( θ ) z z ¯ + W 02 ( θ ) z ¯ 2 2 + ( 1 , α , β ) e i ω 0 τ k θ z + ( 1 , α ¯ , β ¯ ) e i ω 0 τ k θ z ¯ + .
It follows together with (26) that:
g ( z , z ¯ ) = q * ¯ ( 0 ) f 0 ( z , z ¯ ) = q * ¯ ( 0 ) f ( 0 , x t ) = τ k D ¯ ( 1 , υ * ¯ , ν * ¯ ) x 1 t 2 ( 0 ) α x 1 t ( 0 ) x 2 t ( 1 ) β x 1 t ( 1 ) x 2 t ( 0 ) x 2 t 2 ( 0 ) 0 = τ k D ¯ { [ z + z ¯ + W 20 ( 1 ) ( 0 ) z 2 2 + W 11 ( 1 ) ( 0 ) z z ¯ + W 02 ( 1 ) ( 0 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] 2 + α [ z + z ¯ + W 20 ( 1 ) ( 0 ) z 2 2 + W 11 ( 1 ) ( 0 ) z z ¯ + W 02 ( 1 ) ( 0 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] × [ υ e i ω 0 τ k z + υ ¯ e i ω 0 τ k z ¯ + W 20 ( 2 ) ( 1 ) z 2 2 + W 11 ( 2 ) ( 1 ) z z ¯ + W 02 ( 2 ) ( 1 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] υ * ¯ β [ e i ω 0 τ k z + e i ω 0 τ k z ¯ + W 20 ( 1 ) ( 1 ) z 2 2 + W 11 ( 1 ) ( 1 ) z z ¯ + W 02 ( 1 ) ( 1 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] × [ υ z + υ ¯ z ¯ + W 20 ( 2 ) ( 0 ) z 2 2 + W 11 ( 2 ) ( 0 ) z z ¯ + W 02 ( 2 ) ( 0 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] + υ * ¯ [ υ z + υ ¯ z ¯ + W 20 ( 2 ) ( 0 ) z 2 2 + W 11 ( 2 ) ( 0 ) z z ¯ + W 02 ( 2 ) ( 0 ) z ¯ 2 2 + o ( | ( z , z ¯ ) | 3 ) ] 2 } .
Comparing the coefficients with (34), we have:
g 20 = 2 τ k D ¯ ( 1 υ * ¯ υ 2 + υ α e i ω 0 τ k υ * ¯ υ β e i ω 0 τ k ) , g 11 = 2 τ k D ¯ ( 1 υ * ¯ υ υ ¯ + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 υ * ¯ ( υ e i ω 0 τ k + υ e i ω 0 τ k ) 2 ) , g 02 = 2 τ k D ¯ ( 1 υ * ¯ υ ¯ 2 + υ α e i ω 0 τ k υ * ¯ υ ¯ β e i ω 0 τ k ) , g 21 = 2 τ k D ¯ { 2 W 11 ( 1 ) ( 0 ) + W 20 ( 1 ) ( 0 ) + α [ W 11 ( 2 ) ( 1 ) + W 20 ( 2 ) ( 0 ) 2 + e i ω 0 τ k W 20 ( 1 ) ( 0 ) 2 + υ W 11 ( 1 ) ( 0 ) e i ω 0 τ k ] υ * ¯ β [ W 11 ( 2 ) ( 0 ) e i ω 0 τ k + W 20 ( 2 ) ( 0 ) e i ω 0 τ k 2 + υ * ¯ W 11 ( 1 ) ( 1 ) + υ ¯ W 02 ( 1 ) ( 1 ) 2 ] υ W 11 ( 2 ) ( 0 ) 2 υ * ¯ υ W 20 ( 2 ) ( 0 ) } .
Since there are W 20 ( θ ) and W 11 ( θ ) in g 21 , we still need to compute them. From (30) and (32), we have
W ˙ = u ˙ t z ¯ ˙ q z ¯ ˙ q ¯ = A W 2 R e { q * ¯ ( 0 ) f 0 q ( θ ) } , θ [ 1 , 0 ) , A W 2 R e { q * ¯ ( 0 ) f 0 q ( θ ) } + f 0 , θ = 0 , = A W + H ( z , z ¯ , θ ) ,
where
H ( z , z ¯ , θ ) = H 20 ( θ ) z 2 2 + H 11 ( θ ) z z ¯ + H 02 ( θ ) z ¯ 2 2 + .
Substituting the corresponding series into (38) and comparing the coefficients, we obtain
( A 2 i ω 0 τ 0 ) W 20 ( θ ) = H 20 , A W 11 ( θ ) = H 11 , .
From (38), we know that for θ [ 1 , 0 ) ,
H ( z , z ¯ , θ ) = q * ¯ ( 0 ) f 0 q ( θ ) q * ( 0 ) f ¯ 0 q ¯ ( θ ) = g ( z , z ¯ ) q ( θ ) g ( z , z ¯ ) q ¯ ( θ ) .
Comparing the coefficients with (39) gives that
H 20 ( θ ) = g 20 q ( θ ) g ¯ 02 q ¯ ( θ ) ,
and
H 11 ( θ ) = g 11 q ( θ ) g ¯ 11 q ¯ ( θ ) .
From (40), (42) and the definition of A, it follows that
W ˙ 20 = 2 i ω 0 τ k W 20 ( θ ) + g 20 q ( θ ) + g ¯ 02 q ¯ ( θ ) .
Notice that q ( θ ) = ( 1 , υ , ν ) e i ω 0 τ k θ , hence
W 20 ( θ ) = i g 20 ω 0 τ k q ( 0 ) e i ω 0 τ k θ + i g ¯ 02 3 ω 0 τ k q ¯ ( 0 ) e i ω 0 τ k θ + E 1 e 2 i ω 0 τ k θ ,
where E 1 = ( E 1 ( 1 ) , E 1 ( 2 ) , E 1 ( 3 ) ) R 3 is a constant vector. Similarly, from (40) and (43), we obtain
W 11 ( θ ) = i g 11 ω 0 τ k q ( 0 ) e i ω 0 τ k θ + i g ¯ 11 ω 0 τ k q ¯ ( 0 ) e i ω 0 τ k θ + E 2 ,
where E 2 = ( E 2 ( 1 ) , E 2 ( 2 ) , E 2 ( 3 ) ) R 3 is also a constant vector.
In what follows, we shall seek appropriate E 1 and E 2 . From the definition of A and (40), we obtain
1 0 d η ( θ ) W 20 ( θ ) = 2 i ω 0 τ k W 20 ( θ ) H 20 ( θ ) ,
and
1 0 d η ( θ ) W 11 ( θ ) = H 11 ( θ ) ,
where η ( θ ) = η ( 0 , θ ) . By (38), we have
H 20 ( θ ) = g 20 q ( 0 ) g ¯ 20 q ¯ ( 0 ) + 2 τ k 1 + υ α e i ω 0 τ k υ 2 υ β e i ω 0 τ k 0 ,
and
H 11 ( θ ) = g 11 q ( 0 ) g ¯ 11 q ¯ ( 0 ) + 2 τ k 1 + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 υ υ ¯ υ e i ω 0 τ k + υ e i ω 0 τ k 2 0 .
Substituting (44) and (48) into: (46), we obtain
( 2 i ω 0 τ k I 1 0 e 2 i ω 0 τ k θ d η ( θ ) ) E 1 = 2 τ k 1 + υ α e i ω 0 τ k υ 2 υ β e i ω 0 τ k 0 ,
which leads to:
2 i ω 0 + x * α x * e 2 i ω 0 τ k 0 β y * e 2 i ω 0 τ k 2 i ω 0 + r β x * + 2 y * σ 0 c 2 i ω 0 + σ + γ E 1 = 2 1 + υ α e i ω 0 τ k υ 2 υ β e i ω 0 τ k 0 .
It follows that:
E 1 ( 1 ) = 2 Δ 1 1 + υ α e i ω 0 τ k * α x * e 2 i ω 0 τ k 0 υ 2 υ β e i ω 0 τ k 2 i ω 0 + r β x * + 2 y * σ 0 c 2 i ω 0 + σ + γ , E 1 ( 2 ) = 2 Δ 1 2 i ω 0 + x * 1 + υ α e i ω 0 τ k 0 β y * e 2 i ω 0 τ k υ 2 υ β e i ω 0 τ k σ 0 0 2 i ω 0 + σ + γ , E 1 ( 3 ) = 2 Δ 1 2 i ω 0 + x * α x * e 2 i ω 0 τ k 1 + υ α e i ω 0 τ k β y * e 2 i ω 0 τ k 2 i ω 0 + r β x * + 2 y * υ 2 υ β e i ω 0 τ k 0 c 0 ,
where
Δ 1 = 2 i ω 0 + x * α x * e 2 i ω 0 τ k 0 β y * e 2 i ω 0 τ k 2 i ω 0 + r β x * + 2 y * σ 0 c 2 i ω 0 + σ + γ .
Similarly, substituting (45) and (49) into (47), we can obtain:
x * α x * 0 β y * r β x * + 2 y * σ 0 c σ + γ E 2 = 1 + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 υ υ ¯ υ e i ω 0 τ k + υ e i ω 0 τ k 2 0 ,
and hence
E 2 ( 1 ) = 2 Δ 2 1 + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 α x * 0 υ υ ¯ υ e i ω 0 τ k + υ e i ω 0 τ k 2 r β x * + 2 y * σ 0 c σ + γ , E 2 ( 2 ) = 2 Δ 2 x * 1 + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 0 β y * υ υ ¯ υ e i ω 0 τ k + υ e i ω 0 τ k 2 σ 0 0 σ + γ , E 2 ( 3 ) = 2 Δ 2 x * α x * 1 + υ ¯ α e i ω 0 τ k + υ α e i ω 0 τ k 2 β y * r β x * + 2 y * υ υ ¯ υ e i ω 0 τ k + υ e i ω 0 τ k 2 0 c 0 ,
where
Δ 2 = x * α x * 0 β y * r β x * + 2 y * σ 0 c σ + γ .
Thus, we can determine W 20 ( θ ) and W 11 ( θ ) from (44) and (45). Furthermore, g 21 in (37) can be expressed by the parameters and delay. Then, we can compute the following values:
c 1 ( 0 ) = i 2 ω 0 τ k ( g 20 g 11 2 | g 11 | 2 | g 02 | 2 3 ) + g 21 2 , μ 2 = R e { c 1 ( 0 ) } R e { λ ( τ k ) } , β 2 = 2 R e { c 1 ( 0 ) } , T 2 = I m { c 1 ( 0 ) } + μ 2 I m { λ ( τ k ) } ω 0 τ k .
In conclusion, we have the results as following.
Theorem 6.
(1) If μ 2 < 0 ( μ 2 > 0 ), Hopf bifurcation is subcritical (supercritical); (2) If β 2 < 0 ( β 2 > 0 ), the bifurcating periodic solution is stable (unstable); (3) If T 2 > 0 ( T 2 < 0 ), the period increases (decreases).

7. Numerical Simulations

In this section, we provide numerical examples of Theorems 5 and 6 by using Matlab.
Let k = 3 , α = 1 , r = 1 , β = 1.5 , c = 2 , σ = 0.5 , γ = 0.5 . Hence, p = 5 , q = 6.36 , s = 2.16 , l = 3.24 , α + β > 2 , q 2 2 p s l 2 = 8.352 > 0 , τ 0 = 1.3434655328 , E * ( x * , y * , y i * ) = ( 1.2 , 1.8 , 3.6 ) .
(i) τ = 1.2 , τ < τ 0 . In this case, the numerical simulation (see Figure 1) shows that the predator and prey populations spiral toward the equilibrium E * ( 1.2 , 1.8 , 3.6 ) ;
(ii) By calculation, we obtain τ 0 = 1.3434655328 . By Theorem 5.1, a Hopf bifurcation occurs when τ = τ 0 . Select τ = 1.8 . From Figure 2, we can find that both the predator and prey populations reach periodic oscillations around the equilibrium E * ( x * , y * , y i * ) in finite time;
(iii) From our calculations, we can see that ω 0 = 0.6514262661 and λ ( τ 0 ) = 0.6127495603 + 0.4261524984 i . From the formulae (50), we can obtain c 1 ( 0 ) = 10.29371206 73.4416490 i , μ 2 > 0 ,   β 2 < 0 and T 2 > 0 . Hence, the Hopf bifurcation is supercritical, the direction of the bifurcation is τ > τ 0 and these bifurcating periodic solutions are stable (see Figure 2).

8. Discussion

In this paper, a new dynamics for a predator–prey model with staged structure and time delay has been analyzed. We discuss the influence of the parameter τ on the dynamics of the system (3). The system (3) is permanent under some conditions. The local and global stability of the positive equilibrium is presented. By choosing τ as a bifurcation parameter, we prove that the delay loss of stability phenomenon appears under the conditions 1 α β < 0 and q 2 2 p s l 2 > 0 . That is to say, there is a critical value τ 0 of τ such that system (3) is stable in the range τ ( 0 , τ 0 ) at positive equilibrium E * (see Figure 1); when τ = τ 0 , a Hopf bifurcation occurs around E * ( x * , y * , y i * ) ; when τ > τ 0 , the system is unstable (see Figure 2) and there are always Hopf bifurcations near the positive equilibrium E * ( x * , y * , y i * ) when τ takes other critical values. We derive explicit formulae for determining the properties of Hopf bifurcation at the critical value of τ 0 via the ideas of Hassard et al. If we do not consider the staged structure of the predator, system (2) has no periodic solutions, which shows that the staged structure of the predator can severely affect the dynamical behavior of the system. However, for system (3), there is no chaotic behavior of the system (3) by numerical simulations.

Funding

This work is supported by the Natural Science Foundation of Henan (222300420521), Program for Innovative Research Team (in Science and Technology) in University of Henan Province (21IRTSTHN014) and Nanhu Scholars Program for Young Scholars of XYNU.

Data Availability Statement

Not applicable.

Acknowledgments

We would like to thank the anonymous referees for their careful reading of the original manuscript and their many valuable comments and suggestions that greatly improve the presentation of this work.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Conejero, J.A.; Murillo-Arcila, M.; Seoane, J.M.; Seoane-Sepúlveda, J.B. When Does Chaos Appear While Driving? Learning Dynamical Systems via Car-Following Models. Math. Mag. 2022, 95, 302–313. [Google Scholar] [CrossRef]
  2. Zhou, X.; Shi, X.; Wei, M. Dynamical behavior and optimal control of a stochastic mathe-matical model for cholera. Chaos Solitons Fractals 2022, 156, 111854. [Google Scholar] [CrossRef]
  3. Owolabi, K.M. Computational dynamics of predator–prey model with the power-law k-ernel. Results Phys. 2021, 21, 103810. [Google Scholar] [CrossRef]
  4. Cao, Y.; Alamri, S.; Rajhi, A.A.; Anqi, A.E.; Riaz, M.B.; Elagan, S.K.; Jawa, T.M. A novel piece-wise approach to modeling interactions in a food web model Author links open overlay panel. Results Phys. 2021, 31, 104951. [Google Scholar] [CrossRef]
  5. Das, B.K.; Sahoo, D.; Samanta, G.P. Impact of fear in a delay-induced predator–prey system with intraspecific competition within predator species. Math. Comput. Simul. 2022, 191, 134–156. [Google Scholar] [CrossRef]
  6. Liu, Q.; Jiang, D.; Hayat, T. Dynamics of stochastic predator-prey models with distributed delay and stage structure for prey. Int. J. Biomath. 2021, 14, 2150020. [Google Scholar] [CrossRef]
  7. Alsakaji, H.J.; Kundu, S.; Rihan, F.A. Delay differential model of one-predator two-prey system with Monod-Haldane and holling type II functional responses. Appl. Math. Comput. 2021, 397, 125919. [Google Scholar] [CrossRef]
  8. Yousef, F.; Semmar, B.; Al Nasr, K. Dynamics and simulations of discretized Caputo-conformable fractional-order Lotka-Volterra models. Nonlinear Eng. 2022, 11, 100–111. [Google Scholar] [CrossRef]
  9. Yousef, F.; Semmar, B.; Al Nasr, K. Incommensurate conformable-type three-dimensional Lotka–Volterra model: Discretization, stability, and bifurcation. Arab. J. Basic Appl. Sci. 2022, 29, 113–120. [Google Scholar] [CrossRef]
  10. Chen, L.S.; Song, X.Y.; Lu, Z.Y. Mathematica Models and Methods in Ecology; Sichuan Science and Technology: Chengdu, China, 2003. (In Chinese) [Google Scholar]
  11. He, X.Z. The Lyapunov functionals for delay Lotka Volterra type models. SIAM J. Appl. Math. 1998, 58, 1222–1236. [Google Scholar] [CrossRef]
  12. Gopalsamy, K.; He, X.Z. Global stability in n-species competition modelled by “pure-delay type” systems II: Nonautonomous case. Can. Appl. Math. Q. 1998, 6, 17–43. [Google Scholar]
  13. He, X.Z. Global stability in nonautonomous Lotka–Volterra systems of “pure-delay type”. Differ. Integral Equations 1998, 11, 293–310. [Google Scholar]
  14. Wang, W.D.; Ma, Z.E. Harmless delays for uniform peristence. J. Math. Anal. Appl. 1991, 158, 256–268. [Google Scholar]
  15. Cui, J.A.; Chen, L.S.; Wang, W.D. The effect of dispersal on population growth with stage- structure. Comput. Appl. Math. 2000, 39, 91–102. [Google Scholar] [CrossRef] [Green Version]
  16. Song, X.Y.; Cui, J.A. The Stage-structured predator–prey System with Delay and Harvesting. Appl. Anal. 2002, 81, 1127–1142. [Google Scholar] [CrossRef]
  17. Aiello, W.G.; Freedman, H. A time-delay model of single-species growth with stage structure. Math. Biosci. 1990, 101, 139–153. [Google Scholar] [CrossRef]
  18. Cao, Y.; Fan, J.; Gard, T.C. The effects of state-dependent time delay on a stage-structured population growth model. Nonlinear Anal. Theory Methods Appl. 1992, 19, 95–105. [Google Scholar] [CrossRef]
  19. Walter, W. Ordinary Differential Equations. Graduate Texts in Mathematics; Springer: Now York, NY, USA, 1998; Volume 18. [Google Scholar]
  20. Hale, J.K.; Waltman, P. Persistence in inffnite-dimensional systems. SIAM J. Math. Anal. 1989, 20, 388–396. [Google Scholar] [CrossRef]
  21. Kuang, Y. Delay Differential Equations with Applications in Population Dynamics; Academic Press: London, UK, 2004. [Google Scholar]
  22. Hassard, B.D.; Kazarini, N.D.; Wan, Y.H. Theory and Application of Hopf Bifurcation; London Mathematical Society Lecture Note Series, 41; Cambridge University Press: Cambridge, UK, 1981. [Google Scholar]
Figure 1. k = 3 , α = 1 , r = 1 , β = 1.5 , c = 2 , σ = 0.5 , γ = 0.5 , τ = 1.2 . The solution tends to the positive equilibrium. The initial value is constant function x ( t ) 1 , y ( t ) 0.1 , y i ( t ) 0.5 for t [ τ , 0 ] .
Figure 1. k = 3 , α = 1 , r = 1 , β = 1.5 , c = 2 , σ = 0.5 , γ = 0.5 , τ = 1.2 . The solution tends to the positive equilibrium. The initial value is constant function x ( t ) 1 , y ( t ) 0.1 , y i ( t ) 0.5 for t [ τ , 0 ] .
Axioms 11 00575 g001
Figure 2. Trajectory of the system with same parameters as Figure 1. except that τ = 1.8 . The solution tends to the periodic solution. The initial value is constant function x ( t ) 1 , y ( t ) 0.1 , y i ( t ) 0.5 for t [ τ , 0 ] .
Figure 2. Trajectory of the system with same parameters as Figure 1. except that τ = 1.8 . The solution tends to the periodic solution. The initial value is constant function x ( t ) 1 , y ( t ) 0.1 , y i ( t ) 0.5 for t [ τ , 0 ] .
Axioms 11 00575 g002
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhou, X. Stability and Hopf Bifurcation Analysis of a Stage-Structured Predator–Prey Model with Delay. Axioms 2022, 11, 575. https://doi.org/10.3390/axioms11100575

AMA Style

Zhou X. Stability and Hopf Bifurcation Analysis of a Stage-Structured Predator–Prey Model with Delay. Axioms. 2022; 11(10):575. https://doi.org/10.3390/axioms11100575

Chicago/Turabian Style

Zhou, Xueyong. 2022. "Stability and Hopf Bifurcation Analysis of a Stage-Structured Predator–Prey Model with Delay" Axioms 11, no. 10: 575. https://doi.org/10.3390/axioms11100575

APA Style

Zhou, X. (2022). Stability and Hopf Bifurcation Analysis of a Stage-Structured Predator–Prey Model with Delay. Axioms, 11(10), 575. https://doi.org/10.3390/axioms11100575

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