- Email: [email protected]
Global dynamics of a predator–prey model with time delay and stage structure for the prey
Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Global dynamics of a predator–prey model with time delay and stage structure for the prey✩ Rui Xu ∗ Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China
article
info
Article history: Received 11 August 2010 Accepted 29 December 2010 Keywords: Stage structure Time delay Hopf bifurcation LaSalle invariant principle Global stability
abstract A stage-structured predator–prey system with Holling type-II functional response and time delay due to the gestation of predator is investigated. By analyzing the characteristic equations, the local stability of each of feasible equilibria of the system is discussed and the existence of a Hopf bifurcation at the coexistence equilibrium is established. By means of the persistence theory on infinite dimensional systems, it is proven that the system is permanent if the coexistence equilibrium exists. By using Lyapunov functionals and LaSalle invariant principle, it is shown that the trivial equilibrium is globally stable when both the predator-extinction equilibrium and the coexistence equilibrium are not feasible, and that the predator-extinction equilibrium is globally asymptotically stable if the coexistence equilibrium does not exist, and sufficient conditions are derived for the global stability of the coexistence equilibrium. Numerical simulations are carried out to illustrate the main theoretical results. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction In [1], based on experiments, Holling suggested three different kinds of functional responses for different kinds of species to model the phenomena of predation, which made the standard Lotka–Volterra system more realistic. Letting x(t ) and y(t ) represent the densities of the prey and the predator at time t, respectively, a standard predator–prey model with Holling type functional response is of the form x˙ (t ) = xg (x) − φ(x)y, y˙ (t ) = eφ(x)y − c (y).
(1.1)
In (1.1), the functions g (x) and c (y) represent the specific growth rate of the prey in the absence of predation and the mortality rate of the predator, respectively; the function φ(x) is called ‘‘functional response’’ representing the prey consumption per unit time. The most popular functional response used in the modelling of predator–prey systems is Holling type-II with φ(x) = fx/(1 + mx) which takes into account the time a predator uses in handling the prey being captured. There has been a large body of work in the literature on the global dynamics of predator–prey systems with Holling type-II functional responses (see, for example, [2–6] and references cited therein). It is well-known that time delays of one type or another have been incorporated into mathematical models of population dynamics by many researchers. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the
✩ This work was supported by the National Natural Science Foundation of China (Nos. 11071254, 10671209).
∗ Corresponding address: Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, No. 97 Heping West Road, Shijiazhuang 050003, Hebei Province, PR China. Tel.: +86 311 87994035; fax: +86 311 87994007. E-mail address: [email protected]. 1468-1218/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2010.12.029
2152
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
population to fluctuate. We refer to the monographs of Cushing [7], Gopalsamy [8] and Kuang [9] for general delayed biological systems and to Bartlett [10], Beretta and Kuang [11] and Wangersky and Cunningham [12] and references cited therein for studies on delayed predator–prey systems. Time delay due to gestation is a common example, because generally the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Therefore, more realistic models of population interactions should take into account the effect of time delays. In [12], Wangersky and Cunningham proposed and studied the following non-Kolmogorov-type predator–prey model: x˙ (t ) = x(t )(r1 − ax(t ) − a1 y(t )),
(1.2)
y˙ (t ) = a2 x(t − τ )y(t − τ ) − r2 y(t ).
This model assumes that a duration of τ time units elapses when an individual prey is killed and the moment when the corresponding addition is made to the predator population. Stage-structure is a natural phenomenon and represents, for example, the division of a population into immature and mature individuals. As is common, the dynamics-eating habits, susceptibility to predators, etc. are often quite different in these two sub-populations. Hence, it is of ecological importance to investigate the effects of such a subdivision on the interaction of species. There has been a significant body of work on stage-structured population models in last two decades (see, for example, [13–19]). In [18], Xu et al. considered a ratio-dependent predator–prey model with stage structure for the prey. The prey individuals were classified as belonging either the immature or the mature and it was assumed that the immature prey is not at risk of being attacked by the predator. This seems reasonable for a number of mammals, where the immature prey concealed in the mountain cave, are raised by their parents; they do not necessarily go out for seeking food, the rate they are attacked by the predators can be ignored. By constructing appropriate Lyapunov functions, sufficient conditions are obtained for the global asymptotic stability of nonnegative equilibria of the model. Motivated by the work of Chen and Jing [2], Wangersky and Cunningham [12] and Xu et al. [18], in the present paper, we are concerned with the combined effects of stage structure for the prey and time delay due to the gestation of the predator on the global dynamics of a predator–prey model with Holling type-II functional response. To this end, we consider the following delay differential system x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) = bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − y˙ (t ) =
a2 x2 (t − τ )y(t − τ ) 1 + mx2 (t − τ )
a1 x2 (t )y(t ) 1 + mx2 (t )
,
(1.3)
− ry(t ).
In system (1.3), x1 (t ) and x2 (t ) represent the densities of the immature and the mature prey at time t, respectively; y(t ) represents the density of the predator at time t. The model is derived under the following assumptions. (A1) The prey population: the birth rate is proportional to the existing mature population with a proportionality a > 0; the death rate of the immature population and the transformation rate from immature individuals to mature individuals are proportional to the existing immature population with proportionality constants r1 > 0 and b > 0, respectively; b1 > 0 is the intra-specific competition rate of the mature population; the death rate of the mature population is proportional to the existing mature population with a proportionality r2 > 0. (A2) The predator population: the growth of the species obeys a Holling type-II functional response. The predators feed only on the mature prey. r > 0 is the death rate of the predator, a1 > 0 is the capturing rate of the predator, a2 /a1 > 0 is the conversion rate of nutrients into the reproduction of the predator; m > 0 is the half saturation rate of the predator; τ ≥ 0 is a constant delay due to the gestation of the predator, that is, mature adult predators can only contribute to the reproduction of the predators’ biomass. The delay τ is based on the assumption that the change rate of predators depends on the number of the mature prey and of the predators presented at some previous time. The initial conditions for system (1.3) take the form x1 (θ ) = φ1 (θ ), x2 (θ ) = φ2 (θ ), y(θ ) = ψ(θ ), φ1 (θ ) ≥ 0, φ2 (θ ) ≥ 0, ψ(θ ) ≥ 0, θ ∈ [−τ , 0), φ1 (0) > 0, φ2 (0) > 0, ψ(0) > 0,
(1.4)
where (φ1 (θ ), φ2 (θ ), ψ(θ )) ∈ C ([−τ , 0], R3+0 ), the Banach space of continuous functions mapping the interval [−τ , 0] into R3+0 , where R3+0 = {(x1 , x2 , x3 ) : xi ≥ 0, i = 1, 2, 3}. It is well-known by the fundamental theory of functional differential equations [20], that system (1.3) has a unique solution (x1 (t ), x2 (t ), y(t )) satisfying initial conditions (1.4). It is easy to show that all solutions of system (1.3) corresponding to initial conditions (1.4) are defined on [0, +∞) and remain positive for all t ≥ 0. The organization of this paper is as follows. In the next section, by analyzing the corresponding characteristic equations, the local stability of each of the feasible equilibria of system (1.3) is discussed and the existence of a Hopf bifurcation at the coexistence equilibrium is established. In Section 3, by using the persistence theory on infinite dimensional systems
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2153
developed in [21], we prove that system (1.3) is permanent if the coexistence equilibrium exists. In Section 4, by using novel Lyapunov functionals and the LaSalle invariant principle, we show that both the prey and the predator go to extinction if both the predator-extinction equilibrium and the coexistence equilibrium are not feasible, and that the predator species goes to extinction when the coexistence equilibrium does not exist, and sufficient conditions are obtained for the global stability of the coexistence equilibrium. Numerical simulations are carried out in Section 5 to illustrate the main theoretical results. A brief discussion is given in Section 6 to conclude this work. 2. Local stability and Hopf bifurcation In this section, we study the local stability of each of feasible equilibria of system (1.3) and the existence of a Hopf bifurcation at the coexistence equilibrium. System (1.3) always has a trivial equilibrium E0 (0, 0, 0). If the following holds: (H1) ab > r2 (b + r1 ), then (1.3) has a predator-extinction equilibrium E1 (x01 , x02 , 0), where x01 =
a[ab − r2 (b + r1 )] b1 (b + r1 )2
,
x02 =
ab − r2 (b + r1 ) b1 (b + r1 )
.
(2.1)
Further, if (H1) and the following hold: (H2) (a2 − mr )[ab − r2 (b + r1 )] > b1 r (b + r1 ), then system (1.3) has a unique coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ), where ar
,
r
, (b + r1 )(a2 − mr ) a2 − mr a2 {(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 )} y∗ = . a1 (b + r1 )(a2 − mr )2
x∗1 =
x∗2 =
(2.2)
The characteristic equation of system (1.3) at the trivial equilibrium E0 (0, 0, 0) is of the form
(λ + r )[λ2 + (b + r1 + r2 )λ + r2 (b + r1 ) − ab] = 0.
(2.3)
It is readily seen from Eq. (2.3) that if ab < r2 (b + r1 ), then E0 is locally asymptotically stable; if ab > r2 (b + r1 ), then E0 is unstable. The characteristic equation of system (1.3) at the equilibrium E1 (x01 , x02 , 0) takes the form
[λ + (b + r1 + r2 + 2
2b1 x02
[
)λ + ab − r2 (b + r1 )] λ + r −
a2 x02 1 + mx02
e
−λτ
]
= 0.
(2.4)
Clearly, if (H1) holds, then the equation
λ2 + (b + r1 + r2 + 2b1 x02 )λ + ab − r2 (b + r1 ) = 0
(2.5)
has two roots with negative real parts. Other roots of Eq. (2.4) are determined by the equation
λ+r −
a2 x02 1 + mx02
e−λτ = 0.
(2.6)
Denote f (λ) = λ + r −
a2 x02 1 + mx02
e−λτ .
If (H1)–(H2) hold, it is easy to show that, for λ real, f (0) = −
(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) < 0, b1 (b + r1 )(1 + mx02 )
lim f (λ) = +∞.
λ→+∞
Hence, f (λ) = 0 has at least one positive real root. Therefore, if (H1)–(H2) hold, then the equilibrium E1 is unstable. If (H1) and the following hold: (H3) (a2 − mr )[ab − r2 (b + r1 )] < b1 r (b + r1 ),
2154
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
then it is easy to show that r − a2 x02 /(1 + mx02 ) > 0. By Theorem 3.2.1 in Kuang [9], we know that the equilibrium E1 is locally asymptotically stable for all τ ≥ 0. The characteristic equation of system (1.3) at the coexistence equilibrium E ∗ is of the form
λ3 + p2 λ2 + p1 λ + p0 + (q2 λ2 + q1 λ + q0 )e−λτ = 0,
(2.7)
where
[
p0 = r (b + r1 ) r2 + 2b1 x2 + ∗
p1 = r
∗
b + r1 + r2 + 2b1 x2 +
a1 y∗
(1 + mx∗2 )2 ∗
− ab ,
+ (b + r1 ) r2 + 2b1 x2 +
(1 + mx∗2 )2
(1 + mx∗2 )2 ∗ q0 = −r [(b + r1 )(r2 + 2b1 x2 ) − ab], q1 = −r (b + r1 + r2 + 2b1 x∗2 ), q2 = −r .
]
a1 y
a1 y ∗
p2 = b + r1 + r2 + 2b1 x∗2 +
a1 y∗
∗
(1 + mx∗2 )2
− ab, (2.8)
+ r,
When τ = 0, Eq. (2.7) becomes
λ3 + (p2 + q2 )λ2 + (p1 + q1 )λ + p0 + q0 = 0.
(2.9)
Clearly, p2 + q2 > 0. By calculation we derive that p0 + q0 =
a1 r (b + r1 )y∗
(1 + mx∗2 )2
> 0.
(2.10)
Hence, by the Routh–Hurwitz criterion, we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 , then the equilibrium E ∗ is locally asymptotically stable when τ = 0. If iω (ω > 0) is a solution of Eq. (2.7), separating real and imaginary parts, we have
ω3 − p1 ω = (q2 ω2 − q0 ) sin ωτ + q1 ω cos ωτ ,
(2.11)
−p2 ω2 + p0 = (q2 ω2 − q0 ) cos ωτ − q1 ω sin ωτ . Squaring and adding the two equations of (2.11), it follows that
ω6 + (p22 − 2p1 − q22 )ω4 + (p21 − 2p0 p2 + 2q0 q2 − q21 )ω2 + p20 − q20 = 0.
(2.12)
It is easy to show that p22
− 2p1 −
q22
= (b + r1 ) + r2 + 2b1 x2 + 2
p21 − 2p0 p2 + 2q0 q2 − q21 =
(1 + mx∗2 )2
[ (b + r1 ) r2 + 2b1 x∗2 + +
[
2
a1 y∗
∗
a1 r 2 y ∗
[
(1 + mx∗2 )2
p0 − q0 = r (b + r1 ) 2 r2 + 2b1 x∗2 −
ab b + r1
+ 2ab,
a1 y∗
2(r2 + 2b1 x2 ) +
+
− ab
(1 + mx∗2 )2 ∗
a1 y ∗
(1 + mx∗2 )
]2
(2.13)
a1 y ∗
]
(1 + mx∗2 )2 ] . 2
,
Hence, if p0 > q0 , Eq. (2.12) has no positive real roots. Accordingly, by Theorem 3.4.1 in Kuang [9], we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 and p0 > q0 hold, then E ∗ is locally asymptotically stable for all τ ≥ 0. If p0 < q0 , then Eq. (2.12) has a unique positive root ω0 , that is, Eq. (2.7) has a pair of purely imaginary roots of the form ±iω0 . Denote
τ0n =
1
ω0
arccos
q1 ω02 (ω02 − p1 ) + (q2 ω02 − q0 )(p0 − p2 ω02 ) q21 ω02 + (q2 ω02 − q0 )2
+
2nπ
ω0
,
n = 0, 1, 2, . . . .
By Theorem 3.4.1 in Kuang [9] we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 and p0 < q0 hold, then E ∗ remains stable for τ < τ0 := τ00 . We now claim that d(Reλ) dτ
|τ =τ0 > 0.
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2155
This will show that there exists at least one eigenvalue with a positive real part for τ > τ0 . Moreover, the conditions for the existence of a Hopf bifurcation [20] are then satisfied yielding a periodic solution. To this end, differentiating Eq. (2.7) with respect to τ , it follows that
dλ
3λ2 + 2p2 λ + p1
−1 =
dτ
−λ(λ + p2 λ + p1 λ + p0 ) 3
2
+
2q2 λ + q1
λ(q2 λ + q1 λ + q0 ) 2
−
τ . λ
Hence, a direct calculation shows that
sgn
d(Reλ)
= sgn Re
dτ
λ=iω0
dλ
−1
dτ
λ=iω0
(p1 − 3ω02 )(ω02 − p1 ) + 2p2 (p0 − p2 ω02 ) −q21 + 2q2 (q0 − q2 ω02 ) . = sgn − + (q0 − q2 ω02 )2 + q21 ω02 (ω03 − p1 ω0 )2 + (p0 − p2 ω02 )2
We derive from (2.11) that
(ω03 − p1 ω0 )2 + (p0 − p2 ω02 )2 = (q0 − q2 ω02 )2 + q21 ω02 . Hence, it follows that
sgn
d(Reλ)
dτ
= sgn
3ω04 + 2(p22 − 2p1 − q22 )ω02 + p21 − 2p0 p2 + 2q0 q2 − q21
(q0 − q2 ω02 )2 + q21 ω02
λ=iω0
> 0.
Therefore, the transversal condition holds and a Hopf bifurcation occurs at ω = ω0 , τ = τ0 . In conclusion, we have the following results. Theorem 2.1. For system (1.3), we have the following: (i) If ab < r2 (b + r1 ), then the trivial equilibrium E0 (0, 0, 0) is locally asymptotically stable; if ab > r2 (b + r1 ), then E0 is unstable. (ii) If (H1) and (H3) hold, then the predator-extinction equilibrium E1 (x01 , x02 , 0) is locally asymptotically stable; if (H1)–(H2) hold, then E1 is unstable. (iii) Let (H1)–(H2) and (p2 + q2 )(p1 + q1 ) > p0 + q0 hold. If p0 > q0 , then the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) is locally asymptotically stable for all τ ≥ 0; if p0 < q0 , then there exists a positive number τ0 such that E ∗ is locally asymptotically stable if 0 < τ < τ0 and is unstable if τ > τ0 . Further, system (1.3) undergoes a Hopf bifurcation at E ∗ when τ = τ0 . 3. Permanence In this section, we are concerned with the permanence of system (1.3). Definition 3.1. System (1.3) is said to be permanent (uniformly persistent) if there are positive constants m1 , m2 , M1 and M2 such that each positive solution of system (1.3) satisfies m1 ≤ lim inf xi (t ) ≤ lim sup xi (t ) ≤ M1 , t →+∞
t →+∞
i = 1, 2;
m2 ≤ lim inf y(t ) ≤ lim sup y(t ) ≤ M2 . t →+∞
t →+∞
In order to study the permanence of system (1.3), we refer to the persistence theory on infinite dimensional systems developed by Hale and Waltman [21]. Let X be a complete metric space with metric d. Suppose that T is a continuous semiflow on X , i.e., a continuous mapping T : [0, ∞) × X → X with the following properties Tt ◦ Ts = Tt +s ,
t , s ≥ 0,
T0 (x) = x,
x ∈ X,
where Tt denotes the mapping from X to X given by Tt (x) = T (t , x). The distance d(x, Y ) of a point x ∈ X from a subset Y of X is defined by d(x, Y ) = inf d(x, y). y∈Y
Recall that the positive orbit γ + (x) through x is defined as γ + (x) = ∪t ≥0 {T (t )x}, and its ω-limit set is ω(x) = ∩s≥0 ∪t ≥s {T (t )x}. Define W s (A) the strong stable set of a compact invariant set A as W s (A) = {x : x ∈ X , ω(x) ̸= φ, ω(x) ⊂ A}. (C1) Assume that X 0 is open and dense in X and X 0 ∪ X0 = X , X 0 ∩ X0 = φ . Moreover, the C 0 -semigroup T (t ) on X satisfies T (t ) : X 0 → X 0 ,
T (t ) : X0 → X0 .
Let Tb (t ) = T (t )|X0 and Ab be the global attractor for Tb (t ). Define A˜ b = ∪x∈Ab ω(x).
2156
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
Lemma 3.1 (Hale & Waltman [21]). Suppose that T (t ) satisfies (C1) and the following conditions: (i) There is a t0 ≥ 0 such that T (t ) is compact for t > t0 ; (ii) T (t ) is point dissipative in X ; where (iii) A˜ b is isolated and has an acyclic covering M,
= {M ˜ 1, M ˜ 2, . . . , M ˜ n }; M ˜ i ) ∩ X 0 = φ for i = 1, 2, . . . , n. (iv) W s (M Then X0 is a uniform repeller with respect to X 0 , that is, there is an ε > 0 such that for any x ∈ X 0 , lim inft →+∞ d(T (t )x, X0 ) ≥ ε . To study the permanence of system (1.3), we also need the following result. Lemma 3.2. There are positive constants M1 and M2 such that for any positive solution (x1 (t ), x2 (t ), y(t )) of system (1.3), lim sup x1 (t ) < M1 , t →+∞
lim sup x2 (t ) < M1 , t →+∞
lim sup y(t ) < M2 .
(3.1)
t →+∞
Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). Denote r0 = min{r1 , r2 , r }. Define a1 V (t ) = x1 (t − τ ) + x2 (t − τ ) + y(t ). a2 Calculating the derivative of V (t ) along positive solutions of system (1.3), it follows that d dt
V (t ) = ax2 (t − τ ) − r1 x1 (t − τ ) − r2 x2 (t − τ ) − b1 x22 (t − τ ) −
a1 r a2
y(t )
≤ −r0 V (t ) + ax2 (t − τ ) − b1 x22 (t − τ ) ≤ −r0 V (t ) +
a2 4b1
,
which yields lim supt →+∞ V (t ) ≤ a2 /(4b1 r0 ). If we choose M1 = a2 /(4b1 r0 ), M2 = a2 a2 /(4a1 b1 r0 ), then (3.1) follows. This completes the proof. We are now in a position to state and prove our result on the permanence of system (1.3). Theorem 3.1. If (H1)–(H2) hold, then system (1.3) is permanent. Proof. We need only to show that the boundaries of R3+0 repel positive solutions of system (1.3) uniformly. Let C + ([−τ , 0], R3+0 ) denote the space of continuous functions mapping [−τ , 0] into R3+0 . Define C1 = {(φ1 , φ2 , ψ) ∈ C + ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ 0, φ2 (θ ) ≡ 0, θ ∈ [−τ , 0]}, C2 = {(φ1 , φ2 , ψ) ∈ C + ([−τ , 0], R3+0 ) : φ1 (θ ) > 0, φ2 (θ ) > 0, ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}. Denote C0 = C1 ∪ C2 , X = C + ([−τ , 0], R3+0 ) and C 0 = intC + ([−τ , 0], R3+0 ). In the following, we verify that the conditions in Lemma 3.1 are satisfied. By the definition of C 0 and C0 , it is easy to see that C 0 and C0 are positively invariant and the condition (ii) in Lemma 3.1 is clearly satisfied. Noting that the functions in the right side of system (1.3) are in C 1 and the solution of system (1.3) with initial conditions (1.4) is ultimately bounded, using the smoothing property of solutions of delay differential equations introduced in Kuang [9] (Theorem 2.2.8), it follows that condition (i) in Lemma 3.1 is satisfied. Thus, we need only to show that the conditions (iii) and (iv) hold. Clearly, corresponding to x1 (t ) = x2 (t ) = y(t ) = 0 and x1 (t ) = x01 , x2 (t ) = x02 , y(t ) = 0, respectively, there are two constant solutions in C0 : E˜ 0 ∈ C1 , E˜ 1 ∈ C2 satisfying
E˜ 0 = {(φ1 , φ2 , ψ) ∈ ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ 0, φ2 (θ ) ≡ 0, ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}, E˜ 1 = {(φ1 , φ2 , ψ) ∈ ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ x01 , φ2 (θ ) ≡ x02 , ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}. We now verify the condition (iii) of Lemma 3.1. If (x1 (t ), x2 (t ), y(t )) is a solution of system (1.3) initiating from C1 , then y˙ (t ) = −ry(t ), which yields y(t ) → 0 as t → +∞. If (x1 (t ), x2 (t ), y(t )) is a solution of system (1.3) initiating from C2 with x1 (0) > 0, x2 (0) > 0, then we have x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) = bx1 (t ) − r2 x2 (t ) − b1 x22 (t ).
(3.2)
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2157
System (3.2) has two equilibria: A0 (0, 0) and A1 (x01 , x02 ), where x01 and x02 are defined as in (2.1). Define x1
V0 (t ) = x1 − x01 − x01 ln
x01
+
ax02
bx01
x2 − x02 − x02 ln
x2 x02
.
Calculating the derivative of V0 (t ) along positive solutions of system (3.2), it follows that d dt
V0 (t ) = −
ax02 bx01
b1 (x2 (t ) −
) −
x02 2
a
x2 (t ) x1 (t )
x01
(x1 (t ) −
x01
)−
x1 ( t ) x2 ( t )
2 (x2 (t ) −
x02
)
,
which yields xi (t ) → x0i (i = 1, 2) as t → +∞. Noting that C1 ∩ C2 = φ , it follows that the invariant sets E˜ 0 and E˜ 1 are isolated. Hence, {E˜ 0 , E˜ 1 } is isolated and is an acyclic covering satisfying the condition (iii) in Lemma 3.1. We now show that W s (E˜ 1 ) ∩ C 0 = φ , and W s (E˜ 1 ) ∩ C 0 = φ . Here, we only show the second equation holds since the proof of the first equation is simple. Assume W s (E˜ 1 ) ∩ C 0 ̸= φ . Then there is a positive solution (x1 (t ), x2 (t ), y(t )) with limt →+∞ (x1 (t ), x2 (t ), y(t )) = (x01 , x02 , 0). Noting that (H1)–(H2) hold, one can choose ε > 0 sufficiently small such that a2 (x02 − ε) 1 + m(x02 − ε)
> r.
(3.3)
Let t0 be sufficiently large such that, for t > t0 , x02 − ε < x2 (t ) < x02 + ε. For ε > 0 sufficiently small satisfying (3.3), it follows from the third equation of system (1.3) that, for t > t0 + τ , y˙ (t ) >
a2 (x02 − ε) 1 + m(x02 − ε)
y(t − τ ) − ry(t ),
which, together with (3.3), yields limt →+∞ y(t ) = +∞. This contradicts Lemma 3.2. Hence, we have W s (E˜ 1 ) ∩ C 0 = φ . By Lemma 3.1, we are now able to conclude that C0 repels positive solutions of system (1.3) uniformly. Hence, system (1.3) is permanent. The proof is complete. 4. Global stability In this section, we study the global stability of each of feasible equilibria of system (1.3), respectively. The strategy of proofs is to use Lyapunov functionals and the LaSalle invariant principle. Theorem 4.1. If ab < r2 (b + r1 ), then the trivial equilibrium E0 (0, 0, 0) of system (1.3) is globally asymptotically stable. Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). By Theorem 2.1, we see that if ab < r2 (b + r1 ), then E0 is locally asymptotically stable. Define b + r1
V11 (t ) = x1 +
b
x2 +
a1 ( b + r 1 ) a2 b
y.
(4.1)
Calculating the derivative of V11 (t ) along positive solutions of system (1.3), it follows that d dt
V11 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ) +
+ =
a1 (b + r1 )
bx1 (t ) − r2 x2 (t ) −
b
a2 x2 (t − τ )y(t − τ ) 1 + mx2 (t − τ )
a2 b
a−
+
[
[
b + r1
r2 (b + r1 )
b
a1 (b + r1 )
[
x2 (t ) −
b + r1 b
b1 x22 (t ) +
1 + mx2 (t − τ )
(t ) −
a1 x2 (t )y(t )
]
1 + mx2 (t )
]
[
a2 x2 (t − τ )y(t − τ )
a2 b
− ry(t )
b1 x22
a1 x2 (t )y(t )
]
1 + mx2 (t )
]
− ry(t ) .
(4.2)
Define V1 (t ) = V11 (t ) +
a1 ( b + r 1 ) b
∫
t t −τ
x2 (s)y(s) 1 + mx2 (s)
ds.
(4.3)
2158
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
We derive from (4.2) and (4.3) that d dt
V1 ( t ) =
r2 ( b + r1 )
a−
b
x2 (t ) −
b1 (b + r1 ) b
x22 (t ) −
a1 r (b + r1 ) a2 b
y(t ).
(4.4)
If ab < r2 (b + r1 ), it follows from (4.4) that V1′ (t ) ≤ 0. By Theorem 5.3.1 in [20], solutions limit to M , the largest invariant subset of {V1′ (t ) = 0}. Clearly, we see from (4.4) that V1′ (t ) = 0 if and only if x2 = 0, y = 0. Noting that M is invariant, for each element in M , we have x2 (t ) = 0, y(t ) = 0. It therefore follows from the second equation of system (1.3) that 0 = x˙ 2 (t ) = bx1 (t ), which yields x1 (t ) = 0. Hence, V1′ (t ) = 0 if and only if (x1 , x2 , y) = (0, 0, 0). Accordingly, the global asymptotic stability of E0 (0, 0, 0) follows from the LaSalle invariance principle. This completes the proof. Theorem 4.2. If (H1) and (H3) hold, then the predator-extinction equilibrium E1 (x01 , x02 , 0) of system (1.3) is globally asymptotically stable. Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). By Theorem 2.1, we see that if (H1) and (H3) hold, then E1 is locally asymptotically stable. Define V21 (t ) = x1 − x01 − x01 ln
x1 x01
x2 + c1 x2 − x02 − x02 ln 0 + c2 y,
(4.5)
x2
where c1 = ax02 /(bx01 ), c2 = c1 a1 (1 + mx02 )/a2 . Calculating the derivative of V21 (t ) along positive solutions of system (1.3), it follows that d dt
[ ] x0 a1 x2 (t )y(t ) [ax2 (t ) − r1 x1 (t ) − bx1 (t )] + c1 1 − 2 bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − x1 x2 1 + mx2 (t ) [ ] a2 x2 (t − τ )y(t − τ ) + c2 − ry(t ) 1 + mx2 (t − τ ) a x0 = 0 1 − 1 [−x2 (t )(x1 (t ) − x01 ) + x1 (t )(x2 (t ) − x02 )]
V21 (t ) =
1−
x01
x1
x1
+ c1
b
x02
1−
x02
x2
[−x1 (t )(x2 (t ) − x02 ) + x2 (t )(x1 (t ) − x01 )] − c1 b1 (x2 (t ) − x02 )2
[ ] a2 x2 (t − τ )y(t − τ ) x0 (1 + mx2 ) a1 x2 (t )y(t ) − c1 1 + mx02 − 2 + c2 − ry(t ) . x2 1 + mx2 (t ) 1 + mx2 (t − τ )
(4.6)
Define V2 (t ) = V21 (t ) + c2 a2
∫
x2 (s)y(s)
t t −τ
1 + mx2 (s)
ds.
(4.7)
We derive from (4.6) and (4.7) that d dt
V2 ( t ) = −
a x01
x2 (t ) x1 (t )
(x1 (t ) −
x01
− c1 b1 (x2 (t ) − x02 )2 + c1
)−
x1 (t ) x2 (t )
2 (x2 (t ) −
x02
)
a1 [(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 )] a2 b 1 ( b + r 1 )
y(t ).
(4.8)
If (H3) holds, it then follows from (4.8) that V2′ (t ) ≤ 0. By Theorem 5.3.1 in [20], solutions limit to M , the largest invariant subset of {V2′ (t ) = 0}. Clearly, we see from (4.8) that V2′ (t ) = 0 if and only if x2 = x02 , y = 0. Noting that M is invariant, for each element in M , we have x2 (t ) = x02 , y(t ) = 0. It follows from the second equation of system (1.3) that 0 = x˙ 2 (t ) = bx1 (t ) − r2 x02 − b1 (x02 )2 , which yields x1 (t ) = x01 . Hence, V2′ (t ) = 0 if and only if (x1 , x2 , y) = (x01 , x02 , 0). Accordingly, the global asymptotic stability of E1 follows from the LaSalle invariance principle. This completes the proof. We are now in a position to state and prove our result on the global stability of the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) of system (1.3).
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2159
Theorem 4.3. Let (H1)–(H2) hold. Then the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) of system (1.3) is globally asymptotically stable provided that (H4) x2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )]. Here, x > 0 is the persistency constant for x2 satisfying lim inft →+∞ x2 (t ) ≥ x2 . Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). Since x2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )], there is a T > 0 such that x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for all t ≥ T . Accordingly, we have x∗2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )]. In this case, it is easy to show that r2 + 2b1 x∗2 −
ab b + r1
> 0.
(4.9)
Direct calculations show that p1 + q1 =
a1 ry∗
+ (b + r1 ) r2 + 2b1 x∗2 −
(1 + mx∗2 )2
a1 y∗
∗
p2 + q2 = b + r1 + r2 + 2b1 x2 +
(1 + mx∗2 )2
ab b + r1
+
a1 y∗
(1 + mx∗2 )2
, (4.10)
.
Hence, if (H4) holds, it then follows from (2.10), (2.13) and (4.10) that p0 > q0 , (p1 + q1 )(p2 + q2 ) > p0 + q0 . By Theorem 2.1, E ∗ (x∗1 , x∗2 , y∗ ) is locally asymptotically stable for all τ ≥ 0. Define x1 V31 (t ) = x1 − x∗1 − x∗1 ln ∗ + k1 x1
x2 x2 − x∗2 − x∗2 ln ∗ x2
y + k2 y − y∗ − y∗ ln ∗ ,
(4.11)
y
where k1 = ax∗2 /(bx∗1 ), k2 = k1 a1 (1 + mx∗2 )/a2 . Calculating the derivative of V31 (t ) along positive solutions of system (1.3), it follows that d dt
V31 (t ) =
x∗1
[ax2 (t ) − r1 x1 (t ) − bx1 (t )] [ ] x∗2 a1 x2 (t )y(t ) 2 + k1 1 − bx1 (t ) − r2 x2 (t ) − b1 x2 (t ) − x2 1 + mx2 (t ) [ ] ∗ a2 x2 (t − τ )y(t − τ ) y − ry(t ) . + k2 1 − y 1 + mx2 (t − τ ) 1−
x1
(4.12)
Eq. (4.12) can be rewritten as d dt
V31 (t ) =
a
x∗1
b
x∗2
[−x2 (x1 − x1 ) + x1 (x2 − x2 )] + k1 ∗ 1 − [−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )] x2 x2 [ ∗ ] x∗2 x1 a1 x2 (t )y(t ) 2 + k1 1 − b ∗ − r2 x2 (t ) − b1 x2 (t ) − x2 x2 1 + mx2 (t ) ] ∗[ y a2 x2 (t − τ )y(t − τ ) + k2 1 − − ry(t ) . y 1 + mx2 (t − τ ) x∗1
1−
∗
∗
x1
(4.13)
On substituting bx∗1 − r2 x∗2 − b1 x∗2 2 = a1 x∗2 y∗ /(1 + mx∗2 ) into (4.13), we derive that d dt
V31 (t ) =
a
x∗1
b
x∗2
[−x2 (x1 − x1 ) + x1 (x2 − x2 )] + k1 ∗ 1 − [−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )] x2 x2 [ ∗ ∗ ] x∗2 x1 x1 ∗ ∗ + k1 1 − x2 b ∗ − r2 − b1 x2 − x2 b ∗ − r2 − b1 x2 x2 x2 x2 ∗ ∗ ∗ a1 x2 y x a1 x2 (t )y(t ) + k1 1 − 2 − k1 (1 + mx∗2 ) 1 + mx∗2 x2 1 + mx2 (t ) [ ] a2 x2 (t − τ )y(t − τ ) a2 y∗ x2 (t − τ )y(t − τ ) ∗ + k2 − + ry . 1 + mx2 (t − τ ) y(t )(1 + mx2 (t − τ )) x∗1
1−
∗
∗
x1
(4.14)
Define V3 (t ) = V31 (t ) + k2 a2
∫
t t −τ
[
x2 (s)y(s) 1 + mx2 (s)
−
x∗2 y∗ 1 + mx∗2
−
(1 + mx∗2 )x2 (s)y(s) ln ds. 1 + mx∗2 x∗2 y∗ (1 + mx2 (s)) x∗2 y∗
]
(4.15)
2160
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
It follows from (4.14) and (4.15) that d dt
x∗1
a
V3 ( t ) =
1−
x∗1
x1
[−x2 (x1 − x∗1 ) + x1 (x2 − x∗2 )] + k1
(x2 − x∗2 )2
b x∗2
1−
x∗2 x2
[−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )]
] a1 x∗2 y∗ x∗2 − r2 − b1 (x2 + x∗2 ) + k1 1 − x b + r1 1 + mx∗2 x2 ∗ ∗ ∗ x y (1 + mx2 )x2 (t − τ )y(t − τ ) − k2 a2 2 ∗ + k2 ry∗ 1 + mx2 x∗2 y(t )(1 + mx2 (t − τ )) + k1
+ k2 a2
x∗2 y∗ 1 + mx∗2
[
ab
ln
(1 + mx2 (t ))x2 (t − τ )y(t − τ ) . x2 (t )y(t )(1 + mx2 (t − τ ))
(4.16)
Noting that a ∗
x1
= k1
b
∗,
k2 ry∗ = k2 a2
x2
x∗2 y∗ 1 + mx∗2
= k1 a1 x∗2 y∗ ,
and a1 x∗2 y∗
1 + mx∗2
1−
x∗2
∗ ∗
= a1 x2 y
x2
1−
x∗2 (1 + mx2 )
(1 + mx∗2 )x2
,
we derive from (4.16) that d
a
V3 ( t ) = − ∗ dt x1
x2 (t ) x1 (t )
(x1 (t ) − x1 ) − ∗
x1 (t ) x2 (t )
(x2 (t ) − x2 )
+ k1
(x2 − x∗2 )2 x
[
ab b + r1
] − r2 − b1 (x2 + x∗2 )
x (1 + mx2 ) − 1 − ln 2 (1 + mx∗2 )x2 (1 + mx∗2 )x2 [ ] (1 + mx∗2 )x2 (t − τ )y(t − τ ) (1 + mx∗2 )x2 (t − τ )y(t − τ ) x∗ y∗ − k2 a2 2 ∗ − 1 − ln . 1 + mx2 x∗2 y(t )(1 + mx2 (t − τ )) x∗2 y(t )(1 + mx2 (t − τ )) ∗ ∗
− k2 a2
x2 y
[
x2 (1 + mx2 )
2 ∗
∗
∗
]
1 + mx∗2
If x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for t ≥ T , it follows that
(x2 − x∗2 )2 x2
[
ab b + r1
] − r 2 − b1 ( x2 + x2 ) ≤ 0 , ∗
with equality if and only if x2 = x∗2 . This implies that if x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for t ≥ T , V3′ (t ) ≤ 0, √ √ with equality if and only if x2 = x∗2 , x2 (t )/x1 (t )(x1 (t ) − x∗1 ) = x1 (t )/x2 (t )(x2 (t ) − x∗2 ), x∗2 (1 + mx2 )/(x2 (1 + mx∗2 )) = (1 + mx∗2 )x2 (t − τ )y(t − τ )/(x∗2 y(t )(1 + mx2 (t − τ ))) = 1. We now look for the invariant subset M within the set M =
(1 + mx∗2 )x2 (t − τ )y(t − τ ) x∗ (1 + mx2 ) (x1 , x2 , y) : x1 = x∗1 , x2 = x∗2 , 2 = = 1 . (1 + mx∗2 )x2 x∗2 y(t )(1 + mx2 (t − τ ))
Since x1 = x∗1 , x2 = x∗2 on M and consequently, 0 = x˙ 2 (t ) = bx∗1 − r2 x∗2 − b1 x∗2 2 − a1 x∗2 y(t )/(1 + mx∗2 ), which yields y(t ) = y∗ . Hence, the only invariant set in M is M = {(x∗1 , x∗2 , y∗ )}. Using LaSalle’s invariant principle, the global asymptotic stability of E ∗ follows. This completes the proof. 5. Numerical simulations In this section, we give two examples to illustrate the main results in Theorems 2.1 and 4.3. Example 1. In (1.3), let a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 16, r = 1/8, r1 = r2 = 1/8, m = 1/10. It is easy to show that (a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) = 12 621/1024 > 0. Hence, system (1.3) has a unique coexistence equilibrium E ∗ (160/597, 10/199, 105 175/118 803). By calculation, we have p0 − q0 ≈ −0.3963 < 0, (p2 + q2 )(p1 + q1 ) − (p0 + q0 ) ≈ 9.7196, τ0 ≈ 3.1489. By Theorem 2.1, E ∗ is locally asymptotically stable if 0 < τ < τ0 and is unstable if τ > τ0 , and system (1.3) undergoes a Hopf bifurcation at E ∗ when τ = τ0 . A numerical simulation illustrates this fact (see Fig. 1). Example 2. In (1.3), let a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 50, r = 1/8, r1 = 1/2, r2 = 1/8, m = 1/10. A direct calculation shows that (a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) ≈ 5.0836 > 0. By Theorem 3.1, system (1.3) is permanent. It is easy to show that system (1.3) has a unique coexistence equilibrium E ∗ (0.2010, 0.0503, 0.2739). From the proof of Lemma 3.2, we see that lim supt →+∞ y(t ) ≤ M2 := a2 a2 /(4a1 b1 r0 ) ≈ 0.18. Hence, for ε > 0 sufficiently small, there is a constant T > 0 such that if t > T , y(t ) < M2 + ε . We derive from the first and the second equations of system (1.3) that,
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2161
Fig. 1. The temporal solution found by numerical integration of system (1.3) with a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 16, r = 1/8, r1 = r2 = 1/8, m = 1/10; (a) τ = 2 and (b) τ = 3.5, respectively; (φ1 , φ2 , ψ) ≡ (0.1, 0.1, 0.1).
Fig. 2. The temporal solution found by numerical integration of system (1.3) with a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 50, r = 1/8, r1 = 1/2, r2 = 1/8, m = 1/10, τ = 3, (φ1 , φ2 , ψ) ≡ (0.01, 0.01, 0.01).
for t > T , x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) ≥ bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − a1 x2 (t )(M2 + ε),
(5.1)
which yields lim inf x2 (t ) ≥ x2 := t →+∞
ab − r2 (b + r1 ) − a1 (b + r1 )M2 b1 ( b + r 1 )
≈ 0.0595.
By calculation, we have (ab − r2 (b + r1 ))/(2b1 (b + r1 )) ≈ 0.0388 < x2 . By Theorem 4.3, the equilibrium E ∗ is globally asymptotically stable. A numerical simulation illustrates our result (see Fig. 2). 6. Discussion In this paper, the global dynamics of a stage-structured predator–prey model with Holling type-II functional response and time delay due to the gestation of the predator was investigated using Lyapunov functionals and LaSalle’s invariance principle. It has been shown that, under a hypothesis guaranteeing the uniform persistence of the system, a feasible lower bound condition on the density of the mature prey population ensures the global asymptotic stability of the coexistence equilibrium. That is, if the mature prey is always abundant enough, then the system tends to the unique coexistence equilibrium. On the other hand, it was shown that under some conditions, the time delay due to the gestation of the predator may destabilize the coexistence equilibrium of the system. References [1] C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Ent. Sec. Can. 45 (1965) 1–60. [2] L. Chen, Z. Jing, The existence and uniqueness of limit cycles in differential equations modelling the predator–prey interaction, Chin. Sci. Bull. 29 (9) (1984) 521–523.
2162
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
[3] X. Liu, L. Chen, Complex dynamics of Holling type II Lotka–Volterra predator–prey system with impulsive perturbations on the predator, Chaos, Solitons Fractals 16 (2003) 311–320. [4] R. Xu, M.A.J. Chaplain, F.A. Davidson, Periodic solutions for a predator–prey model with Holling-type functional response and time delays, Appl. Math. Comput. 161 (2005) 637–654. [5] Wonlyul Ko, Kimun Ryu, Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a prey refuge, J. Differential Equations 231 (2006) 534–550. [6] G. Jiang, Q. Lu, L. Qian, Complex dynamics of a Holling type II prey-predator system with state feedback control, Chaos, Solitons Fractals 31 (2007) 448–461. [7] J.M. Cushing, Integro-differential equations and delay models in population dynamics, in: Lecture Notes in Biomathematics, vol. 20, Springer-Verlag, Berlin, Heidelberg, New York, 1977. [8] K. Gopalsamy, Stability and Oscillations in Delay Differential Equations of Population Dynamics, Kluwer Academic, Dordrecht, Norwell, MA, 1992. [9] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, 1993. [10] M.S. Bartlett, On theoretical models for competitive and predatory biological systems, Biometrika 44 (1957) 27–42. [11] E. Beretta, Y. Kuang, Global analyses in some delayed ratio-dependent predator–prey systems, Nonlinear Anal. TMA 32 (1998) 381–408. [12] P.J. Wangersky, W.J. Cunningham, Time lag in prey-predator population models, Ecology 38 (1957) 136–139. [13] W.G. Aiello, H.I. Freedman, A time delay model of single species growth with stage structure, Math. Biosci. 101 (1990) 139–156. [14] X. Song, L. Chen, Optimal harvesting and stability for a two-species competitive system with stage structure, Math. Biosci. 170 (2001) 173–186. [15] W. Wang, L. Chen, A predator–prey system with stage structure for predator, Comput. Math. Appl. 33 (1997) 83–91. [16] W. Wang, Global dynamics of a population model with stage structure for predator, in: L. Chen, et al. (Eds.), Advanced Topics in Biomathematics, Proceeding of the international conference on mathematical biology, World Scientific Publishing Co. Pte. Ltd, 1997, pp. 253–257. [17] Y. Xiao, L. Chen, Global stability of a predator–prey system with stage structure for the predator, Acta Math. Sin. Engl. Ser. 20 (2004) 63–70. [18] R. Xu, M.A.J. Chaplain, F. Davidson, Persistence and global stability of a ratio-dependent predator–prey model with stage structure, Appl. Math. Comput. 158 (2004) 729–744. [19] X. Zhang, L. Chen, Avidan U. Neumann, The stage-structured predator–prey model and optimal havesting policy, Math. Biosci. 168 (2000) 201–210. [20] J.K. Hale, Theory of Functional Differential Equations, Springer, New York, 1976. [21] J. Hale, P. Waltman, Persistence in infinite-dimensional systems, SIAM J. Math. Anal. 20 (1989) 388–395.
Contents lists available at ScienceDirect
Nonlinear Analysis: Real World Applications journal homepage: www.elsevier.com/locate/nonrwa
Global dynamics of a predator–prey model with time delay and stage structure for the prey✩ Rui Xu ∗ Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, Shijiazhuang 050003, PR China
article
info
Article history: Received 11 August 2010 Accepted 29 December 2010 Keywords: Stage structure Time delay Hopf bifurcation LaSalle invariant principle Global stability
abstract A stage-structured predator–prey system with Holling type-II functional response and time delay due to the gestation of predator is investigated. By analyzing the characteristic equations, the local stability of each of feasible equilibria of the system is discussed and the existence of a Hopf bifurcation at the coexistence equilibrium is established. By means of the persistence theory on infinite dimensional systems, it is proven that the system is permanent if the coexistence equilibrium exists. By using Lyapunov functionals and LaSalle invariant principle, it is shown that the trivial equilibrium is globally stable when both the predator-extinction equilibrium and the coexistence equilibrium are not feasible, and that the predator-extinction equilibrium is globally asymptotically stable if the coexistence equilibrium does not exist, and sufficient conditions are derived for the global stability of the coexistence equilibrium. Numerical simulations are carried out to illustrate the main theoretical results. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction In [1], based on experiments, Holling suggested three different kinds of functional responses for different kinds of species to model the phenomena of predation, which made the standard Lotka–Volterra system more realistic. Letting x(t ) and y(t ) represent the densities of the prey and the predator at time t, respectively, a standard predator–prey model with Holling type functional response is of the form x˙ (t ) = xg (x) − φ(x)y, y˙ (t ) = eφ(x)y − c (y).
(1.1)
In (1.1), the functions g (x) and c (y) represent the specific growth rate of the prey in the absence of predation and the mortality rate of the predator, respectively; the function φ(x) is called ‘‘functional response’’ representing the prey consumption per unit time. The most popular functional response used in the modelling of predator–prey systems is Holling type-II with φ(x) = fx/(1 + mx) which takes into account the time a predator uses in handling the prey being captured. There has been a large body of work in the literature on the global dynamics of predator–prey systems with Holling type-II functional responses (see, for example, [2–6] and references cited therein). It is well-known that time delays of one type or another have been incorporated into mathematical models of population dynamics by many researchers. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the
✩ This work was supported by the National Natural Science Foundation of China (Nos. 11071254, 10671209).
∗ Corresponding address: Institute of Applied Mathematics, Shijiazhuang Mechanical Engineering College, No. 97 Heping West Road, Shijiazhuang 050003, Hebei Province, PR China. Tel.: +86 311 87994035; fax: +86 311 87994007. E-mail address: [email protected]. 1468-1218/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.nonrwa.2010.12.029
2152
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
population to fluctuate. We refer to the monographs of Cushing [7], Gopalsamy [8] and Kuang [9] for general delayed biological systems and to Bartlett [10], Beretta and Kuang [11] and Wangersky and Cunningham [12] and references cited therein for studies on delayed predator–prey systems. Time delay due to gestation is a common example, because generally the consumption of prey by the predator throughout its past history governs the present birth rate of the predator. Therefore, more realistic models of population interactions should take into account the effect of time delays. In [12], Wangersky and Cunningham proposed and studied the following non-Kolmogorov-type predator–prey model: x˙ (t ) = x(t )(r1 − ax(t ) − a1 y(t )),
(1.2)
y˙ (t ) = a2 x(t − τ )y(t − τ ) − r2 y(t ).
This model assumes that a duration of τ time units elapses when an individual prey is killed and the moment when the corresponding addition is made to the predator population. Stage-structure is a natural phenomenon and represents, for example, the division of a population into immature and mature individuals. As is common, the dynamics-eating habits, susceptibility to predators, etc. are often quite different in these two sub-populations. Hence, it is of ecological importance to investigate the effects of such a subdivision on the interaction of species. There has been a significant body of work on stage-structured population models in last two decades (see, for example, [13–19]). In [18], Xu et al. considered a ratio-dependent predator–prey model with stage structure for the prey. The prey individuals were classified as belonging either the immature or the mature and it was assumed that the immature prey is not at risk of being attacked by the predator. This seems reasonable for a number of mammals, where the immature prey concealed in the mountain cave, are raised by their parents; they do not necessarily go out for seeking food, the rate they are attacked by the predators can be ignored. By constructing appropriate Lyapunov functions, sufficient conditions are obtained for the global asymptotic stability of nonnegative equilibria of the model. Motivated by the work of Chen and Jing [2], Wangersky and Cunningham [12] and Xu et al. [18], in the present paper, we are concerned with the combined effects of stage structure for the prey and time delay due to the gestation of the predator on the global dynamics of a predator–prey model with Holling type-II functional response. To this end, we consider the following delay differential system x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) = bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − y˙ (t ) =
a2 x2 (t − τ )y(t − τ ) 1 + mx2 (t − τ )
a1 x2 (t )y(t ) 1 + mx2 (t )
,
(1.3)
− ry(t ).
In system (1.3), x1 (t ) and x2 (t ) represent the densities of the immature and the mature prey at time t, respectively; y(t ) represents the density of the predator at time t. The model is derived under the following assumptions. (A1) The prey population: the birth rate is proportional to the existing mature population with a proportionality a > 0; the death rate of the immature population and the transformation rate from immature individuals to mature individuals are proportional to the existing immature population with proportionality constants r1 > 0 and b > 0, respectively; b1 > 0 is the intra-specific competition rate of the mature population; the death rate of the mature population is proportional to the existing mature population with a proportionality r2 > 0. (A2) The predator population: the growth of the species obeys a Holling type-II functional response. The predators feed only on the mature prey. r > 0 is the death rate of the predator, a1 > 0 is the capturing rate of the predator, a2 /a1 > 0 is the conversion rate of nutrients into the reproduction of the predator; m > 0 is the half saturation rate of the predator; τ ≥ 0 is a constant delay due to the gestation of the predator, that is, mature adult predators can only contribute to the reproduction of the predators’ biomass. The delay τ is based on the assumption that the change rate of predators depends on the number of the mature prey and of the predators presented at some previous time. The initial conditions for system (1.3) take the form x1 (θ ) = φ1 (θ ), x2 (θ ) = φ2 (θ ), y(θ ) = ψ(θ ), φ1 (θ ) ≥ 0, φ2 (θ ) ≥ 0, ψ(θ ) ≥ 0, θ ∈ [−τ , 0), φ1 (0) > 0, φ2 (0) > 0, ψ(0) > 0,
(1.4)
where (φ1 (θ ), φ2 (θ ), ψ(θ )) ∈ C ([−τ , 0], R3+0 ), the Banach space of continuous functions mapping the interval [−τ , 0] into R3+0 , where R3+0 = {(x1 , x2 , x3 ) : xi ≥ 0, i = 1, 2, 3}. It is well-known by the fundamental theory of functional differential equations [20], that system (1.3) has a unique solution (x1 (t ), x2 (t ), y(t )) satisfying initial conditions (1.4). It is easy to show that all solutions of system (1.3) corresponding to initial conditions (1.4) are defined on [0, +∞) and remain positive for all t ≥ 0. The organization of this paper is as follows. In the next section, by analyzing the corresponding characteristic equations, the local stability of each of the feasible equilibria of system (1.3) is discussed and the existence of a Hopf bifurcation at the coexistence equilibrium is established. In Section 3, by using the persistence theory on infinite dimensional systems
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2153
developed in [21], we prove that system (1.3) is permanent if the coexistence equilibrium exists. In Section 4, by using novel Lyapunov functionals and the LaSalle invariant principle, we show that both the prey and the predator go to extinction if both the predator-extinction equilibrium and the coexistence equilibrium are not feasible, and that the predator species goes to extinction when the coexistence equilibrium does not exist, and sufficient conditions are obtained for the global stability of the coexistence equilibrium. Numerical simulations are carried out in Section 5 to illustrate the main theoretical results. A brief discussion is given in Section 6 to conclude this work. 2. Local stability and Hopf bifurcation In this section, we study the local stability of each of feasible equilibria of system (1.3) and the existence of a Hopf bifurcation at the coexistence equilibrium. System (1.3) always has a trivial equilibrium E0 (0, 0, 0). If the following holds: (H1) ab > r2 (b + r1 ), then (1.3) has a predator-extinction equilibrium E1 (x01 , x02 , 0), where x01 =
a[ab − r2 (b + r1 )] b1 (b + r1 )2
,
x02 =
ab − r2 (b + r1 ) b1 (b + r1 )
.
(2.1)
Further, if (H1) and the following hold: (H2) (a2 − mr )[ab − r2 (b + r1 )] > b1 r (b + r1 ), then system (1.3) has a unique coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ), where ar
,
r
, (b + r1 )(a2 − mr ) a2 − mr a2 {(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 )} y∗ = . a1 (b + r1 )(a2 − mr )2
x∗1 =
x∗2 =
(2.2)
The characteristic equation of system (1.3) at the trivial equilibrium E0 (0, 0, 0) is of the form
(λ + r )[λ2 + (b + r1 + r2 )λ + r2 (b + r1 ) − ab] = 0.
(2.3)
It is readily seen from Eq. (2.3) that if ab < r2 (b + r1 ), then E0 is locally asymptotically stable; if ab > r2 (b + r1 ), then E0 is unstable. The characteristic equation of system (1.3) at the equilibrium E1 (x01 , x02 , 0) takes the form
[λ + (b + r1 + r2 + 2
2b1 x02
[
)λ + ab − r2 (b + r1 )] λ + r −
a2 x02 1 + mx02
e
−λτ
]
= 0.
(2.4)
Clearly, if (H1) holds, then the equation
λ2 + (b + r1 + r2 + 2b1 x02 )λ + ab − r2 (b + r1 ) = 0
(2.5)
has two roots with negative real parts. Other roots of Eq. (2.4) are determined by the equation
λ+r −
a2 x02 1 + mx02
e−λτ = 0.
(2.6)
Denote f (λ) = λ + r −
a2 x02 1 + mx02
e−λτ .
If (H1)–(H2) hold, it is easy to show that, for λ real, f (0) = −
(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) < 0, b1 (b + r1 )(1 + mx02 )
lim f (λ) = +∞.
λ→+∞
Hence, f (λ) = 0 has at least one positive real root. Therefore, if (H1)–(H2) hold, then the equilibrium E1 is unstable. If (H1) and the following hold: (H3) (a2 − mr )[ab − r2 (b + r1 )] < b1 r (b + r1 ),
2154
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
then it is easy to show that r − a2 x02 /(1 + mx02 ) > 0. By Theorem 3.2.1 in Kuang [9], we know that the equilibrium E1 is locally asymptotically stable for all τ ≥ 0. The characteristic equation of system (1.3) at the coexistence equilibrium E ∗ is of the form
λ3 + p2 λ2 + p1 λ + p0 + (q2 λ2 + q1 λ + q0 )e−λτ = 0,
(2.7)
where
[
p0 = r (b + r1 ) r2 + 2b1 x2 + ∗
p1 = r
∗
b + r1 + r2 + 2b1 x2 +
a1 y∗
(1 + mx∗2 )2 ∗
− ab ,
+ (b + r1 ) r2 + 2b1 x2 +
(1 + mx∗2 )2
(1 + mx∗2 )2 ∗ q0 = −r [(b + r1 )(r2 + 2b1 x2 ) − ab], q1 = −r (b + r1 + r2 + 2b1 x∗2 ), q2 = −r .
]
a1 y
a1 y ∗
p2 = b + r1 + r2 + 2b1 x∗2 +
a1 y∗
∗
(1 + mx∗2 )2
− ab, (2.8)
+ r,
When τ = 0, Eq. (2.7) becomes
λ3 + (p2 + q2 )λ2 + (p1 + q1 )λ + p0 + q0 = 0.
(2.9)
Clearly, p2 + q2 > 0. By calculation we derive that p0 + q0 =
a1 r (b + r1 )y∗
(1 + mx∗2 )2
> 0.
(2.10)
Hence, by the Routh–Hurwitz criterion, we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 , then the equilibrium E ∗ is locally asymptotically stable when τ = 0. If iω (ω > 0) is a solution of Eq. (2.7), separating real and imaginary parts, we have
ω3 − p1 ω = (q2 ω2 − q0 ) sin ωτ + q1 ω cos ωτ ,
(2.11)
−p2 ω2 + p0 = (q2 ω2 − q0 ) cos ωτ − q1 ω sin ωτ . Squaring and adding the two equations of (2.11), it follows that
ω6 + (p22 − 2p1 − q22 )ω4 + (p21 − 2p0 p2 + 2q0 q2 − q21 )ω2 + p20 − q20 = 0.
(2.12)
It is easy to show that p22
− 2p1 −
q22
= (b + r1 ) + r2 + 2b1 x2 + 2
p21 − 2p0 p2 + 2q0 q2 − q21 =
(1 + mx∗2 )2
[ (b + r1 ) r2 + 2b1 x∗2 + +
[
2
a1 y∗
∗
a1 r 2 y ∗
[
(1 + mx∗2 )2
p0 − q0 = r (b + r1 ) 2 r2 + 2b1 x∗2 −
ab b + r1
+ 2ab,
a1 y∗
2(r2 + 2b1 x2 ) +
+
− ab
(1 + mx∗2 )2 ∗
a1 y ∗
(1 + mx∗2 )
]2
(2.13)
a1 y ∗
]
(1 + mx∗2 )2 ] . 2
,
Hence, if p0 > q0 , Eq. (2.12) has no positive real roots. Accordingly, by Theorem 3.4.1 in Kuang [9], we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 and p0 > q0 hold, then E ∗ is locally asymptotically stable for all τ ≥ 0. If p0 < q0 , then Eq. (2.12) has a unique positive root ω0 , that is, Eq. (2.7) has a pair of purely imaginary roots of the form ±iω0 . Denote
τ0n =
1
ω0
arccos
q1 ω02 (ω02 − p1 ) + (q2 ω02 − q0 )(p0 − p2 ω02 ) q21 ω02 + (q2 ω02 − q0 )2
+
2nπ
ω0
,
n = 0, 1, 2, . . . .
By Theorem 3.4.1 in Kuang [9] we see that if (p2 + q2 )(p1 + q1 ) > p0 + q0 and p0 < q0 hold, then E ∗ remains stable for τ < τ0 := τ00 . We now claim that d(Reλ) dτ
|τ =τ0 > 0.
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2155
This will show that there exists at least one eigenvalue with a positive real part for τ > τ0 . Moreover, the conditions for the existence of a Hopf bifurcation [20] are then satisfied yielding a periodic solution. To this end, differentiating Eq. (2.7) with respect to τ , it follows that
dλ
3λ2 + 2p2 λ + p1
−1 =
dτ
−λ(λ + p2 λ + p1 λ + p0 ) 3
2
+
2q2 λ + q1
λ(q2 λ + q1 λ + q0 ) 2
−
τ . λ
Hence, a direct calculation shows that
sgn
d(Reλ)
= sgn Re
dτ
λ=iω0
dλ
−1
dτ
λ=iω0
(p1 − 3ω02 )(ω02 − p1 ) + 2p2 (p0 − p2 ω02 ) −q21 + 2q2 (q0 − q2 ω02 ) . = sgn − + (q0 − q2 ω02 )2 + q21 ω02 (ω03 − p1 ω0 )2 + (p0 − p2 ω02 )2
We derive from (2.11) that
(ω03 − p1 ω0 )2 + (p0 − p2 ω02 )2 = (q0 − q2 ω02 )2 + q21 ω02 . Hence, it follows that
sgn
d(Reλ)
dτ
= sgn
3ω04 + 2(p22 − 2p1 − q22 )ω02 + p21 − 2p0 p2 + 2q0 q2 − q21
(q0 − q2 ω02 )2 + q21 ω02
λ=iω0
> 0.
Therefore, the transversal condition holds and a Hopf bifurcation occurs at ω = ω0 , τ = τ0 . In conclusion, we have the following results. Theorem 2.1. For system (1.3), we have the following: (i) If ab < r2 (b + r1 ), then the trivial equilibrium E0 (0, 0, 0) is locally asymptotically stable; if ab > r2 (b + r1 ), then E0 is unstable. (ii) If (H1) and (H3) hold, then the predator-extinction equilibrium E1 (x01 , x02 , 0) is locally asymptotically stable; if (H1)–(H2) hold, then E1 is unstable. (iii) Let (H1)–(H2) and (p2 + q2 )(p1 + q1 ) > p0 + q0 hold. If p0 > q0 , then the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) is locally asymptotically stable for all τ ≥ 0; if p0 < q0 , then there exists a positive number τ0 such that E ∗ is locally asymptotically stable if 0 < τ < τ0 and is unstable if τ > τ0 . Further, system (1.3) undergoes a Hopf bifurcation at E ∗ when τ = τ0 . 3. Permanence In this section, we are concerned with the permanence of system (1.3). Definition 3.1. System (1.3) is said to be permanent (uniformly persistent) if there are positive constants m1 , m2 , M1 and M2 such that each positive solution of system (1.3) satisfies m1 ≤ lim inf xi (t ) ≤ lim sup xi (t ) ≤ M1 , t →+∞
t →+∞
i = 1, 2;
m2 ≤ lim inf y(t ) ≤ lim sup y(t ) ≤ M2 . t →+∞
t →+∞
In order to study the permanence of system (1.3), we refer to the persistence theory on infinite dimensional systems developed by Hale and Waltman [21]. Let X be a complete metric space with metric d. Suppose that T is a continuous semiflow on X , i.e., a continuous mapping T : [0, ∞) × X → X with the following properties Tt ◦ Ts = Tt +s ,
t , s ≥ 0,
T0 (x) = x,
x ∈ X,
where Tt denotes the mapping from X to X given by Tt (x) = T (t , x). The distance d(x, Y ) of a point x ∈ X from a subset Y of X is defined by d(x, Y ) = inf d(x, y). y∈Y
Recall that the positive orbit γ + (x) through x is defined as γ + (x) = ∪t ≥0 {T (t )x}, and its ω-limit set is ω(x) = ∩s≥0 ∪t ≥s {T (t )x}. Define W s (A) the strong stable set of a compact invariant set A as W s (A) = {x : x ∈ X , ω(x) ̸= φ, ω(x) ⊂ A}. (C1) Assume that X 0 is open and dense in X and X 0 ∪ X0 = X , X 0 ∩ X0 = φ . Moreover, the C 0 -semigroup T (t ) on X satisfies T (t ) : X 0 → X 0 ,
T (t ) : X0 → X0 .
Let Tb (t ) = T (t )|X0 and Ab be the global attractor for Tb (t ). Define A˜ b = ∪x∈Ab ω(x).
2156
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
Lemma 3.1 (Hale & Waltman [21]). Suppose that T (t ) satisfies (C1) and the following conditions: (i) There is a t0 ≥ 0 such that T (t ) is compact for t > t0 ; (ii) T (t ) is point dissipative in X ; where (iii) A˜ b is isolated and has an acyclic covering M,
= {M ˜ 1, M ˜ 2, . . . , M ˜ n }; M ˜ i ) ∩ X 0 = φ for i = 1, 2, . . . , n. (iv) W s (M Then X0 is a uniform repeller with respect to X 0 , that is, there is an ε > 0 such that for any x ∈ X 0 , lim inft →+∞ d(T (t )x, X0 ) ≥ ε . To study the permanence of system (1.3), we also need the following result. Lemma 3.2. There are positive constants M1 and M2 such that for any positive solution (x1 (t ), x2 (t ), y(t )) of system (1.3), lim sup x1 (t ) < M1 , t →+∞
lim sup x2 (t ) < M1 , t →+∞
lim sup y(t ) < M2 .
(3.1)
t →+∞
Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). Denote r0 = min{r1 , r2 , r }. Define a1 V (t ) = x1 (t − τ ) + x2 (t − τ ) + y(t ). a2 Calculating the derivative of V (t ) along positive solutions of system (1.3), it follows that d dt
V (t ) = ax2 (t − τ ) − r1 x1 (t − τ ) − r2 x2 (t − τ ) − b1 x22 (t − τ ) −
a1 r a2
y(t )
≤ −r0 V (t ) + ax2 (t − τ ) − b1 x22 (t − τ ) ≤ −r0 V (t ) +
a2 4b1
,
which yields lim supt →+∞ V (t ) ≤ a2 /(4b1 r0 ). If we choose M1 = a2 /(4b1 r0 ), M2 = a2 a2 /(4a1 b1 r0 ), then (3.1) follows. This completes the proof. We are now in a position to state and prove our result on the permanence of system (1.3). Theorem 3.1. If (H1)–(H2) hold, then system (1.3) is permanent. Proof. We need only to show that the boundaries of R3+0 repel positive solutions of system (1.3) uniformly. Let C + ([−τ , 0], R3+0 ) denote the space of continuous functions mapping [−τ , 0] into R3+0 . Define C1 = {(φ1 , φ2 , ψ) ∈ C + ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ 0, φ2 (θ ) ≡ 0, θ ∈ [−τ , 0]}, C2 = {(φ1 , φ2 , ψ) ∈ C + ([−τ , 0], R3+0 ) : φ1 (θ ) > 0, φ2 (θ ) > 0, ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}. Denote C0 = C1 ∪ C2 , X = C + ([−τ , 0], R3+0 ) and C 0 = intC + ([−τ , 0], R3+0 ). In the following, we verify that the conditions in Lemma 3.1 are satisfied. By the definition of C 0 and C0 , it is easy to see that C 0 and C0 are positively invariant and the condition (ii) in Lemma 3.1 is clearly satisfied. Noting that the functions in the right side of system (1.3) are in C 1 and the solution of system (1.3) with initial conditions (1.4) is ultimately bounded, using the smoothing property of solutions of delay differential equations introduced in Kuang [9] (Theorem 2.2.8), it follows that condition (i) in Lemma 3.1 is satisfied. Thus, we need only to show that the conditions (iii) and (iv) hold. Clearly, corresponding to x1 (t ) = x2 (t ) = y(t ) = 0 and x1 (t ) = x01 , x2 (t ) = x02 , y(t ) = 0, respectively, there are two constant solutions in C0 : E˜ 0 ∈ C1 , E˜ 1 ∈ C2 satisfying
E˜ 0 = {(φ1 , φ2 , ψ) ∈ ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ 0, φ2 (θ ) ≡ 0, ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}, E˜ 1 = {(φ1 , φ2 , ψ) ∈ ([−τ , 0], R3+0 ) : φ1 (θ ) ≡ x01 , φ2 (θ ) ≡ x02 , ψ(θ ) ≡ 0, θ ∈ [−τ , 0]}. We now verify the condition (iii) of Lemma 3.1. If (x1 (t ), x2 (t ), y(t )) is a solution of system (1.3) initiating from C1 , then y˙ (t ) = −ry(t ), which yields y(t ) → 0 as t → +∞. If (x1 (t ), x2 (t ), y(t )) is a solution of system (1.3) initiating from C2 with x1 (0) > 0, x2 (0) > 0, then we have x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) = bx1 (t ) − r2 x2 (t ) − b1 x22 (t ).
(3.2)
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2157
System (3.2) has two equilibria: A0 (0, 0) and A1 (x01 , x02 ), where x01 and x02 are defined as in (2.1). Define x1
V0 (t ) = x1 − x01 − x01 ln
x01
+
ax02
bx01
x2 − x02 − x02 ln
x2 x02
.
Calculating the derivative of V0 (t ) along positive solutions of system (3.2), it follows that d dt
V0 (t ) = −
ax02 bx01
b1 (x2 (t ) −
) −
x02 2
a
x2 (t ) x1 (t )
x01
(x1 (t ) −
x01
)−
x1 ( t ) x2 ( t )
2 (x2 (t ) −
x02
)
,
which yields xi (t ) → x0i (i = 1, 2) as t → +∞. Noting that C1 ∩ C2 = φ , it follows that the invariant sets E˜ 0 and E˜ 1 are isolated. Hence, {E˜ 0 , E˜ 1 } is isolated and is an acyclic covering satisfying the condition (iii) in Lemma 3.1. We now show that W s (E˜ 1 ) ∩ C 0 = φ , and W s (E˜ 1 ) ∩ C 0 = φ . Here, we only show the second equation holds since the proof of the first equation is simple. Assume W s (E˜ 1 ) ∩ C 0 ̸= φ . Then there is a positive solution (x1 (t ), x2 (t ), y(t )) with limt →+∞ (x1 (t ), x2 (t ), y(t )) = (x01 , x02 , 0). Noting that (H1)–(H2) hold, one can choose ε > 0 sufficiently small such that a2 (x02 − ε) 1 + m(x02 − ε)
> r.
(3.3)
Let t0 be sufficiently large such that, for t > t0 , x02 − ε < x2 (t ) < x02 + ε. For ε > 0 sufficiently small satisfying (3.3), it follows from the third equation of system (1.3) that, for t > t0 + τ , y˙ (t ) >
a2 (x02 − ε) 1 + m(x02 − ε)
y(t − τ ) − ry(t ),
which, together with (3.3), yields limt →+∞ y(t ) = +∞. This contradicts Lemma 3.2. Hence, we have W s (E˜ 1 ) ∩ C 0 = φ . By Lemma 3.1, we are now able to conclude that C0 repels positive solutions of system (1.3) uniformly. Hence, system (1.3) is permanent. The proof is complete. 4. Global stability In this section, we study the global stability of each of feasible equilibria of system (1.3), respectively. The strategy of proofs is to use Lyapunov functionals and the LaSalle invariant principle. Theorem 4.1. If ab < r2 (b + r1 ), then the trivial equilibrium E0 (0, 0, 0) of system (1.3) is globally asymptotically stable. Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). By Theorem 2.1, we see that if ab < r2 (b + r1 ), then E0 is locally asymptotically stable. Define b + r1
V11 (t ) = x1 +
b
x2 +
a1 ( b + r 1 ) a2 b
y.
(4.1)
Calculating the derivative of V11 (t ) along positive solutions of system (1.3), it follows that d dt
V11 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ) +
+ =
a1 (b + r1 )
bx1 (t ) − r2 x2 (t ) −
b
a2 x2 (t − τ )y(t − τ ) 1 + mx2 (t − τ )
a2 b
a−
+
[
[
b + r1
r2 (b + r1 )
b
a1 (b + r1 )
[
x2 (t ) −
b + r1 b
b1 x22 (t ) +
1 + mx2 (t − τ )
(t ) −
a1 x2 (t )y(t )
]
1 + mx2 (t )
]
[
a2 x2 (t − τ )y(t − τ )
a2 b
− ry(t )
b1 x22
a1 x2 (t )y(t )
]
1 + mx2 (t )
]
− ry(t ) .
(4.2)
Define V1 (t ) = V11 (t ) +
a1 ( b + r 1 ) b
∫
t t −τ
x2 (s)y(s) 1 + mx2 (s)
ds.
(4.3)
2158
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
We derive from (4.2) and (4.3) that d dt
V1 ( t ) =
r2 ( b + r1 )
a−
b
x2 (t ) −
b1 (b + r1 ) b
x22 (t ) −
a1 r (b + r1 ) a2 b
y(t ).
(4.4)
If ab < r2 (b + r1 ), it follows from (4.4) that V1′ (t ) ≤ 0. By Theorem 5.3.1 in [20], solutions limit to M , the largest invariant subset of {V1′ (t ) = 0}. Clearly, we see from (4.4) that V1′ (t ) = 0 if and only if x2 = 0, y = 0. Noting that M is invariant, for each element in M , we have x2 (t ) = 0, y(t ) = 0. It therefore follows from the second equation of system (1.3) that 0 = x˙ 2 (t ) = bx1 (t ), which yields x1 (t ) = 0. Hence, V1′ (t ) = 0 if and only if (x1 , x2 , y) = (0, 0, 0). Accordingly, the global asymptotic stability of E0 (0, 0, 0) follows from the LaSalle invariance principle. This completes the proof. Theorem 4.2. If (H1) and (H3) hold, then the predator-extinction equilibrium E1 (x01 , x02 , 0) of system (1.3) is globally asymptotically stable. Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). By Theorem 2.1, we see that if (H1) and (H3) hold, then E1 is locally asymptotically stable. Define V21 (t ) = x1 − x01 − x01 ln
x1 x01
x2 + c1 x2 − x02 − x02 ln 0 + c2 y,
(4.5)
x2
where c1 = ax02 /(bx01 ), c2 = c1 a1 (1 + mx02 )/a2 . Calculating the derivative of V21 (t ) along positive solutions of system (1.3), it follows that d dt
[ ] x0 a1 x2 (t )y(t ) [ax2 (t ) − r1 x1 (t ) − bx1 (t )] + c1 1 − 2 bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − x1 x2 1 + mx2 (t ) [ ] a2 x2 (t − τ )y(t − τ ) + c2 − ry(t ) 1 + mx2 (t − τ ) a x0 = 0 1 − 1 [−x2 (t )(x1 (t ) − x01 ) + x1 (t )(x2 (t ) − x02 )]
V21 (t ) =
1−
x01
x1
x1
+ c1
b
x02
1−
x02
x2
[−x1 (t )(x2 (t ) − x02 ) + x2 (t )(x1 (t ) − x01 )] − c1 b1 (x2 (t ) − x02 )2
[ ] a2 x2 (t − τ )y(t − τ ) x0 (1 + mx2 ) a1 x2 (t )y(t ) − c1 1 + mx02 − 2 + c2 − ry(t ) . x2 1 + mx2 (t ) 1 + mx2 (t − τ )
(4.6)
Define V2 (t ) = V21 (t ) + c2 a2
∫
x2 (s)y(s)
t t −τ
1 + mx2 (s)
ds.
(4.7)
We derive from (4.6) and (4.7) that d dt
V2 ( t ) = −
a x01
x2 (t ) x1 (t )
(x1 (t ) −
x01
− c1 b1 (x2 (t ) − x02 )2 + c1
)−
x1 (t ) x2 (t )
2 (x2 (t ) −
x02
)
a1 [(a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 )] a2 b 1 ( b + r 1 )
y(t ).
(4.8)
If (H3) holds, it then follows from (4.8) that V2′ (t ) ≤ 0. By Theorem 5.3.1 in [20], solutions limit to M , the largest invariant subset of {V2′ (t ) = 0}. Clearly, we see from (4.8) that V2′ (t ) = 0 if and only if x2 = x02 , y = 0. Noting that M is invariant, for each element in M , we have x2 (t ) = x02 , y(t ) = 0. It follows from the second equation of system (1.3) that 0 = x˙ 2 (t ) = bx1 (t ) − r2 x02 − b1 (x02 )2 , which yields x1 (t ) = x01 . Hence, V2′ (t ) = 0 if and only if (x1 , x2 , y) = (x01 , x02 , 0). Accordingly, the global asymptotic stability of E1 follows from the LaSalle invariance principle. This completes the proof. We are now in a position to state and prove our result on the global stability of the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) of system (1.3).
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2159
Theorem 4.3. Let (H1)–(H2) hold. Then the coexistence equilibrium E ∗ (x∗1 , x∗2 , y∗ ) of system (1.3) is globally asymptotically stable provided that (H4) x2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )]. Here, x > 0 is the persistency constant for x2 satisfying lim inft →+∞ x2 (t ) ≥ x2 . Proof. Let (x1 (t ), x2 (t ), y(t )) be any positive solution of system (1.3) with initial conditions (1.4). Since x2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )], there is a T > 0 such that x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for all t ≥ T . Accordingly, we have x∗2 > [ab − r2 (b + r1 )]/[2b1 (b + r1 )]. In this case, it is easy to show that r2 + 2b1 x∗2 −
ab b + r1
> 0.
(4.9)
Direct calculations show that p1 + q1 =
a1 ry∗
+ (b + r1 ) r2 + 2b1 x∗2 −
(1 + mx∗2 )2
a1 y∗
∗
p2 + q2 = b + r1 + r2 + 2b1 x2 +
(1 + mx∗2 )2
ab b + r1
+
a1 y∗
(1 + mx∗2 )2
, (4.10)
.
Hence, if (H4) holds, it then follows from (2.10), (2.13) and (4.10) that p0 > q0 , (p1 + q1 )(p2 + q2 ) > p0 + q0 . By Theorem 2.1, E ∗ (x∗1 , x∗2 , y∗ ) is locally asymptotically stable for all τ ≥ 0. Define x1 V31 (t ) = x1 − x∗1 − x∗1 ln ∗ + k1 x1
x2 x2 − x∗2 − x∗2 ln ∗ x2
y + k2 y − y∗ − y∗ ln ∗ ,
(4.11)
y
where k1 = ax∗2 /(bx∗1 ), k2 = k1 a1 (1 + mx∗2 )/a2 . Calculating the derivative of V31 (t ) along positive solutions of system (1.3), it follows that d dt
V31 (t ) =
x∗1
[ax2 (t ) − r1 x1 (t ) − bx1 (t )] [ ] x∗2 a1 x2 (t )y(t ) 2 + k1 1 − bx1 (t ) − r2 x2 (t ) − b1 x2 (t ) − x2 1 + mx2 (t ) [ ] ∗ a2 x2 (t − τ )y(t − τ ) y − ry(t ) . + k2 1 − y 1 + mx2 (t − τ ) 1−
x1
(4.12)
Eq. (4.12) can be rewritten as d dt
V31 (t ) =
a
x∗1
b
x∗2
[−x2 (x1 − x1 ) + x1 (x2 − x2 )] + k1 ∗ 1 − [−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )] x2 x2 [ ∗ ] x∗2 x1 a1 x2 (t )y(t ) 2 + k1 1 − b ∗ − r2 x2 (t ) − b1 x2 (t ) − x2 x2 1 + mx2 (t ) ] ∗[ y a2 x2 (t − τ )y(t − τ ) + k2 1 − − ry(t ) . y 1 + mx2 (t − τ ) x∗1
1−
∗
∗
x1
(4.13)
On substituting bx∗1 − r2 x∗2 − b1 x∗2 2 = a1 x∗2 y∗ /(1 + mx∗2 ) into (4.13), we derive that d dt
V31 (t ) =
a
x∗1
b
x∗2
[−x2 (x1 − x1 ) + x1 (x2 − x2 )] + k1 ∗ 1 − [−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )] x2 x2 [ ∗ ∗ ] x∗2 x1 x1 ∗ ∗ + k1 1 − x2 b ∗ − r2 − b1 x2 − x2 b ∗ − r2 − b1 x2 x2 x2 x2 ∗ ∗ ∗ a1 x2 y x a1 x2 (t )y(t ) + k1 1 − 2 − k1 (1 + mx∗2 ) 1 + mx∗2 x2 1 + mx2 (t ) [ ] a2 x2 (t − τ )y(t − τ ) a2 y∗ x2 (t − τ )y(t − τ ) ∗ + k2 − + ry . 1 + mx2 (t − τ ) y(t )(1 + mx2 (t − τ )) x∗1
1−
∗
∗
x1
(4.14)
Define V3 (t ) = V31 (t ) + k2 a2
∫
t t −τ
[
x2 (s)y(s) 1 + mx2 (s)
−
x∗2 y∗ 1 + mx∗2
−
(1 + mx∗2 )x2 (s)y(s) ln ds. 1 + mx∗2 x∗2 y∗ (1 + mx2 (s)) x∗2 y∗
]
(4.15)
2160
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
It follows from (4.14) and (4.15) that d dt
x∗1
a
V3 ( t ) =
1−
x∗1
x1
[−x2 (x1 − x∗1 ) + x1 (x2 − x∗2 )] + k1
(x2 − x∗2 )2
b x∗2
1−
x∗2 x2
[−x1 (x2 − x∗2 ) + x2 (x1 − x∗1 )]
] a1 x∗2 y∗ x∗2 − r2 − b1 (x2 + x∗2 ) + k1 1 − x b + r1 1 + mx∗2 x2 ∗ ∗ ∗ x y (1 + mx2 )x2 (t − τ )y(t − τ ) − k2 a2 2 ∗ + k2 ry∗ 1 + mx2 x∗2 y(t )(1 + mx2 (t − τ )) + k1
+ k2 a2
x∗2 y∗ 1 + mx∗2
[
ab
ln
(1 + mx2 (t ))x2 (t − τ )y(t − τ ) . x2 (t )y(t )(1 + mx2 (t − τ ))
(4.16)
Noting that a ∗
x1
= k1
b
∗,
k2 ry∗ = k2 a2
x2
x∗2 y∗ 1 + mx∗2
= k1 a1 x∗2 y∗ ,
and a1 x∗2 y∗
1 + mx∗2
1−
x∗2
∗ ∗
= a1 x2 y
x2
1−
x∗2 (1 + mx2 )
(1 + mx∗2 )x2
,
we derive from (4.16) that d
a
V3 ( t ) = − ∗ dt x1
x2 (t ) x1 (t )
(x1 (t ) − x1 ) − ∗
x1 (t ) x2 (t )
(x2 (t ) − x2 )
+ k1
(x2 − x∗2 )2 x
[
ab b + r1
] − r2 − b1 (x2 + x∗2 )
x (1 + mx2 ) − 1 − ln 2 (1 + mx∗2 )x2 (1 + mx∗2 )x2 [ ] (1 + mx∗2 )x2 (t − τ )y(t − τ ) (1 + mx∗2 )x2 (t − τ )y(t − τ ) x∗ y∗ − k2 a2 2 ∗ − 1 − ln . 1 + mx2 x∗2 y(t )(1 + mx2 (t − τ )) x∗2 y(t )(1 + mx2 (t − τ )) ∗ ∗
− k2 a2
x2 y
[
x2 (1 + mx2 )
2 ∗
∗
∗
]
1 + mx∗2
If x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for t ≥ T , it follows that
(x2 − x∗2 )2 x2
[
ab b + r1
] − r 2 − b1 ( x2 + x2 ) ≤ 0 , ∗
with equality if and only if x2 = x∗2 . This implies that if x2 (t ) > [ab − r2 (b + r1 )]/[2b1 (b + r1 )] for t ≥ T , V3′ (t ) ≤ 0, √ √ with equality if and only if x2 = x∗2 , x2 (t )/x1 (t )(x1 (t ) − x∗1 ) = x1 (t )/x2 (t )(x2 (t ) − x∗2 ), x∗2 (1 + mx2 )/(x2 (1 + mx∗2 )) = (1 + mx∗2 )x2 (t − τ )y(t − τ )/(x∗2 y(t )(1 + mx2 (t − τ ))) = 1. We now look for the invariant subset M within the set M =
(1 + mx∗2 )x2 (t − τ )y(t − τ ) x∗ (1 + mx2 ) (x1 , x2 , y) : x1 = x∗1 , x2 = x∗2 , 2 = = 1 . (1 + mx∗2 )x2 x∗2 y(t )(1 + mx2 (t − τ ))
Since x1 = x∗1 , x2 = x∗2 on M and consequently, 0 = x˙ 2 (t ) = bx∗1 − r2 x∗2 − b1 x∗2 2 − a1 x∗2 y(t )/(1 + mx∗2 ), which yields y(t ) = y∗ . Hence, the only invariant set in M is M = {(x∗1 , x∗2 , y∗ )}. Using LaSalle’s invariant principle, the global asymptotic stability of E ∗ follows. This completes the proof. 5. Numerical simulations In this section, we give two examples to illustrate the main results in Theorems 2.1 and 4.3. Example 1. In (1.3), let a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 16, r = 1/8, r1 = r2 = 1/8, m = 1/10. It is easy to show that (a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) = 12 621/1024 > 0. Hence, system (1.3) has a unique coexistence equilibrium E ∗ (160/597, 10/199, 105 175/118 803). By calculation, we have p0 − q0 ≈ −0.3963 < 0, (p2 + q2 )(p1 + q1 ) − (p0 + q0 ) ≈ 9.7196, τ0 ≈ 3.1489. By Theorem 2.1, E ∗ is locally asymptotically stable if 0 < τ < τ0 and is unstable if τ > τ0 , and system (1.3) undergoes a Hopf bifurcation at E ∗ when τ = τ0 . A numerical simulation illustrates this fact (see Fig. 1). Example 2. In (1.3), let a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 50, r = 1/8, r1 = 1/2, r2 = 1/8, m = 1/10. A direct calculation shows that (a2 − mr )[ab − r2 (b + r1 )] − b1 r (b + r1 ) ≈ 5.0836 > 0. By Theorem 3.1, system (1.3) is permanent. It is easy to show that system (1.3) has a unique coexistence equilibrium E ∗ (0.2010, 0.0503, 0.2739). From the proof of Lemma 3.2, we see that lim supt →+∞ y(t ) ≤ M2 := a2 a2 /(4a1 b1 r0 ) ≈ 0.18. Hence, for ε > 0 sufficiently small, there is a constant T > 0 such that if t > T , y(t ) < M2 + ε . We derive from the first and the second equations of system (1.3) that,
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
2161
Fig. 1. The temporal solution found by numerical integration of system (1.3) with a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 16, r = 1/8, r1 = r2 = 1/8, m = 1/10; (a) τ = 2 and (b) τ = 3.5, respectively; (φ1 , φ2 , ψ) ≡ (0.1, 0.1, 0.1).
Fig. 2. The temporal solution found by numerical integration of system (1.3) with a = 6, a1 = 5, a2 = 5/2, b = 1, b1 = 50, r = 1/8, r1 = 1/2, r2 = 1/8, m = 1/10, τ = 3, (φ1 , φ2 , ψ) ≡ (0.01, 0.01, 0.01).
for t > T , x˙ 1 (t ) = ax2 (t ) − r1 x1 (t ) − bx1 (t ), x˙ 2 (t ) ≥ bx1 (t ) − r2 x2 (t ) − b1 x22 (t ) − a1 x2 (t )(M2 + ε),
(5.1)
which yields lim inf x2 (t ) ≥ x2 := t →+∞
ab − r2 (b + r1 ) − a1 (b + r1 )M2 b1 ( b + r 1 )
≈ 0.0595.
By calculation, we have (ab − r2 (b + r1 ))/(2b1 (b + r1 )) ≈ 0.0388 < x2 . By Theorem 4.3, the equilibrium E ∗ is globally asymptotically stable. A numerical simulation illustrates our result (see Fig. 2). 6. Discussion In this paper, the global dynamics of a stage-structured predator–prey model with Holling type-II functional response and time delay due to the gestation of the predator was investigated using Lyapunov functionals and LaSalle’s invariance principle. It has been shown that, under a hypothesis guaranteeing the uniform persistence of the system, a feasible lower bound condition on the density of the mature prey population ensures the global asymptotic stability of the coexistence equilibrium. That is, if the mature prey is always abundant enough, then the system tends to the unique coexistence equilibrium. On the other hand, it was shown that under some conditions, the time delay due to the gestation of the predator may destabilize the coexistence equilibrium of the system. References [1] C.S. Holling, The functional response of predator to prey density and its role in mimicry and population regulation, Mem. Ent. Sec. Can. 45 (1965) 1–60. [2] L. Chen, Z. Jing, The existence and uniqueness of limit cycles in differential equations modelling the predator–prey interaction, Chin. Sci. Bull. 29 (9) (1984) 521–523.
2162
R. Xu / Nonlinear Analysis: Real World Applications 12 (2011) 2151–2162
[3] X. Liu, L. Chen, Complex dynamics of Holling type II Lotka–Volterra predator–prey system with impulsive perturbations on the predator, Chaos, Solitons Fractals 16 (2003) 311–320. [4] R. Xu, M.A.J. Chaplain, F.A. Davidson, Periodic solutions for a predator–prey model with Holling-type functional response and time delays, Appl. Math. Comput. 161 (2005) 637–654. [5] Wonlyul Ko, Kimun Ryu, Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a prey refuge, J. Differential Equations 231 (2006) 534–550. [6] G. Jiang, Q. Lu, L. Qian, Complex dynamics of a Holling type II prey-predator system with state feedback control, Chaos, Solitons Fractals 31 (2007) 448–461. [7] J.M. Cushing, Integro-differential equations and delay models in population dynamics, in: Lecture Notes in Biomathematics, vol. 20, Springer-Verlag, Berlin, Heidelberg, New York, 1977. [8] K. Gopalsamy, Stability and Oscillations in Delay Differential Equations of Population Dynamics, Kluwer Academic, Dordrecht, Norwell, MA, 1992. [9] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, 1993. [10] M.S. Bartlett, On theoretical models for competitive and predatory biological systems, Biometrika 44 (1957) 27–42. [11] E. Beretta, Y. Kuang, Global analyses in some delayed ratio-dependent predator–prey systems, Nonlinear Anal. TMA 32 (1998) 381–408. [12] P.J. Wangersky, W.J. Cunningham, Time lag in prey-predator population models, Ecology 38 (1957) 136–139. [13] W.G. Aiello, H.I. Freedman, A time delay model of single species growth with stage structure, Math. Biosci. 101 (1990) 139–156. [14] X. Song, L. Chen, Optimal harvesting and stability for a two-species competitive system with stage structure, Math. Biosci. 170 (2001) 173–186. [15] W. Wang, L. Chen, A predator–prey system with stage structure for predator, Comput. Math. Appl. 33 (1997) 83–91. [16] W. Wang, Global dynamics of a population model with stage structure for predator, in: L. Chen, et al. (Eds.), Advanced Topics in Biomathematics, Proceeding of the international conference on mathematical biology, World Scientific Publishing Co. Pte. Ltd, 1997, pp. 253–257. [17] Y. Xiao, L. Chen, Global stability of a predator–prey system with stage structure for the predator, Acta Math. Sin. Engl. Ser. 20 (2004) 63–70. [18] R. Xu, M.A.J. Chaplain, F. Davidson, Persistence and global stability of a ratio-dependent predator–prey model with stage structure, Appl. Math. Comput. 158 (2004) 729–744. [19] X. Zhang, L. Chen, Avidan U. Neumann, The stage-structured predator–prey model and optimal havesting policy, Math. Biosci. 168 (2000) 201–210. [20] J.K. Hale, Theory of Functional Differential Equations, Springer, New York, 1976. [21] J. Hale, P. Waltman, Persistence in infinite-dimensional systems, SIAM J. Math. Anal. 20 (1989) 388–395.