Global Stability for a Delayed Predator-Prey System with Stage Structure for the Predator

Abstract

A delayed predator-prey system with stage structure for the predator is investigated. By analyzing the corresponding characteristic equations, the local stability of equilibria of the system is discussed. The existence of Hopf bifurcation at the positive equilibrium is established. By using an iteration technique and comparison argument, respectively, sufficient conditions are derived for the global stability of the positive equilibrium and two boundary equilibria of the system. Numerical simulations are carried out to illustrate the theoretical results.

1. Introduction

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, and so forth. are often quite different in these two subpopulations. Hence, it is of ecological importance to investigate the effects of such a subdivision on the interaction of species. In [1], Chen et al. introduced the following stage-structured single-species population model:

𝑁̇_𝑖(𝑡)=𝐵(𝑡)−𝐷_𝑖𝑁̇(𝑡)−𝑊(𝑡),_𝑚(𝑡)=𝛼𝑊(𝑡)−𝐷_𝑚(𝑡), (1.1) where 𝑁_𝑖(𝑡) and 𝑁_𝑚(𝑡) denote the immature and mature population densities at time 𝑡, respectively; 𝐵(𝑡) is the birth rate of the immature population at time 𝑡; 𝐷_𝑖(𝑡) and 𝐷_𝑚(𝑡) are the death rates of the immature and mature at time 𝑡; 𝑊(𝑡) represents the transformation rate of the immature into the mature; 𝛼 is the probability of the successful transformation of the immature into the mature. If the birth rate of model (1.1) obeys the Malthus rule, that is, 𝐵(𝑡)=𝑎𝑁_𝑚, the death rates of the immature and mature populations are logistic, that is,

𝐷_𝑖(𝑡)=𝑟_𝑖𝑁_𝑖(𝑡)+𝑏_𝑖𝑁²_𝑖(𝑡),𝐷_𝑚(𝑡)=𝑟_𝑚𝑁_𝑚(𝑡)+𝑏_𝑚𝑁²_𝑚(𝑡), (1.2) and the transformation rate of the immature into mature is proportional to the immature population, that is, 𝑊(𝑡)=𝑏𝑁_𝑖(𝑡), then model (1.1) becomes

𝑁̇_𝑖(𝑡)=𝑎𝑁_𝑚−𝑟_𝑖𝑁_𝑖(𝑡)−𝑏_𝑖𝑁²_𝑖(𝑡)−𝑏𝑁_𝑖𝑁̇(𝑡),_𝑚(𝑡)=𝛼𝑏𝑁_𝑖(𝑡)−𝑟_𝑚𝑁_𝑚(𝑡)−𝑏_𝑚𝑁²_𝑚(𝑡). (1.3) Based on the idea above, many authors studied different kinds of stage-structured models, and a significant body of work has been carried out (see, for example, [2–8]).

In [3], Gao et al. considered the following predator-prey model with stage structure:

𝑥̇(𝑡)=𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂𝑦₁(𝑡−𝜏₁𝑦̇)),₁(𝑡)=𝑦₁(𝑡)[−𝑟₁+𝑎₂₁𝑥(𝑡−𝜏₂)−𝑎₂₂𝑦₁(𝑡)]+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝛼𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡), (1.4) where 𝑥(𝑡) represents the density of the prey at time 𝑡; 𝑦₁(𝑡) and 𝑦₂(𝑡) represent the densities of the mature and the immature predator at time 𝑡, respectively. The parameters 𝑟,𝑟₁,𝑟₂,𝑎₁₁,𝑎₁₂,𝑎₂₁,𝑎₂₂,𝑏,𝜃 are positive constants in which 𝑟 is the intrinsic growth of the prey, 𝑟₁ is the death rate of the mature predator population, 𝑟₂ is the death rate of the immature predator population, 𝑎₁₁ is the intra-specific competition rate of the prey population, 𝑎₁₂ is the capturing rate of the predator population, 𝑎₂₁/𝑎₁₂ is the conversion rate of nutrients into the reproduction of the predator, 𝑎₂₂ is the intra-specific competition rate of the mature predator, 𝑏 is the birth rate of the immature predator, 𝜃 is the transformation rate from the immature predator individuals to mature predator individuals. The predation decreases the average growth rate of prey linearly with a certain time delay 𝜏₁≥0, this assumption corresponds to the fact that predators cannot hunt prey when the predators are infant; predators have to mature for a duration of 𝜏₁ units of time before they are capable of decreasing the average growth rate of the prey species; 𝜏₂≥0 is the time delay due to gestation, the delay in time for prey biomass to increase predator number. In [3], Gao et al. studied the global stability of the positive equilibrium and boundary equilibria of model (1.4) by constructing Liapunov functionals and comparison argument, respectively.

We note that most of the predator-prey models with time delays studied in the literature are all of the Kolmogorov-type. In [9], Wangersky and Cunningham proposed delayed predator-prey models that are not of the Kolmogorov-type. They considered the following delayed system:

𝑥̇(𝑡)=𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂̇𝑦(𝑡)),𝑦(𝑡)=𝑎₂₁𝑥(𝑡−𝜏)𝑦(𝑡−𝜏)−𝑟₁𝑦(𝑡), (1.5) where the delay 𝜏 is a constant based on the assumption that the change rate of predators depends on the number of both the prey and the predators present at some previous time.

Motivated by the work of Gao et al. [3] and Wangersky and Cunningham [9], in the present paper, we consider the following predator-prey model with stage structure and time delay:

𝑥̇(𝑡)=𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂𝑦₁𝑦̇(𝑡)),₁(𝑡)=𝑎₂₁𝑥(𝑡−𝜏)𝑦₁(𝑡−𝜏)−𝑟₁𝑦₁(𝑡)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (1.6) The meanings of the positive parameters 𝑟,𝑟₁,𝑟₂,𝑎₁₁,𝑎₁₂,𝑎₂₁,𝑎₂₂,𝑏,𝜃 are the same as those in system (1.4). The meaning of time delay 𝜏≥0 is the same as in system (1.5).

The initial conditions for system (1.6) take the form

𝑥(𝜃)=𝜙(𝜃),𝑦₁(𝜃)=𝜓₁(𝜃),𝑦₂(𝜃)=𝜓₂(𝜃),𝜙(𝜃)≥0,𝜓₁(𝜃)≥0,𝜓₂(𝜃)≥0,𝜃∈[−𝜏,0],𝜙(0)>0,𝜓₁(0)>0,𝜓₂(0)>0, (1.7) where (𝜙₁(𝜃),𝜓₁(𝜃),𝜓₂(𝜃))∈𝐶([−𝜏,0],𝑅³₊₀), the Banach space of continuous functions mapping the interval [−𝜏,0] into 𝑅³₊₀, where 𝑅³₊₀={(𝑥₁,𝑥₂,𝑥₃)∶𝑥_𝑖≥0,𝑖=1,2,3}.

This paper is organized as follows. In the next section, we introduce some notations and state several lemmas which will be essential to our proofs. In Section 3, we discuss the local stability of a positive equilibrium and boundary equilibria of system (1.6). The existence of Hopf bifurcation is studied. In Section 4, by means of an iterative technique and comparison argument, sufficient conditions are derived for the global stability of the positive equilibrium and boundary equilibria of system (1.6). Some numerical examples are given to illustrate the results above. A brief discussion is given in Section 6 to end this work.

2. Preliminaries

In this section, we introduce some notations and state several results which will be useful in next section. Let ℝ^𝑛₊ be the cone of nonnegative vectors in ℝ^𝑛. If 𝑥,𝑦∈ℝ^𝑛, we write 𝑥≤𝑦(𝑥<𝑦) if 𝑥_𝑖≤𝑦_𝑖(𝑥_𝑖<𝑦_𝑖) for 1≤𝑖≤𝑛. Let 𝑒₁,𝑒₂,…,𝑒_𝑛 denote the standard basis in ℝ^𝑛. Suppose 𝑟≥0 and let 𝐶=𝐶([−𝑟,0],ℝ^𝑛) be the Banach space of continuous functions mapping the interval [−𝑟,0] into ℝ^𝑛 with supremum norm. If 𝜙,𝜓∈𝐶, we write 𝜙≤𝜓(𝜙<𝜓) when the indicated inequality holds at each point of [−𝑟,0]. Let 𝐶⁺={𝜙∈𝐶∶𝜙≥0} and let denote the inclusion ℝ^𝑛→𝐶([−𝑟,0],ℝ^𝑛) by 𝑥→𝑥,𝑥(𝜃)=𝑥,𝜃∈[−𝑟,0]. Denote the space of functions of bounded variation on [−𝑟,0] by BV[−𝑟,0]. If 𝑡₀∈ℝ^𝑛,𝐴≥0, and 𝑥∈𝐶([−𝑡₀−𝑟,𝑡₀+𝐴],ℝ^𝑛), then for any 𝑡∈[𝑡₀,𝑡₀+𝐴], we let 𝑥_𝑡∈𝐶 be defined by 𝑥_𝑡(𝜃)=𝑥(𝑡+𝜃),−𝑟≤𝜃≤0.

We now consider

𝑥̇(𝑡)=𝑓(𝑡,𝑥_𝑡). (2.1) We assume throughout this section that 𝑓∶ℝ×𝐶→ℝ^𝑛 is continuous; 𝑓(𝑡,𝜙) is continuously differentiable in 𝜙; 𝑓(𝑡+𝑇,𝜙)=𝑓(𝑡,𝜙) for all (𝑡,𝜙)∈ℝ×𝐶⁺, and some 𝑇>0. Then by [10], there exists a unique solution of (2.1) through (𝑡₀,𝜙) for 𝑡₀∈ℝ,𝜙∈𝐶⁺. This solution will be denoted by 𝑥(𝑡,𝑡₀,𝜙) if we consider the solution in ℝ^𝑛, or by 𝑥_𝑡(𝑡₀,𝜙) if we work in the space C. Again by [10], 𝑥(𝑡,𝑡₀,𝜙)(𝑥_𝑡(𝑡₀,𝜙)) is continuously differentiable in 𝜙. In the following, the notation 𝑥_𝑡0=𝜙 will be used as the condition of the initial data of (2.1), by which we mean that we consider the solution 𝑥(𝑡) of (2.1) which satisfies 𝑥(𝑡₀+𝜃)=𝜙(𝜃),𝜃∈[−𝑟,0].

To proceed further, we need the following results from [11, 12]. Let 𝑟=(𝑟₁,𝑟₂,…,𝑟_𝑛)∈ℝ^𝑛,|𝑟|=max_𝑖{𝑟_𝑖}, and define

𝐶_𝑟=^𝑛_𝑖=1𝐶([−𝑟_𝑖,0],ℝ). (2.2) We write 𝜙=(𝜙₁,𝜙₂,…,𝜙_𝑛) for a generic point of 𝐶_𝑟. Let 𝐶⁺_𝑟={𝜙∈𝐶_𝑟∶𝜙≥0}. Due to the ecological applications, we choose 𝐶⁺_𝑟 as the state space of (2.1) in the following discussions.

Fix 𝜙₀∈𝐶⁺_𝑟 arbitrarily. Then we set 𝐿(𝑡,⋅)=𝐷_𝜙0𝑓(𝑡,𝜙₀), and 𝐷_𝜙0𝑓(𝑡,𝜙₀) denotes the Frechet derivation of 𝑓 with respect to 𝜙₀. It is convenient to have the standard representation of 𝐿=(𝐿₁,𝐿₂,…,𝐿_𝑛) as

𝐿_𝑖(𝑡,𝜙)=^𝑛_𝑗=1⁰_−𝑟𝑗𝜙_𝑗(𝜃)𝑑_𝜃𝜂_𝑖𝑗(𝜃,𝑡)(1≤𝑖≤𝑛), (2.3) in which 𝜂_𝑖𝑗∶ℝ×ℝ→ℝ satisfies

𝜂_𝑖𝑗(𝜃,𝑡)=𝜂_𝑖𝑗(0,𝑡),𝜃≤−𝑟_𝑗,𝜂_𝑖𝑗𝜂(𝜃,𝑡)=0,𝜃≥0,_𝑖𝑗(⋅,𝑡)∈BV[−𝑟_𝑗,0], (2.4) where 𝜂_𝑖𝑗(⋅,𝑡) is continuous from the left in (−𝑟_𝑗,0).

We make the following assumptions for (2.1).

(h0)If 𝜙,𝜓∈𝐶⁺,𝜙≤𝜑, and 𝜙_𝑖(0)=𝜑_𝑖(0) for some 𝑖, then 𝑓_𝑖(𝑡,𝜙)≤𝑓_𝑖(𝑡,𝜓).(h1)For all 𝜙∈𝐶⁺_𝑟 with 𝜙_𝑖(0)=0,𝐿_𝑖(𝑡,𝜙)≥0 for 𝑡∈ℝ.(h2)The matrix 𝐴(𝑡) defined by 𝐴(𝑡)=col(𝐿(𝑡,𝑒₁̂),𝐿(𝑡,𝑒₂̂),…,𝐿(𝑡,𝑒_𝑛̂))=(𝜂_𝑖𝑗(0,𝑡)), (2.5) is irreducible for each 𝑡∈ℝ. (h3)For each 𝑗, for which 𝑟_𝑗>0, there exist 𝑖 such that for all 𝑡∈ℝ and for positive constant 𝜀 sufficiently small, 𝜂_𝑖𝑗(−𝑟_𝑗+𝜀,𝑡)>0.(h4)If 𝜙=0, then 𝑥(𝑡,𝑡₀,𝜙)≡0 for all 𝑡≥𝑡₀.

The following result was established by Wang et al. [12].

Lemma 2.1. Let (h1)–(h4) hold. Then the hypothesis (h0) is valid and the following.(i)If 𝜙 and 𝜓 are distinct elements of 𝐶⁺_𝑟 with 𝜙≤𝜓 and [𝑡₀,𝑡₀+𝜎) with 𝑛|𝑟|<𝜎≤∞ is the intersection of the maximal intervals of existence of 𝑥(𝑡,𝑡₀,𝜙) and 𝑥(𝑡,𝑡₀,𝜓), then 0≤𝑥(𝑡,𝑡₀,𝜙)≤𝑥(𝑡,𝑡₀,𝜓)for𝑡₀≤𝑡<𝑡₀+𝜎,0≤𝑥(𝑡,𝑡₀,𝜙)<𝑥(𝑡,𝑡₀,𝜓)for𝑡₀+𝑛|𝑟|≤𝑡<𝑡₀+𝜎. (2.6)(ii)If 𝜙∈𝐶⁺_𝑟,𝜙≠0,𝑡₀∈ℝ and 𝑥(𝑡,𝑡₀,𝜙) is defined on [𝑡₀,𝑡₀+𝜎) with 𝜎>𝑛|𝑟|, then 0<𝑥(𝑡,𝑡₀,𝜙)for𝑡₀+𝑛|𝑟|≤𝑡<𝑡₀+𝜎. (2.7)

Lemma 2.1 shows that if (h1)–(h4) hold, then the positivity of solutions of (2.1) follows.

The following definitions and results are useful in proving our lemma.

Definition 2.2. System (2.1) is cooperative if 𝜕𝑓_𝑖/𝜕𝑥_𝑗≥0 whenever 𝑖≠𝑗.

Definition 2.3. A square matrix 𝐴 is said to be a reducible matrix if and only if for some permutation matrix 𝑃 the matrix 𝑃^𝑇𝐴𝑃 is block upper triangular. If a square matrix is not reducible, it is said to be an irreducible matrix. System (2.1) is called irreducible if the Jacobian matrix ((𝜕𝑓_𝑖/𝜕𝑥_𝑗)) is irreducible.

Lemma 2.4 (Smith [11]). If (2.1) is cooperative and irreducible in 𝐷, where 𝐷 is an open subset of 𝐶, and the solution with positive initial data is bounded, then the trajectory of (2.1) tends to some single equilibrium.

We now consider the following delay differential system:

𝑦̇₁(𝑡)=−𝑎𝑦₁(𝑡)+𝑏𝑦₁(𝑡−𝜏)−𝑐𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑑𝑦₁(𝑡)−(𝜃+𝑒)𝑦₂(𝑡), (2.8) with initial conditions

𝑦_𝑖(𝑠)=𝜓_𝑖(𝑠)≥0,𝑠∈[−𝜏,0),𝜓_𝑖∈𝐶([−𝜏,0),ℝ₊),𝜓_𝑖(0)>0,(𝑖=1,2). (2.9)

System (2.8) always has a trivial equilibrium 𝐸₀(0,0). If 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), then system (2.8) has a unique positive equilibrium 𝐸₊(𝑦^∗₁,𝑦^∗₂), where

𝑦^∗₁=(−𝑎+𝑏)(𝜃+𝑒)+𝜃𝑑𝑐(𝜃+𝑒),𝑦^∗₂=𝑑(−𝑎+𝑏)(𝜃+𝑒)+𝜃𝑑²𝑐(𝜃+𝑒)². (2.10)

The characteristic equation of system (1.6) at the equilibrium 𝐸₀ is of the form

𝜆²+(𝑎+𝑒+𝜃)𝜆+𝑎(𝜃+𝑒)−𝑑𝜃−(𝑏𝜆+𝑏(𝜃+𝑒))𝑒^−𝜆𝜏=0. (2.11) Let

𝑓₁(𝜆)=𝜆²+(𝑎+𝑒+𝜃)𝜆+𝑎(𝜃+𝑒)−𝑑𝜃−(𝑏𝜆+𝑏(𝜃+𝑒))𝑒^−𝜆𝜏. (2.12) If 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), then it is easy to see that, for 𝜆 real,

𝑓₁(0)=𝑎(𝜃+𝑒)−𝑑𝜃−𝑏(𝜃+𝑒)<0,lim_𝜆→+∞𝑓₁(𝜆)=+∞. (2.13) Hence, 𝑓₁(𝜆)=0 has a positive real root. Therefore, the equilibrium 𝐸₀ is unstable. If 𝑎(𝜃+𝑒)>𝑑𝜃+𝑏(𝜃+𝑒), when 𝜏=0, it is easy to see that the equilibrium 𝐸₀ is stable. Therefore, if 𝑎(𝜃+𝑒)−𝑑𝜃>𝑏(𝜃+𝑒), by Kuang and So [13, Lemma ?B], we see that the equilibrium 𝐸₀ is locally stable for all 𝜏>0.

If 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), the characteristic equation of the positive equilibrium 𝐸₊ takes the form

𝜆²+(𝑎+𝑒+𝜃+2𝑐𝑦^∗₁)𝜆+(𝑎+2𝑐𝑦^∗₁)(𝜃+𝑒)−𝑑𝜃−(𝑏𝜆+𝑏(𝜃+𝑒))𝑒^−𝜆𝜏=0. (2.14) When 𝜏=0, it is easy to see that the equilibrium 𝐸₊ is stable. Therefore, if 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), by Kuang and So [13, Lemma ?B], we see that the equilibrium 𝐸₊ is locally stable for all 𝜏>0.

Lemma 2.5. For system (2.8), one has the following. (i)If 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), then the positive equilibrium 𝐸₊ of system (2.8) is globally stable.(ii)If 𝑎(𝜃+𝑒)>𝑑𝜃+𝑏(𝜃+𝑒), the equilibrium 𝐸₀ of system (2.8) is globally stable.

Proof. We represent the right-hand side of (2.8) by 𝑓(𝑡,𝑥_𝑡)=(𝑓₁(𝑡,𝑥_𝑡),𝑓₂(𝑡,𝑥_𝑡)) and set 𝐿(𝑡,⋅)=𝐷_𝜙𝑓(𝑡,𝜙). By direct calculation we have 𝐿₁(𝑡,ℎ)=−𝑎ℎ₁(0)+𝑏ℎ₁(−𝜏)−2𝑐𝜙₁(0)ℎ₁(0)+𝜃ℎ₂𝐿(0),₂(𝑡,ℎ)=𝑑ℎ₁(0)−(𝜃+𝑒)ℎ₂(0). (2.15) We now claim that the hypotheses (h1)–(h4) hold for system (2.8). It is easily seen that (h1) and (h4) hold for system (2.8). We need only to verify that (h2) and (h3) hold.
The matrix 𝐴(𝑡) takes the form
−𝑎+𝑏−2𝑐𝜙₁(0)𝜃𝑑−(𝜃+𝑒). (2.16) Clearly, the matrix 𝐴(𝑡) is irreducible for each 𝑡∈ℝ.
From the definition of 𝐴(𝑡) and 𝜂_𝑖𝑗, it is readily seen that 𝜂₁₂(𝜃,𝑡)=𝜂₁₂(0,𝑡)=𝜃,𝜂₂₁(𝜃,𝑡)=𝜂₂₁(0,𝑡)=𝑑, for 𝜃≥0; and 𝜂_𝑖𝑗(𝜃,𝑡)=0,𝑖≠𝑗 for 𝜃≤−𝜏; and 𝜂_𝑖𝑗(⋅,𝑡)∈BV[−𝜏,0], where 𝜂_𝑖𝑗 is a positive Borel measure on [−𝜏,0]. Therefore, 𝜂_𝑖𝑗(⋅,𝑡)>0. Thus, for each 𝑗, there is 𝑖≠𝑗 such that 𝜂_𝑖𝑗(−𝑟_𝑗+𝜀,𝑡)=𝜂_𝑖𝑗(−𝜏+𝜀,𝑡)>0 for all 𝑡∈ℝ and 𝜀>0 sufficiently small, 𝑖=1,2. Hence, (h3) holds.
Thus, the conditions of Lemma 2.1 are satisfied. Therefore, the positivity of solutions to system (2.8) follows. It is easy to see that system (2.8) is cooperative. By Lemma 2.4 we see that any solution starting from 𝐷=𝐶⁺_𝜏 converges to some single equilibrium. However, system (2.8) has only two equilibria: 𝐸₀ and 𝐸₊. Note that if 𝑎(𝜃+𝑒)<𝑑𝜃+𝑏(𝜃+𝑒), the equilibrium 𝐸₊ is locally stable. Hence, any solution starting from 𝐷 converges to 𝐸⁺. Using a similar argument one can show the global stability of the equilibrium 𝐸₀ when 𝑎(𝜃+𝑒)>𝑑𝜃+𝑏(𝜃+𝑒). This completes the proof.

By a similar argument one can show that all solutions of system (1.6) with initial conditions (1.7) are defined on [0,+∞) and remain positive for all 𝑡≥0.

3. Local Stability

In this section, we discuss the local stability of each equilibria and the existence of Hopf bifurcation of system (1.6).

It is easy to show that system (1.6) always has two equilibria 𝐸₀(0,0,0) and 𝐸₁(𝑟/𝑎₁₁,0,0). If the following holds:

(H1)𝜃𝑏>𝑟₁(𝜃+𝑟₂),

then system (1.6) has another boundary equilibrium 𝐸₂(0,𝑦₁,𝑦₂), where 𝑦₁=1𝑎₂₂(𝜃𝑏𝜃+𝑟₂−𝑟₁),𝑦₂=𝑏𝑎₂₂(𝜃+𝑟₂)(𝜃𝑏𝜃+𝑟₂−𝑟₁). (3.1) Further, if the following holds:

(H2)𝑟(𝜃+𝑟₂)(𝑎₁₂𝑎₂₁+𝑎₁₁𝑎₂₂)>𝑎₁₂[𝜃𝑏𝑎₁₁+(𝑎₂₁𝑟−𝑎₁₁𝑟₁)(𝜃+𝑟₂)]>0,

then system (1.6) has a unique positive equilibrium 𝐸^∗=(𝑥^∗,𝑦^∗₁,𝑦^∗₂), where

𝑥^∗=𝑟𝑎₁₁−𝑎₁₂[𝜃𝑏𝑎₁₁+(𝜃+𝑟₂)(𝑎₂₁𝑟−𝑎₁₁𝑟₁)](𝜃+𝑟₂)𝑎₁₁(𝑎₁₂𝑎₂₁+𝑎₁₁𝑎₂₂),𝑦^∗₁=𝜃𝑏𝑎₁₁+(𝜃+𝑟₂)(𝑎₂₁𝑟−𝑎₁₁𝑟₁)(𝜃+𝑟₂)(𝑎₁₂𝑎₂₁+𝑎₁₁𝑎₂₂),𝑦^∗₂=𝑏[𝜃𝑏𝑎₁₁+(𝜃+𝑟₂)(𝑎₂₁𝑟−𝑎₁₁𝑟₁)](𝜃+𝑟₂)²(𝑎₁₂𝑎₂₁+𝑎₁₁𝑎₂₂). (3.2)

We now study the local stability of each of the nonnegative equilibrium of system (1.6).

Let 𝐸=(𝑥,𝑦₁,𝑦₂) be any arbitrary equilibrium. Then the characteristic equation of system (1.6) at the equilibrium 𝐸 is given by

𝑟−2𝑎₁₁𝑥−𝑎₁₂𝑦₁−𝜆−𝑎₁₂𝑎𝑥0₂₁𝑦₁𝑒^−𝜆𝜏−𝑟₁−2𝑎₂₂𝑦₁+𝑎₂₁𝑥𝑒^−𝜆𝜏−𝜆𝜃0𝑏−𝜃−𝑟₂−𝜆=0. (3.3) The characteristic equation of system (1.6) at the equilibrium 𝐸₀(0,0,0) reduces

(𝜆−𝑟)[(𝑟₁+𝜆)(𝜃+𝑟₂+𝜆)−𝑏𝜃]=0. (3.4) Clearly, 𝜆=𝑟 is a positive real root. Hence, 𝐸₀(0,0,0) is always unstable.

The characteristic equation of system (1.6) at the equilibrium 𝐸₁(𝑟/𝑎₁₁,0,0) reduces

(𝜆+𝑟)[𝜆²+(𝑟₁+𝑟₂+𝜃)𝜆+𝑟₁(𝜃+𝑟₂𝑎)−𝑏𝜃−²¹𝑟𝑎₁₁(𝜆+𝜃+𝑟₂)𝑒^−𝜆𝜏]=0. (3.5) Clearly, 𝜆=−𝑟 is a negative real root of (3.5). All other roots are give by the roots of equation

𝜆²+(𝑟₁+𝑟₂+𝜃)𝜆+𝑟₁(𝜃+𝑟₂𝑎)−𝑏𝜃−²¹𝑟𝑎₁₁(𝜆+𝜃+𝑟₂)𝑒^−𝜆𝜏=0. (3.6) Let

𝑓(𝜆)=𝜆²+(𝑟₁+𝑟₂+𝜃)𝜆+𝑟₁(𝜃+𝑟₂𝑎)−𝑏𝜃−²¹𝑟𝑎₁₁(𝜆+𝜃+𝑟₂)𝑒^−𝜆𝜏. (3.7) If (𝜃+𝑟₂)(𝑟₁−𝑎₂₁𝑟/𝑎₁₁)<𝑏𝜃, then it is easy to see that for 𝜆 real,

𝑓(0)=(𝜃+𝑟₂𝑟)(¹−𝑎₂₁𝑟𝑎₁₁)−𝑏𝜃<0,lim_𝜆→+∞𝑓(𝜆)=+∞. (3.8) Hence, 𝑓(𝜆)=0 has a positive real root. Therefore, the equilibrium 𝐸₁(𝑟/𝑎₁₁,0,0) is unstable. If (𝜃+𝑟₂)(𝑟₁−𝑎₂₁𝑟/𝑎₁₁)>𝑏𝜃, when 𝜏=0, it is easy to see that the equilibrium 𝐸₁(𝑟/𝑎₁₁,0,0) is stable. Therefore, if (𝜃+𝑟₂)(𝑟₁−𝑎₂₁𝑟/𝑎₁₁)>𝑏𝜃, by Kuang and So [13, Lemma ?B], we see that the equilibrium 𝐸₁(𝑟/𝑎₁₁,0,0) is locally stable for all 𝜏>0.

The characteristic equation of system (1.6) at the equilibria 𝐸₂(0,𝑦₁,𝑦₂) reduces to

(𝜆−𝑟+𝑎₁₂𝑦₁)[𝜆²+(𝑟₁+2𝑎₂₂𝑦₁+𝜃+𝑟₂)𝜆+(𝑟₁+2𝑎₂₂𝑦₁)(𝜃+𝑟₂)−𝑏𝜃]=0. (3.9) Clearly, if 𝑟>𝑎₁₂𝑦₁, 𝜆=𝑟−𝑎₁₂𝑦₁ is a positive real root. Hence, if 𝑟>𝑎₁₂𝑦₁, then 𝐸₂(0,𝑦₁,𝑦₂) is unstable. Noting that

(𝑟₁+2𝑎₂₂𝑦₁)(𝜃+𝑟₂)−𝑏𝜃=𝑏𝜃−𝑟₁(𝜃+𝑟₂)>0, (3.10) if 𝑟<𝑎₁₂𝑦₁, then (3.9) only has negative real root, and 𝐸₂(0,𝑦₁,𝑦₂) is stable.

The characteristic equation of system (1.6) at the positive equilibria 𝐸^∗(𝑥^∗,𝑦^∗₁,𝑦^∗₂) is

𝜆³+𝑝₂𝜆²+𝑝₁𝜆+𝑝₀+(𝑞₂𝜆²+𝑞₁𝜆+𝑞₀)𝑒^−𝜆𝜏=0, (3.11) where

𝑝₀=𝑎₁₁𝑥^∗(𝜃+𝑟₂)(𝑎₂₁𝑥^∗+𝑎₂₂𝑦^∗₁𝑝),₁=𝑎₁₁𝑥^∗(𝑟₁+2𝑎₂₂𝑦^∗₁)+(𝜃+𝑟₂)(𝑎₂₁𝑥^∗+𝑎₂₂𝑦^∗₁+𝑎₁₁𝑥^∗𝑝),₂=𝑎₁₁𝑥^∗+2𝑎₂₂𝑦^∗₁+𝑟₁+𝑟₂𝑞+𝜃,₀=𝑎₂₁(𝜃+𝑟₂)𝑥^∗(𝑎₁₂𝑦^∗₁−𝑎₁₁𝑥^∗𝑞),₁=𝑎₂₁𝑥^∗(𝑎₁₂𝑦^∗₁−𝑎₁₁𝑥^∗−𝜃−𝑟₂𝑞),₂=−𝑎₂₁𝑥^∗. (3.12) It is easy to show that

𝑝₀+𝑞₀=(𝑎₁₁𝑎₂₂+𝑎₂₁𝑎₁₂)(𝜃+𝑟₂)𝑥^∗𝑦^∗₁𝑝>0,₁+𝑞₁=(𝑎₁₁𝑎₂₂+𝑎₁₂𝑎₂₁)𝑥^∗𝑦^∗₁+(𝜃+𝑟₂)(𝑎₂₂𝑦^∗₁+𝑎₁₁𝑥^∗)+𝑏𝜃𝑎₁₁𝑥^∗(𝜃+𝑟₂)𝑝>0,₂+𝑞₂=𝑎₁₁𝑥^∗+𝑎₂₂𝑦^∗₁+𝑟₂+𝜃+𝑏𝜃(𝜃+𝑟₂)>0. (3.13) It is easy to see that (𝑝₂+𝑞₂)(𝑝₁+𝑞₁)>𝑝₀+𝑞₀. Hence, by the Routh-Hurwitz theorem, when 𝜏=0, the positive equilibrium 𝐸^∗ of system (1.6) is locally asymptotically stable.

If i𝜔(𝜔>0) is a root of (3.11), separating the real and imaginary parts, we obtain

𝑝₀−𝑝₂𝜔²=𝑞₂𝜔²cos𝜔𝜏−𝑞₁𝜔sin𝜔𝜏−𝑞₀cos𝜔𝜏,−𝜔³+𝑝₁𝜔=−𝑞₂𝜔²sin𝜔𝜏−𝑞₁𝜔cos𝜔𝜏+𝑞₀sin𝜔𝜏. (3.14) Squaring and adding the two equations of (3.14), it follows that

𝜔⁶+𝑄₁𝜔⁴+𝑄₂𝜔²+𝑄₃=0, (3.15) where

𝑄₁=𝑝²₂−2𝑝₁−𝑞²₂,𝑄₂=𝑝²₁+2𝑞₂𝑞₀−2𝑝₂𝑝₀−𝑞²₁,𝑄₃=𝑝²₀−𝑞²₀. (3.16) It is easy to show that

𝑄₁=𝑝²₂−𝑞²₂−2𝑝₁=𝑎²₁₁(𝑥^∗)²+(𝑟₂+𝜃)²+(𝑟₁+2𝑎₂₂𝑦^∗₁+𝑎₂₁𝑥^∗)(𝑎₂₂𝑦^∗₁+𝑏𝜃𝜃+𝑟₂𝑄)+2𝑏𝜃>0,₂=𝑝²₁+2𝑞₂𝑞₀−2𝑝₂𝑝₀−𝑞²₁=[(𝑎₁₁𝑎₂₂+𝑎₂₁𝑎₁₂)𝑥^∗𝑦^∗₁+𝑎₁₁𝑏𝜃𝑥^∗(𝑟₂+𝜃)][2𝑎₁₁𝑎₂₁(𝑥^∗)²−(𝑎₁₂𝑎₂₁−𝑎₁₁𝑎₂₂)𝑥^∗𝑦^∗₁+𝑎₁₁𝑏𝜃𝑥^∗(𝑟₂]+𝜃)+(𝜃+𝑟₂)²[𝑎²₂₂(𝑦^∗₁)²+𝑎²₁₁(𝑥^∗)²]+2𝑎₂₁𝑎₂₂(𝜃+𝑟₂)²𝑥^∗+2𝑏𝜃𝑎²₁₁(𝑥^∗)². (3.17)

Noting that if 𝑝₀>𝑞₀, 2𝑎₁₁𝑎₂₁𝑥^∗>(𝑎₁₂𝑎₂₁−𝑎₁₁𝑎₂₂)𝑦^∗₁, hence 𝑄₂>0,𝑄₃>0. The positive equilibrium 𝐸^∗ of system (1.6) is locally asymptotically stable for all 𝜏>0. If 𝑝₀<𝑞₀, we know that 𝑄₃<0, (3.15) has a unique positive root 𝜔₀. Define

𝜏_𝑗=1𝜔₀arccos(𝑞₁−𝑞₂𝑝₂)𝜔⁴₀+(𝑝₀𝑞₂+𝑝₂𝑞₀−𝑝₁𝑞₁)𝜔²₀−𝑞₀𝑝₀𝑞²₂𝜔⁴₀+(𝑞²₁−2𝑞₀𝑞₂)𝜔²₀+𝑞²₀+2𝑗𝜋𝜔₀,𝑗=0,1,…. (3.18) Then (𝜏_𝑗,𝜔₀) solves (3.11). This means that when 𝜏=𝜏_𝑗, (3.11) has a pair of purely imaginary roots ±i𝜔₀. Noting that the positive equilibrium 𝐸^∗ is locally stable for 𝜏=0, by the general theory on characteristic equations of delay differential equations from [14, Theorem ?4.1], 𝐸^∗ remains stable for 𝜏<𝜏₀.

Let 𝑝₀<𝑞₀ and 𝜏₀ be defined in (3.18). Denoting

𝜆(𝜏)=𝛼(𝜏)+i𝜔(𝜏), (3.19) the root of (3.11) is such that

𝛼(𝜏₀)=0,𝜔(𝜏₀)=𝜔₀. (3.20) In the following we claim that

𝑑(Re𝜆)|𝑑𝜏_𝜏=𝜏0>0. (3.21) This will signify that there exists at least one eigenvalue with positive real part for 𝜏>𝜏₀. Moreover, the conditions for the existence of a Hopf bifurcation [10] are then satisfied yielding a periodic solution. To this end, differentiating equation (3.11) with respect 𝜏, we obtain that

(𝑑𝜆(𝜏))d𝜏⁻¹=−3𝜆²+2𝑝₂𝜆+𝑝₁𝜆(𝜆³+𝑝₂𝜆²+𝑝₁𝜆+𝑝₀)+2𝑞₂𝜆+𝑞₁𝜆(𝑞₂𝜆²+𝑞₁𝜆+𝑞₀)−𝜏𝜆, (3.22) which leads to

sign{Re(𝑑𝜆𝑑𝜏)|_{𝜏=𝜏𝑗}}=sign{Re(𝑑𝜆𝑑𝜏)|⁻¹_{𝜏=𝜏𝑗}}=sign{3𝜔⁴₀+(2𝑝²₂−4𝑝₁−2𝑞²₂)𝜔²₀+𝑝²₁+2𝑞₂𝑞₀−2𝑝₂𝑝₀−𝑞²₁}=sign{3𝜔⁴₀+2𝑄₁𝜔²₀+𝑄₂}. (3.23) If 𝑄₂>0, 3𝜔⁴₀+2𝑄₁𝜔²₀+𝑄₂>0. Therefore, system (1.6) undergoes a Hopf bifurcation.

We therefore obtain the following results.

Theorem 3.1. For system (1.6), let 𝜏₀ be defined as in (3.18), one has the following.(i)The positive equilibrium 𝐸₀ of system (1.6) is always unstable.(ii)If (𝜃+𝑟₂)(𝑟₁−𝑎₂₁𝑟/𝑎₁₁)>𝑏𝜃, the equilibrium 𝐸₁(𝑟/𝑎₁₁,0,0) of system (1.6) is locally stable for all 𝜏>0; and if (𝜃+𝑟₂)(𝑟₁−𝑎₂₁𝑟/𝑎₁₁)<𝑏𝜃, 𝐸₁ is unstable for all 𝜏.(iii)Let (𝐻1) hold. If 𝑟>𝑎₁₂𝑦₁, the equilibrium 𝐸₂(0,𝑦₁,𝑦₂) of system (1.6) is unstable; if 𝑟<𝑎₁₂𝑦₁, 𝐸₂ is stable for all 𝜏≥0.(iv)Let (H2) hold. If 𝑝₀>𝑞₀, then the positive equilibrium 𝐸^∗(𝑥^∗,𝑦^∗₁,𝑦^∗₂) of system (1.6) is locally asymptotically stable for all 𝜏>0. If 𝑝₀<𝑞₀, then 𝐸^∗ is locally stable for 𝜏<𝜏₀; and 𝐸^∗ is unstable for 𝜏>𝜏₀; if 𝑄₂>0 system (1.6) undergoes a Hopf Bifurcation at the positive equilibrium 𝐸^∗ when 𝜏=𝜏₀.

4. Global Stability

In this section, we are concerned with the global stability of the equilibria 𝐸₁,𝐸₂,𝐸^∗ of system (1.6). The strategy of proofs is to use an iteration technique and comparison arguments, respectively.

Theorem 4.1. Let (H2) hold. Then the positive equilibrium 𝐸^∗ of system (1.6) is globally asymptotically stable provided that (H3)𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁.

Proof. Let (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) be any positive solution of system (1.6) with initial conditions (1.7). Let 𝑈=limsup_𝑡→+∞𝑥(𝑡),𝑉=liminf_𝑡→+∞𝑈𝑥(𝑡),_𝑖=limsup_𝑡→+∞𝑦_𝑖(𝑡),𝑉_𝑖=liminf_𝑡→+∞𝑦_𝑖(𝑡)(𝑖=1,2). (4.1) We now claim that 𝑈=𝑉=𝑥^∗,𝑈₁=𝑉₁=𝑦^∗₁,𝑈₂=𝑉₂=𝑦^∗₂. The strategy of the proof is to use an iteration technique.
We derive from the first equation of system (1.6) that
𝑥̇(𝑡)≤𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)). (4.2) A standard comparison argument shows that 𝑈=limsup_𝑡→+∞𝑟𝑥(𝑡)≤𝑎₁₁∶=𝑀^𝑥₁. (4.3) Hence, for 𝜀>0 sufficiently small there exists a 𝑇₁>0 such that if 𝑡>𝑇₁, 𝑥(𝑡)≤𝑀^𝑥₁+𝜀.
We derive from the second and the third equations of system (1.6) that for 𝑡>𝑇₁+𝜏,
𝑦̇₁(𝑡)≤−𝑟₁𝑦₁(𝑡)+𝑎₂₁(𝑀^𝑥₁+𝜀)𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.4) Consider the following auxiliary equations: 𝑢̇₁(𝑡)=−𝑟₁𝑢₁(𝑡)+𝑎₂₁(𝑀^𝑥₁+𝜀)𝑢₁(𝑡−𝜏)−𝑎₂₂𝑢²₁(𝑡)+𝜃𝑢₂𝑢̇(𝑡),₂(𝑡)=𝑏𝑢₁(𝑡)−(𝜃+𝑟₂)𝑢₂(𝑡). (4.5) Since (H2) holds, by Lemma 2.5 it follows from (4.5) that lim_𝑡→+∞𝑢₁(𝑡)=(−𝑟₁+𝑎₂₁(𝑀^𝑥₁+𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),lim_𝑡→+∞𝑢₂(𝑡)=𝑏(−𝑟₁+𝑎₂₁(𝑀^𝑥₁+𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.6) By comparison, we obtain that 𝑈₁=limsup_𝑡→+∞𝑦₁(𝑡)≤(−𝑟₁+𝑎₂₁(𝑀^𝑥₁+𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑈₂=limsup_𝑡→+∞𝑦₂(𝑡)≤𝑏(−𝑟₁+𝑎₂₁(𝑀^𝑥₁+𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.7) Since these inequalities are true for arbitrary 𝜀>0 sufficiently small, it follows that 𝑈₁≤𝑀^𝑦1₁,𝑈₂≤𝑀^𝑦2₁, where 𝑀^𝑦1₁=(−𝑟₁+𝑎₂₁𝑀^𝑥₁)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑀^𝑦2₁=𝑏(−𝑟₁+𝑎₂₁𝑀^𝑥₁)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.8) Hence, for 𝜀>0 sufficiently small, there is a 𝑇₂≥𝑇₁+𝜏 such that if 𝑡>𝑇₂, 𝑦₁(𝑡)≤𝑀^𝑦1₁+𝜀, 𝑦₂(𝑡)≤𝑀^𝑦2₁+𝜀.
For 𝜀>0 sufficiently small, we derive from the first equation of system (1.6) that, for 𝑡>𝑇₂,
𝑥̇(𝑡)≥𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂(𝑀^𝑦1₁+𝜀)). (4.9) By comparison it follows that 𝑉=liminf_𝑡→+∞𝑥(𝑡)≥𝑟−𝑎₁₂(𝑀^𝑦1₁+𝜀)𝑎₁₁. (4.10) Since this is true for arbitrary 𝜀>0 sufficiently small, we conclude that 𝑉≥𝑁^𝑥₁, where 𝑁^𝑥₁=𝑟−𝑎₁₂𝑀^𝑦1₁𝑎₁₁. (4.11) Therefore, for 𝜀>0 sufficiently small, there is a 𝑇₃≥𝑇₂ such that if 𝑡>𝑇₃, 𝑥(𝑡)≥𝑁^𝑥₁−𝜀.
For 𝜀>0 sufficiently small, we derive from the second and the third equations of system (1.6) that, for 𝑡>𝑇₃+𝜏,
𝑦̇₁(𝑡)≥−𝑟₁𝑦₁(𝑡)+𝑎₂₁(𝑁^𝑥₁−𝜀)𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.12) Consider the following auxiliary equations: 𝑢̇₁(𝑡)=−𝑟₁𝑢₁(𝑡)+𝑎₂₁(𝑁^𝑥₁−𝜀)𝑢₁(𝑡−𝜏)−𝑎₂₂𝑢²₁(𝑡)+𝜃𝑢₂𝑢̇(𝑡),₂(𝑡)=𝑏𝑢₁(𝑡)−(𝜃+𝑟₂)𝑢₂(𝑡). (4.13) Since (H2) and (H3) hold, by Lemma 2.5 it follows from (4.13) that lim_𝑡→+∞𝑢₁(𝑡)=(−𝑟₁+𝑎₂₁(𝑁^𝑥₁−𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),lim_𝑡→+∞𝑢₂(𝑡)=𝑏(−𝑟₁+𝑎₂₁(𝑁^𝑥₁−𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.14) By comparison, we obtain that 𝑉₁=liminf_𝑡→+∞𝑦₁(𝑡)≥(−𝑟₁+𝑎₂₁(𝑁^𝑥₁−𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑉₂=liminf_𝑡→+∞𝑦₂(𝑡)≥𝑏(−𝑟₁+𝑎₂₁(𝑁^𝑥₁−𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.15) Since these two inequalities hold for arbitrary 𝜀>0 sufficiently small, we conclude that 𝑉₁≥𝑁^𝑦1₁,𝑉₂≥𝑁^𝑦2₁, where 𝑁^𝑦1₁=(−𝑟₁+𝑎₂₁𝑁^𝑥₁)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑁^𝑦2₁=𝑏(−𝑟₁+𝑎₂₁𝑁^𝑥₁)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.16) Therefore, for 𝜀>0 sufficiently small, there exists a 𝑇₄≥𝑇₃+𝜏 such that if 𝑡>𝑇₄, 𝑦₁(𝑡)≥𝑁^𝑦1₁−𝜀, 𝑦₂(𝑡)≥𝑁^𝑦2₁−𝜀.
For 𝜀>0 sufficiently small, it follows from the first equation of system (1.6) that, for 𝑡>𝑇₄,
𝑥̇(𝑡)≤𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂(𝑁^𝑦1₁−𝜀)). (4.17) A comparison argument yields 𝑈=limsup_𝑡→+∞𝑥(𝑡)≤𝑟−𝑎₁₂(𝑁^𝑦1₁−𝜀)𝑎₁₁. (4.18) Since this is true for arbitrary 𝜀>0, we conclude that 𝑈≤𝑀^𝑥₂, where 𝑀^𝑥₂=𝑟−𝑎₁₂𝑁^𝑦1₁𝑎₁₁. (4.19) Hence, for 𝜀>0 sufficiently small, there exists a 𝑇₅≥𝑇₄ such that if 𝑡>𝑇₅, 𝑥(𝑡)≤𝑀^𝑥₂+𝜀.
Again, we derive from the second and the third equations of system (1.6) that for 𝑡>𝑇₅+𝜏,
𝑦̇₁(𝑡)≤−𝑟₁𝑦₁(𝑡)+𝑎₂₁(𝑀^𝑥₂+𝜀)𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.20) Since (H2) and (H3) hold, by Lemma 2.5, a comparison argument shows that 𝑈₁=limsup_𝑡→+∞𝑦₁(𝑡)≤(−𝑟₁+𝑎₂₁(𝑀^𝑥₂+𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑈₂=limsup_𝑡→+∞𝑦₂(𝑡)≤𝑏(−𝑟₁+𝑎₂₁(𝑀^𝑥₂+𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.21) Since these inequalities are true for arbitrary 𝜀>0 sufficiently small, we conclude that 𝑈₁≤𝑀^𝑦1₂,𝑈₂≤𝑀^𝑦2₂, where 𝑀^𝑦1₂=(−𝑟₁+𝑎₂₁𝑀^𝑥₂)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑀^𝑦2₂=𝑏(−𝑟₁+𝑎₂₁𝑀^𝑥₂)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.22) Hence, for 𝜀>0 sufficiently small, there is a 𝑇₆≥𝑇₅+𝜏 such that if 𝑡>𝑇₆, 𝑦₁(𝑡)≤𝑀^𝑦1₂+𝜀, 𝑦₂(𝑡)≤𝑀^𝑦2₂+𝜀.
For 𝜀>0 sufficiently small, it follows from the first equation of system (1.6) that, for 𝑡>𝑇₆,
𝑥̇(𝑡)≥𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂(𝑀^𝑦1₂+𝜀)). (4.23) By comparison we obtain that 𝑉=liminf_𝑡→+∞𝑥(𝑡)≥𝑟−𝑎₁₂(𝑀^𝑦1₂+𝜀)𝑎₁₁. (4.24) Since this is true for arbitrary 𝜀>0 sufficiently small, we conclude that 𝑉≥𝑁^𝑥₂, where 𝑁^𝑥₂=𝑟−𝑎₁₂𝑀^𝑦1₂𝑎₁₁. (4.25) Therefore, for 𝜀>0 sufficiently small, there is a 𝑇₇≥𝑇₆ such that if 𝑡>𝑇₇, 𝑥(𝑡)≥𝑁^𝑥₂−𝜀.
For 𝜀>0 sufficiently small, we derive from the second and the third equations of system (1.6) that for 𝑡>𝑇₇+𝜏,
𝑦̇₁(𝑡)≥−𝑟₁𝑦₁(𝑡)+𝑎₂₁(𝑁^𝑥₂−𝜀)𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.26) Since (H2) and (H3) hold, by Lemma 2.5 and by comparison, it follows from (4.26) that 𝑉₁=liminf_𝑡→+∞𝑦₁(𝑡)≥(−𝑟₁+𝑎₂₁(𝑁^𝑥₂−𝜀))(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑉₂=liminf_𝑡→+∞𝑦₂(𝑡)≥𝑏(−𝑟₁+𝑎₂₁(𝑁^𝑥₂−𝜀))(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.27) Since these two inequalities hold for arbitrary 𝜀>0 sufficiently small, we conclude that 𝑉₁≥𝑁^𝑦1₂,𝑉₂≥𝑁^𝑦2₂, where 𝑁^𝑦1₂=(−𝑟₁+𝑎₂₁𝑁^𝑥₂)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑁^𝑦2₂=𝑏(−𝑟₁+𝑎₂₁𝑁^𝑥₂)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.28) Continuing this process, we derive six sequences 𝑀^𝑥_𝑛,𝑀^𝑦1_𝑛,𝑀^𝑦2_𝑛,𝑁^𝑥_𝑛,𝑁^𝑦1_𝑛,𝑁^𝑦2_𝑛(𝑛=1,2,…) such that, for 𝑛≥2, 𝑀^𝑥_𝑛=𝑟−𝑎₁₂𝑁^𝑦1_𝑛−1𝑎₁₁,𝑀^𝑦1_𝑛=(−𝑟₁+𝑎₂₁𝑀^𝑥_𝑛)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑀^𝑦2_𝑛=𝑏(−𝑟₁+𝑎₂₁𝑀^𝑥_𝑛)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)²,𝑁^𝑥_𝑛=𝑟−𝑎₁₂𝑀^𝑦1_𝑛𝑎₁₁,𝑁^𝑦1_𝑛=(−𝑟₁+𝑎₂₁𝑁^𝑥_𝑛)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),𝑁^𝑦2_𝑛=𝑏(−𝑟₁+𝑎₂₁𝑁^𝑥_𝑛)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)². (4.29) It is readily seen that 𝑁^𝑥_𝑛≤𝑉≤𝑈≤𝑀^𝑥_𝑛,𝑁^𝑦𝑖_𝑛≤𝑉_𝑖≤𝑈_𝑖≤𝑀^𝑦𝑖_𝑛(𝑖=1,2). (4.30) We derive (4.29) that 𝑀^𝑦1_𝑛+1=𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁𝑎²₁₁𝑎²₂₂(𝑎₁₁𝑏𝜃𝜃+𝑟₂+𝑎₂₁𝑟−𝑎₁₁𝑟₁𝑎)+²¹²𝑎²₂₁𝑎²₁₁𝑎²₂₂𝑀^𝑦1_𝑛. (4.31) Noting that 𝑀^𝑦1_𝑛≥𝑦^∗₁ and 𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁, it follows from (4.31) that 𝑀^𝑦1_𝑛+1−𝑀^𝑦1_𝑛=𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁𝑎²₁₁𝑎²₂₂(𝑎₁₁𝑏𝜃𝜃+𝑟₂+𝑎₂₁𝑟−𝑎₁₁𝑟₁𝑎)+(²¹²𝑎²₂₁𝑎²₁₁𝑎²₂₂−1)𝑀^𝑦1_𝑛≤𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁𝑎²₁₁𝑎²₂₂(𝑎₁₁𝑏𝜃𝜃+𝑟₂+𝑎₂₁𝑟−𝑎₁₁𝑟₁𝑎)+(²¹²𝑎²₂₁𝑎²₁₁𝑎²₂₂−1)𝑦^∗₁=0. (4.32) Thus, the sequence 𝑀^𝑦1_𝑛 is nonincreasing. Hence, lim_𝑛→+∞𝑀^𝑦1_𝑛 exists. Taking 𝑛→+∞, it follows from (4.31) that lim_𝑛→+∞𝑀^𝑦1_𝑛=(𝜃𝑏𝑎₁₁/(𝜃+𝑟₂))+𝑎₂₁𝑟−𝑎₁₁𝑟₁𝑎₁₂𝑎₂₁+𝑎₁₁𝑎₂₂=𝑦^∗₁. (4.33) We therefore obtain from (4.30) and (4.33) that lim_𝑛→+∞𝑀^𝑦2_𝑛=𝑦^∗₂,lim_𝑛→+∞𝑁^𝑦1_𝑛=𝑦^∗₁,lim_𝑛→+∞𝑁^𝑦2_𝑛=𝑦^∗₂,lim_𝑛→+∞𝑀^𝑥_𝑛=𝑥^∗,lim_𝑛→+∞𝑁^𝑥_𝑛=𝑥^∗. (4.34) It follows from (4.30), (4.33) and (4.34) that 𝑈=𝑉=𝑥^∗,𝑈₁=𝑉₁=𝑦^∗₁,𝑈₂=𝑉₂=𝑦^∗₂. (4.35) We therefore have lim_𝑛→+∞𝑥(𝑡)=𝑥^∗,lim_𝑛→+∞𝑦₁(𝑡)=𝑦^∗₁,lim_𝑛→+∞𝑦₂(𝑡)=𝑦^∗₂. (4.36) Noting that if 𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁, then 𝑝₀−𝑞₀>2(𝜃+𝑟₂)𝑎₁₁𝑎₂₁(𝑥^∗)²>0. By Theorem 3.1, the positive equilibrium 𝐸^∗ is locally stable. We therefore conclude that 𝐸^∗ is globally stable. The proof is complete.

Theorem 4.2. If 𝑟₁(𝜃+𝑟₂)>𝑏𝜃+𝑎₂₁𝑟/𝑎₁₁(𝜃+𝑟₂), the equilibrium 𝐸₁ of system (1.6) is globally asymptotically stable.

Proof. Let (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) be any positive solution of system (1.6) with initial conditions (1.7). It follows from the first equation of system (1.6) that 𝑥̇(𝑡)≤𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)). (4.37) A standard comparison argument shows that limsup_𝑡→+∞𝑟𝑥(𝑡)≤𝑎₁₁. (4.38) Choose 𝜀>0 sufficiently small such that 𝑟₁(𝜃+𝑟₂)>𝑏𝜃+𝑎₂₁(𝑟𝑎₁₁+𝜀)(𝜃+𝑟₂). (4.39) Hence, for 𝜀>0 sufficiently small satisfying (4.39), there is a 𝑇₁>0 such that if 𝑡>𝑇₁, then 𝑥(𝑡)≤𝑟/𝑎₁₁+𝜀. We derive from the second and the third equations of system (1.6) that for 𝑡>𝑇₁+𝜏, 𝑦̇₁(𝑡)≤−𝑟₁𝑦₁(𝑡)+𝑎₂₁(𝑟𝑎₁₁+𝜀)𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.40) Consider the following auxiliary equations: 𝑢̇₁(𝑡)=−𝑟₁𝑢₁(𝑡)+𝑎₂₁(𝑟𝑎₁₁+𝜀)𝑢₁(𝑡−𝜏)−𝑎₂₂𝑢²₁(𝑡)+𝜃𝑢₂𝑢̇(𝑡),₂(𝑡)=𝑏𝑢₁(𝑡)−(𝜃+𝑟₂)𝑢₂(𝑡). (4.41) Since (4.39) holds, by Lemma 2.5 it follows from (4.41) that lim_𝑡→+∞𝑢₁(𝑡)=0,lim_𝑡→+∞𝑢₂(𝑡)=0. (4.42) By comparison, we obtain that lim_𝑡→+∞𝑦₁(𝑡)=0,lim_𝑡→+∞𝑦₂(𝑡)=0. (4.43) Hence, for 𝜀>0 sufficiently small satisfying (4.39), there is a 𝑇₂>𝑇₁ such that if 𝑡>𝑇₂,𝑦₁(𝑡)<𝑎₁₁𝜀.
For 𝜀>0 sufficiently small satisfying (4.39) and 𝑎₁₁𝜀<𝑟, we derive from the first equation of system (1.6) that for 𝑡>𝑇₂
𝑥̇(𝑡)≥𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₁𝜀). (4.44) By comparison, we obtain that liminf_𝑡→+∞𝑟𝑥(𝑡)≥𝑎₁₁−𝜀, (4.45) which, together with (4.38), yields lim_𝑡→+∞𝑟𝑥(𝑡)=𝑎₁₁. (4.46)
By Theorem 3.1, if 𝑟₁(𝜃+𝑟₂)>𝑏𝜃+𝑎₂₁𝑟/𝑎₁₁(𝜃+𝑟₂), the boundary equilibrium 𝐸₁ is locally stable. We therefore conclude that 𝐸₁ is globally stable in this case. This completes the proof.

Theorem 4.3. Let (𝐻1) hold. If (𝜃+𝑟₂)(𝑎₂₂𝑟+𝑎₁₂𝑟₁)<𝑎₁₂𝜃𝑏, then the equilibrium 𝐸₂ of system (1.6) is globally asymptotically stable.

Proof. Let (𝑥(𝑡),𝑦₁(𝑡),𝑦₂(𝑡)) be any positive solution of system (1.6) with initial conditions (1.7). We derive from the second and the third equations of system (1.6) that, for 𝑡>𝑇₁+𝜏, 𝑦̇₁(𝑡)≥−𝑟₁𝑦₁(𝑡)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.47) Consider the following auxiliary equations: 𝑢̇₁(𝑡)=−𝑟₁𝑢₁(𝑡)−𝑎₂₂𝑢²₁(𝑡)+𝜃𝑢₂𝑢̇(𝑡),₂(𝑡)=𝑏𝑢₁(𝑡)−(𝜃+𝑟₂)𝑢₂(𝑡). (4.48) Since (H1) holds, by Lemma 2.5 it follows from (4.48) that lim_𝑡→+∞𝑢₁1(𝑡)=𝑎₂₂(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₁,lim_𝑡→+∞𝑢₂𝑏(𝑡)=𝑎₂₂(𝜃+𝑟₂)(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₂. (4.49) By comparison we derive that lim_𝑡→+∞𝑦₁1(𝑡)≥𝑎₂₂(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₁,lim_𝑡→+∞𝑦₂𝑏(𝑡)≥𝑎₂₂(𝜃+𝑟₂)(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₂. (4.50) Hence, for 𝜀>0 sufficiently small, there is a 𝑇₁ such that if 𝑡>𝑇₁, 𝑦₁(𝑡)≥𝑦₁−𝜀, 𝑦₂(𝑡)≥𝑦₂−𝜀.
For 𝜀>0 sufficiently small, it follows from the first equation of system (1.6) that for 𝑡>𝑇₁,
𝑥̇(𝑡)≤𝑥(𝑡)(𝑟−𝑎₁₁𝑥(𝑡)−𝑎₁₂(𝑦₁−𝜀)). (4.51) Since (𝜃+𝑟₂)(𝑎₂₂𝑟+𝑎₁₂𝑟₁)<𝑎₁₂𝜃𝑏 and 𝜀>0 is sufficiently small, by comparison we derive that limsup_𝑡→+∞𝑥(𝑡)≤0. (4.52) We therefore have lim_𝑡→+∞𝑥(𝑡)=0. Hence, for 𝜀>0 sufficiently small, there is a 𝑇₂>𝑇₁+𝜏 such that if 𝑡>𝑇₂, 𝑥(𝑡)≤𝜀.
For 𝜀>0 sufficiently small, it follows from the first equation of system (1.6) that for 𝑡>𝑇₂+𝜏
𝑦̇₁(𝑡)≤−𝑟₁𝑦₁(𝑡)+𝑎₂₁𝜀𝑦₁(𝑡−𝜏)−𝑎₂₂𝑦²₁(𝑡)+𝜃𝑦₂𝑦̇(𝑡),₂(𝑡)=𝑏𝑦₁(𝑡)−(𝜃+𝑟₂)𝑦₂(𝑡). (4.53) Since (H1) holds, by Lemma 2.5, a comparison argument shows that lim_𝑡→+∞𝑦₁(𝑡)≤(−𝑟₁+𝑎₂₁𝜀)(𝜃+𝑟₂)+𝑏𝜃𝑎₂₂(𝜃+𝑟₂),lim_𝑡→+∞𝑦₂(𝑡)≤𝑏(−𝑟₁+𝑎₂₁𝜀)(𝜃+𝑟₂)+𝑏²𝜃𝑎₂₂(𝜃+𝑟₂)², (4.54) Together with (4.49), yields lim_𝑡→+∞𝑦₁1(𝑡)=𝑎₂₂(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₁,lim_𝑡→+∞𝑦₂𝑏(𝑡)=𝑎₂₂(𝜃+𝑟₂)(𝜃𝑏𝜃+𝑟₂−𝑟₁)=𝑦₂. (4.55) Noting that if (H1) holds and (𝜃+𝑟₂)(𝑎₂₂𝑟+𝑎₁₂𝑟₁)<𝑎₁₂𝜃𝑏, the equilibrium 𝐸₂ of system (1.6) is locally stable. We therefore conclude that 𝐸₂ is globally stable. The proof is complete.

5. Numerical Examples

In this section, we give some examples to illustrate the main results.

Example 5.1. In system (1.6), we let 𝑟=2,𝑟₁=2,𝑟₂=1,𝑎₁₁=4,𝑎₁₂=3,𝑎₂₁=2,𝑎₂₂=4,𝑏=1,𝜃=2,and𝜏=1. It is easy to show that 𝑟₁(𝜃+𝑟₂)−𝑏𝜃−𝑎₂₁𝑟(𝜃+𝑟₂)/𝑎₁₁=1>0. By Theorem 4.2 we see that the equilibrium 𝐸₁(0.5,0,0) of system (1.6) is globally stable. Numerical simulation illustrates our result (see Figure 1).

Figure 1 

The temporal solution found by numerical integration of system (1.6) with 𝑟=2,𝑟1=2,𝑟2=1,𝑎11=4,𝑎12=3,𝑎21=2,𝑎22=4,𝑏=1,𝜃=2,𝜏=1 and initial value is (0.6,0.2,0.2).

Example 5.2. In system (1.6), we let 𝑟=1,𝑟₁=1,𝑟₂=1,𝑎₁₁=4,𝑎₁₂=3,𝑎₂₁=2,𝑎₂₂=1,𝑏=4,𝜃=1,and𝜏=1. It is easy to show that 𝑏𝜃−𝑟₁(𝜃+𝑟₂)=2>0,and𝑟−𝑎12𝑦₁=−2<0. By Theorem 4.3 we see that the equilibrium 𝐸₂(0,1,2) of system (1.6) is globally stable, as depicted in Figure 2.

Figure 2 

The temporal solution found by numerical integration of system (1.6) with 𝑟=1,𝑟1=1,𝑟2=1,𝑎11=4,𝑎12=3,𝑎21=2,𝑎22=1,𝑏=4,𝜃=1,𝜏=1 and initial value is (0.6,0.2,0.2).

Example 5.3. In system (1.6), we let 𝑟=4,𝑟₁=1,𝑟₂=1,𝑎₁₁=2,𝑎₁₂=1,𝑎₂₁=2,𝑎₂₂=1.5,𝑏=4,𝜃=1,and𝜏=1. System (1.6) with above coefficients has a unique positive equilibrium 𝐸^∗(1.5,2,4). It is easy to show that 𝑎₁₁𝑎₂₂−𝑎₁₂𝑎₂₁=1>0. By Theorem 4.1 we see that the positive equilibrium 𝐸^∗(1.5,2,4) of system (1.6) is globally stable. Numerical simulation illustrates our result (see Figure 3).

Figure 3 

The temporal solution found by numerical integration of system (1.6) with 𝑟=4,𝑟1=1,𝑟2=1,𝑎11=2,𝑎12=1,𝑎21=2,𝑎22=1.5,𝑏=4,𝜃=1,𝜏=1 and initial value is (0.8,1.5,3.6).

Example 5.4. In system (1.6), we let 𝑟=4,𝑟₁=1,𝑟₂=1,𝑎₁₁=1,𝑎₁₂=2,𝑎₂₁=4,𝑎₂₂=2,𝑏=4,and𝜃=1. System (1.6) with above coefficients has a unique positive equilibrium 𝐸^∗(0.6,1.7,3.4). By computing, we can get 𝑝₀−𝑞₀=−6.48<0 and 𝜔₀=0.9326,𝜏₀=1.664. By Theorem 3.1 we see that the equilibrium 𝐸^∗ is locally stable for 𝜏<𝜏₀; and 𝐸^∗ is unstable for 𝜏>𝜏₀. And system (1.6) undergoes a Hopf Bifurcation at the positive equilibrium 𝐸^∗ when 𝜏=𝜏₀. See Figures 4 and 5.

Figure 4 

When 𝜏=1<𝜏0, the positive equilibrium 𝐸∗ is asymptotically stable. Here 𝑟=4,𝑟1=1,𝑟2=1,𝑎11=1,𝑎12=2,𝑎21=4,𝑎22=2,𝑏=4,𝜃=1 and initial value is (0.8,1.5,3.6).

Figure 5 

When 𝜏=2>𝜏0, the positive equilibrium 𝐸∗ is asymptotically stable. Here 𝑟=4,𝑟1=1,𝑟2=1,𝑎11=1,𝑎12=2,𝑎21=4,𝑎22=2,𝑏=4,𝜃=1 and the initial value is (0.8,1.5,3.6).

6. Discussion

In this paper, we considered a delayed predator-prey model with stage structure for the predator. By using the iteration technique and comparison argument, respectively, sufficient conditions were established for the global stability of the positive equilibrium and two boundary equilibria of system (1.6). By Theorems 4.1, 4.2 and 4.3, we see that: (i) If (H2) holds and 𝑎₁₁𝑎₂₂>𝑎₁₂𝑎₂₁, system (1.6) is permanent. (ii) If 𝑟₁(𝜃+𝑟₂)>𝑏𝜃+𝑎₂₁𝑟/𝑎₁₁(𝜃+𝑟₂), the prey species is persistent but the predator becomes extinct. (iii) If (H1) holds and (𝜃+𝑟₂)(𝑎₂₂𝑟+𝑎₁₂𝑟₁)<𝑎₁₂𝜃𝑏, the predator species is persistent but the prey species becomes extinct.

Acknowledgments

The authors wish to thank the reviewers for their valuable comments and suggestions that greatly improved the presentation of this work. This work was supported by the National Natural Science Foundation of China (no. 10671209) and the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.