DYNAMICS OF A STAGE-STRUCTURED POPULATION MODEL WITH A STATE-DEPENDENT DELAY

. This paper is devoted to the dynamics of a predator-prey model with stage structure for prey and state-dependent maturation delay. Firstly, positivity and boundedness of solutions are addressed to describe the population survival and the natural restriction of limited resources. Secondly, the existence, uniqueness, and local asymptotical stability of (boundary and coexisting) equilibria are investigated by means of degree theory and Routh-Hurwitz criteria. Thirdly, the explicit bounds for the eventual behaviors of the mature population are obtained. Finally, by means of comparison principle and two auxiliary systems, it observed that the local asymptotical stability of either of the positive interior equilibrium and the positive boundary equilibrium implies that it is also globally asymptotical stable if the derivative of the delay function around this equilibrium is small enough.


1.
Introduction. Time-delay in a natural ecosystem has been widely considered for a long time (see e.g. [23,35,36,37]). Recently, stage-structured models have received great attention. For example, Gurney, Blythe, and Nisbet [14] proposed a time delay growth model of blowflies and verified solutions of the time delay model are in accordance with the data in blowflies growth experiments by Nicholson (see [28]). This forces researchers to introduce time delay and stage structure into population research (see e.g. [8,9,10,15,16,21,22,31,38,40]). The work of Aiello and Freedman [2] on a single species growth model with stage structure represents a mathematically more careful and biologically meaningful formulation approach. This model assumes an average age to maturity which appears as a constant time delay reflecting a delayed birth of immatures and a reduced survival of immatures to their maturity. The model takes the forṁ v(t) = αu(t) − γv(t) − αe −γτ u(t − τ ), where v(t) denotes the immature population density, u(t) represents the mature population density, α > 0 represents the birth rate, γ > 0 is the immature death rate, β > 0 is the mature death and overcrowding rate, τ is the time to maturity. The term αe −γτ u(t − τ ) represents the immatures who were born at time t − τ and survive at time t (with the immature death rate γ), and therefore represents the transformation of immatures to matures. Following the way of Aiello and Freedman [2], a number of authors studied different kinds of stage-structured models and carried out some significant works (see e.g. [3,4,5,12,30,32,33,34]). Due to the influence of circumstances such as resources and interaction, however, the constant time delay is not reasonable any more (see e.g. [1,5,11,13,18,42,43]). Motivated by their knowledge about whale and seal populations, Aiello, Freedman and Wu [3] considered a system with a state-dependent delay and two stages, and proposed the following state-dependent delay model where the state-dependent time delay τ (v + u) is taken to be an increasing differentiable function of the total population (v + u) so that τ (v + u) ≥ 0 and τ m ≤ τ (v + u) ≤ τ M with τ (0) = τ m and τ (+∞) = τ M . These assumptions imply that the maturation time for the species depends on the total number of them (matures plus immatures) around. The term αe −γτ (v+u) u(t − τ (v + u)) appearing in both equations of system (2) represents the density of individuals survive to leave the immature and just enter the mature class. Lotka-Volterra type predator-prey systems are very important in the models of multi-species populations interactions and have been studied widely. For a number of animals, it seems reasonable to assume that the predator population feed on the mature prey because immature prey population are concealed in the mountain caves and are raised by their parents, and that the rate of predators attack at immature preys can be ignored (see e.g. [26]). Song and Chen [33] considered the exploitation of a predator-prey population with stage structure and harvesting for prey, and proposed the following system of delay differential equations: where X 1 (t) and X 2 (t) represent the immature population density and mature prey population density, respectively; Y (t) represents the density of predator population; a is the transformation coefficient of mature predator population; τ represents the transformation of immatures to matures; αe −γτ X 2 (t − τ ) represents the immatures who were born at t − τ and survive at t with the immature death rate γ; α is the birth rate of the immature prey population; β represents the mature death and overcrowding rate; E is the harvesting effort. Song and Chen [33] investigated the global asymptotical stability of three nonnegative equilibria and a threshold of harvesting for the mature prey population. The effect of the delay on the populations at positive equilibrium and the optimal harvesting of the mature prey population are also considered. Since then, more and more researchers have worked on two-species predator-prey systems with different stage structure on predator or prey, see, for example, [19,20,24,25,30,32,34,39]. Furthermore, Qu and Wei [30] investigated the stability and Hopf bifurcation of the interior equilibrium. Considering a window in maturation time delay parameter that generates sustainable oscillatory dynamics, Gourley and Kuang [12] formulated a general and robust predator-prey model with stage-structure and constant maturation time delay. And then by following Gourley and Kuang [12], further results have been investigated (see e.g. [4]).
Motivated by the works of Aiello, Freedman, and Wu [3] and Song and Chen [33], in this paper, we are concerned with the effect of stage structure for prey and state-dependent delay on Lotka-Volterra type predator-prey system. To do so, we study the following differential systeṁ where α, β, γ, a, b, c, r are all positive constant, the state-dependent time delay τ (X) is taken to be an increasing, concave downwards, twice differentiable function on [0, +∞) of the total population X = X 1 + X 2 , so that τ (X) ≥ 0 and τ (X) ≤ 0 for all X ≥ 0, and we shall also assume that which is the number of immatures that have survived to time t = 0. Here, τ s is the maturation time of the prey at t = 0, and the lower limit on the integral is −τ s because any prey born before that time will have matured before time t = 0. It follows that τ s is given by τ s = τ (X 1 (0) + X 2 (0)), i.e., Note that τ s appears on both the left-and right-hand sides of the above equation, so that τ s is determined implicitly. For our model to make sense, i.e., to exclude the possibility of adults becoming immatures except by birth, we need to find conditions ensuring that t − τ (X(t)) is an increasing function of t as t increases. Namely, we need τ (X 1 + X 2 )(Ẋ 1 +Ẋ 2 ) < 1, which is equivalent to Compared with [3] and [33], our theoretical results are more accurate and our method is completely different. More precisely, we shall employ degree theory and Routh-Hurwitz criteria to investigate the existence and linearized stability of equilibria. Moreover, we find out the relationship among uniqueness, local asymptotical stability and global asymptotical stability of equilibria. Two auxiliary systems and comparison principles are introduced to prove the global asymptotical stability of nonnegative equilibria. This paper is organized as follows: the positivity and boundedness of all solutions of (4) are obtained in section 2. In section 3, we first investigate the existence of positive coexisting equilibria by means of degree theory. Then we discuss the uniqueness of positive coexisting equilibria. Section 4 is devoted to the stability of equilibria, especially the positive ones, through Routh-Hurwitz criteria. In section 5, we investigate the state bounds on the eventual behaviour of X 2 (t) and also introduce two auxiliary systems in order to overcome the difficulties caused by the presence of state-dependent delay during the investigation of the global stability of system (4). Section 6 is devoted to the global asymptotical stability. Moreover, we illustrate our results with some numerical simulations in Section 7.
2. Positivity and boundedness. Since the solutions of system (4) represent populations and we also anticipate that limited resources will place a natural restriction to how many individuals can survive, we need address positivity and boundedness of the solution of the solution of the system.
Then we complete the proof of this theorem.
It follows from (4) thaṫ ) is a solution of the following system with the initial value x 1 (s) = X 1 (s) and x 2 (0) = X 2 (0) for −τ M ≤ s ≤ 0, It follows from Theorem 5.1 in Aiello, Freedman and Wu [3] that x 1 (t) and x 2 (t) are bounded, and that we have the following theorem immediately. Proof. It follows from Theorem 2.
It follows from Theorem 2.3 that for every sufficiently small positive constant ε there exists T > 0 such that If βr > αce −γτm then ε can chosen small enough such that A ε ≤ 0 and hence lim t→∞ Y (t) = 0. If βr ≤ αce −γτm then A ε > 0 for all ε, and hence lim sup t→∞ Y (t) ≤ A ε /b. Since this is true for any ε > 0, it follows that lim sup t→∞ Y (t) ≤ A 0 /b. This completes the proof.
Proof. It follows from Theorem 2.1 that X 2 (t) > 0 for all t > 0. Now, we show that there exists a positive constant δ depending on ϕ such that X 2 (t) > δ for all t ≥ 0. Assume on the contrary that lim inf t→∞ X 2 (t) = 0. One case is that X 2 (t) is eventually monotonic decreasing and tends to 0 as t → ∞. In this case, for all t > T + τ M . This is a contraction. The other case is that X 2 (t) is eventually oscillatory and lim inf t→∞ X 2 (t) = 0. Suppose that there exists {t k } such that t k − t k−1 > τ M for all k ∈ N, that X 2 (t) achieves a local minimum and X 2 (t) ≥ X 2 (t k ) for all t ∈ (t 1 , t k ), and that lim k→∞ X 2 (t k ) = 0. If βr > αce −γτm then there exists K > 0 such that X 2 (t k ) < α 2β e −γτ M and Y (t k ) < α 2a e −γτ M for all k ≥ K and hence for all k > K. This is a contraction. If βr ≤ αce −γτm then for every 0 < ε < This is a contraction. Thus, we complete the proof of this theorem. Theorem 2.5 implies that for each given positive initial function ϕ, the mature population X 2 (t) is uniformly bounded away from zero. We have proved that X 2 remains positive and bounded with the given initial condition. Similar to Aiello, Freedman and Wu [3], we shall give the positivity of X 1 , which depend on the initial condition and our having a strictly positive lower bound δ and an upper bound ∆ for X 2 .

Remark 1.
To prove X 1 (t) > 0 for values of t, it seems to be impossible without placing additional restrictions on either the initial conditions or on the delay τ (X). For example, if τ (X) ≡ 0, it has been shown in [2] that X 1 (t) is positive for all t.
In Theorem 2.6, we give a set of initial conditions on τ (X), while maintaining the essential character of the state-dependent time delay. The following is a corollary of Theorem 2.6, both of them are only sufficient conditions for the positivity of X 1 (t).
3. Existence of equilibria. The purpose of this section is to investigate the existence of equilibria (X 1 , X 2 , Y ) of system (4), which satisfy It follows from the first two equations of (5) that X 1 = (X 2 , Y ), where : R 2 → R is defined as for all x, y ∈ R. Thus, system (5) can be reduced to Our main interest is the existence and uniqueness of positive equilibria E(X 1 , X 2 , Y ) with Y = (cX 2 −r)/b and (X 1 , X 2 ) satisfying X 1 = g(X 2 ) and f (X 2 ) = 0, where g: R → R and f : [0, ∞) → R are defined by respectively. Thus, it suffices to investigate the existence and uniqueness of positive zero points of f . Note that Thus, each zero x * of f satisfies 0 < x * ≤ bαe −γτm +ar bβ+ac . Consider a set Ω = x ∈ R 0 < x < bα (1 + e −γτm ) + ar bβ + ac .

DYNAMICS OF A STAGE-STRUCTURED POPULATION MODEL 3529
Note that for all x ∈ Ω. This implies that the graph of the curve y = f (x) is concave upwards on Ω. In what follows, we shall employ the Brouwer degree to prove the existence of positive zero points of function f .
It turns out that H(t, x) is an Ω-admissible homotopy. By the homotopy invariance, Note that deg(G, Ω) = −1. This implies that f has at least one zero point X * 2 in Ω. Finally, it follows from f (x) > 0 for all x ∈ Ω that f has exactly one or two different zeros in Ω. In fact, if f has exactly two different zeros in Ω, then deg(f, Ω) = 0, which contradicts deg(f, Ω) = −1. Therefore, f has exactly one zero in Ω.
In view of Lemma 3.1, system (4) has at least one nontrivial equilibrium ( Theorem 3.2. Assume that βr < αce −γτ M , then system (4) has exactly one posi- where Finally, we consider the existence of boundary equilibria. For an equilibrium (X 1 , X 2 , Y ) of system (4), if X 2 = 0 then X 1 = 0. It is clear that system (4) has at least two boundary equilibria E 0 (0, 0, 0) and E 1 (X 10 , X 20 , 0) with X 10 = (X 20 , 0) and It is easy to see that each zero This implies that the graph of the curve y = h(u) is concave upwards on Ω 1 . Using a similar argument as the proof of Lemma 3.1, we can conclude that h has exactly one zero X 20 in Ω 0 , and hence that system (4) has a nontrivial equilibrium (X 10 , X 20 , 0), where X 10 = (X 20 , 0). It follows from βX 20 = αe −γτ (X10+X20) < α that 0 < X 20 < α/β, which implies that X 10 > 0. Thus, we have the following result.
In view of h(X 20 ) = 0, using a similar argument as the proof of Theorem 3.1, we have h (X 20 ) < 0 and hence 4. Linearized stability. Linearizing (4) is not completely straightforward because the delay is a function depending on the state variables X 1 and X 2 . It was shown in [7,17] that generically the behaviour of the state-dependent delay except for its value has no effect on the stability of an equilibrium, and that a local linearization is valid by treating the delay function as a constant at the equilibrium point. Hence to study the local stability of an equilibrium E(X 0 1 , X 0 2 , Y 0 ) of (4), we linearize (4) at E(X 0 1 , X 0 2 , Y 0 ) by treating the delay τ (X 1 (t) + X 2 (t)) as τ (X 0 1 + X 0 2 ). The resulting linear system is a differential equation with a constant delay: for This leads to the following characteristic equation For the extinction equilibrium E 0 (0, 0, 0), (12) reduces to Clearly λ = −γ and λ = −r are two of the eigenvalues. All other eigenvalues λ satisfy the equation λ = αe −(γ+λ)τm , which always has a solution with positive real part. Hence E 0 is unstable.
For the boundary equilibrium E 1 (X 10 , X 20 , 0) with its existence guaranteed by Theorem 3.3, the eigenvalues λ are the roots of the equation 20 and X 0 = X 10 + X 20 . We have the following result.
Theorem 4.1. The boundary equilibrium E 1 (X 10 , X 20 , 0) is locally asymptotically stable if and only if X 20 < r/c.
Proof. Since one eigenvalue is λ = cX 20 − r, we only need to investigate the eigenvalues given by It suffices to show that all zeros of (13) have negative real parts.

SHANGZHI LI AND SHANGJIANG GUO
For such ζ to exist, (15) must have real roots v. After substituting for ζ and rearranging, we see that v is a zero of the following functioñ It follows from (10) that ρ < 1 and hence thath(v) > 0 for all v ∈ R. Therefore, (15) has no real solutions. Thus, (13) has no purely imaginary zeros when τ (X 0 ) > 0. Note that all zeros of (13) have negative real parts when τ (X 0 ) = 0. Then by continuity, all zeros of (13) have negative real parts when τ (X 0 ) ≥ 0. (13) is exactly the characteristic equation described by Aiello, Freedman and Wu [3], who obtained a set of sufficient conditions ensuring that all zeros of (13) have negative real parts. In view of the proof of Theorem 4.1, we see that the sufficient conditions proposed by Aiello, Freedman and Wu [3] are unnecessary and can be removed.

Remark 2. Note that
It is very interesting to investigate the stability of the boundary equilibrium E 1 (X 10 , X 20 , 0) when βr < αce −γτ M . In this case, we have ch(r/c) > αce −γτ M − βr > 0 and hence X 20 > r/c. In view of Theorem 4.1, we have the following result.
Corollary 2. Assume that βr < αce −γτ M , then the boundary equilibrium E 1 (X 10 , X 20 , 0) is unstable. Now, we consider the stability at the positive equilibrium E * (X * 1 , X * 2 , Y * ) under the assumption that βr < αce −γτ M . The characteristic equation (12) can be rewritten as It is helpful to consider τ (X * ) (or equivalently, ξ * ) as a parameter. If the stability of the positive equilibrium point switches as the value of τ (X * ) (or equivalently, ξ * ) increases starting from 0, there must be a value of τ (X * ) (respectively, ξ * ) at which there are eigenvalues on the imaginary axis.
Note that E * is locally asymptotically stable when τ (X * ) = 0 (or equivalently, ξ * ). As τ (X * ) (or equivalently, ξ * ) increases, the stability of the steady state can be lost only if purely imaginary roots of (16) appear. Substituting λ = i √ w (w > 0) into (16) and then separating the real and imaginary parts, we have Similarly, we square and add the above two equations and obtain that where then we have F (w, 0) = (w + γ 2 )(w 2 + a * w + b * ). It follows from the discussion about (18) that F (w, 0) > 0 for all w ≥ 0. Note that It follows from Routh-Hurwitz criteria [27] that no positive solution to (20) exists. Thus no stability switches can occur. This completes the proof of the theorem.

Global behaviors.
In order to investigate the global asymptotical stability of the positive interior equilibrium E * (X * 1 , X * 2 , Y * ) and the positive boundary equilibrium E 1 (X 10 , X 20 , 0) of (4), we first investigate the state bounds on the eventual behaviour of X 2 (t). We introduce the following two quantities depending on a: It is easy to see that X − a is increasing (respectively, decreasing) with respect to a if βr > αce −γτ M (respectively, βr < αce −γτ M ), and that X + a is increasing (respectively, decreasing) with respect to a if βr > αce −γτm (respectively, βr < αce −γτm ).
The following result gives the state bounds on the eventual behaviour of X 2 (t), independent of admissible initial conditions.
Then 0 < ξ < ∞ and lim sup t→∞ X 2 (t) = ξ. We claim that ξ ≤ X + a . In fact, if cξ ≤ r then it follows from (25) that ξ ≤ r/c < X − a ≤ X + a . If cξ > r then lim k→∞ Y (t k ) = (−r + cξ)/b. Suppose on the contrary that ξ > X + a . We now choose a subsequence of {t k }, relabelled as {t k } such that t k+1 ≥ t k + τ M and X 2 (t k ) → ξ as k → ∞. Let X = lim sup k→∞ X(t k ), where X(t) = X 1 (t) + X 2 (t). We then choose a further subsequence of {t k }, again relabelled by {t k } such that lim k→∞ X(t k ) = X. Now let ξ = lim sup k→∞ X 2 (t k − τ (X(t k ))) for this subsequence {t k }. We choose a final subsequence of {t k }, once again relabelled by {t k }, such that lim k→∞ X 2 (t k − τ (X(t k ))) = ξ . It follows from the definition of ξ that ξ ≤ ξ. Then from (4) and Lemma 5.7, taking the limit as k → ∞, we obtain a contradiction. Therefore, ξ ≤ X + a . This completes the proof. In order to overcome the difficulties caused by the presence of state-dependent delay during the investigation of the global stability of system (4), we introduce two auxiliary systems, one of which takes the form Here, τ (·) is the same to that in system (4), k ≥ 0, α, β, γ are positive constants, s, d ∈ [k − , k + ], and Note that (26) is a mixed quasi-monotone system (see [29,41]), then we have the following observations. Proof. Obviously, the equilibria (X 1 , X 2 ) of system (26) satisfy where (·, k): Note that for all x ∈ [0, k − ], and Thus, there exists some X 2 (s, d, k) ∈ (k − , k + ) such that F( X 2 (s, d, k), s, d, k) = 0. It follows from (27) and hence Therefore, system (26) has a positive equilibrium point ( X 1 (s, d, k), X 2 (s, d, k)).

SHANGZHI LI AND SHANGJIANG GUO
then lim n→∞ B u n = lim n→∞ B l n = η(k). Proof. It follows from Lemma 5.4(i)(iv) that k − < B l n < B u n < k + for all n ∈ N. Let . This completes the proof of this lemma.
The following result is trivial and the proof can be found in Chen [6].
Lemma 5.7. Consider the following logistic equation The other auxiliary system is where γ are positive constants, k ≥ 0, and function τ (·) is the same to that of system (4). We have the following result on the existence, uniqueness, and global attractivity of a positive equilibrium point of system (30). Assume further that (29) holds, then the positive equilibrium point (ζ(k), η(k)) is unique and attracts all of the positive solutions of system (30). In particular, ζ(0) = X 10 and η(0) = X 20 .
Proof. It follows from the proof of Lemma 5.4 that we can obtain the existence and uniqueness of the positive equilibrium point (ζ(k), η(k)). In what follows, we only need to prove the global attractivity of the positive equilibrium point (ζ(k), η(k)).
Using a similar argument as that of [3], we see that then for all large enough t, we have Thus, we have where (x, y, x, y) is a solution to the following system It follows from (5.3) that Thus, the positive equilibrium point (ζ(k), η(k)) attracts all of the positive solutions of system (30). This completes the proof of this lemma. 6. Global asymptotical stability. In this section, we shall investigate the global asymptotical stability of the positive equilibrium E * (X * 1 , X * 2 , Y * ) and the positive boundary equilibrium E 1 (X 10 , X 20 , 0) of (4). Theorem 3.2 says that system (4) has exactly one positive interior equilibrium E * if βr < αce −γτ M . In Theorem 4.2, we see that this equilibrium E * is locally asymptotically stable. Furthermore, the following theorem says that this equilibrium E * is also globally asymptotically stable if we further assume that the derivative of the delay function is small enough. Theorem 6.1. Assume that βr < αce −γτ M and αe −γτ (X − a ) τ (X − a )X + a < 1 2 , then positive equilibrium E * (X * 1 , X * 2 , Y * ) is globally asymptotically stable.
Thus, there exists t 1 > t 0 + τ M such that for all t > t 1 , . By comparison principle, X 2 (t) ≥ y(t) and X 1 (t) ≥ x(t) for all t > t 1 , where (x(t), y(t)) is the solution to the following equatioṅ with the initial values x(t 1 ) = X 1 (t 1 ) and . Thus, it follows from Lemma 5.8 that X l 1 ≥ ζ(ε) and X l 2 ≥ η(ε). Since this is true for any sufficiently small ε, it follows that It follows from Theorem 2.2 thaṫ X 2 (t) < αe −γτ (X) X 2 (t − τ (X)) − βX 2 2 (t). Also, by comparison principle and Lemma 5.8, we get It follows from (36) and (37) that This completes the proof. Theorem 6.2 means that the boundary equilibrium E 1 is globally asymptotically stable if βr > αce −γτm and the derivative of the delay function is small enough. It would be very interesting to see what happens to the case where αce −γτ M < βr ≤ αce −γτm . In this case, we have X − 0 ≤ X − a ≤ X + a ≤ X + 0 and X − 0 < r c ≤ X + 0 . We have the following results. τ (X − 0 )X + 0 < 1 2 , then positive boundary equilibrium E 1 (X 10 , X 20 , 0) is globally asymptotically stable.
Proof. It follows from Theorems 3.3 and 4.1 that E 1 is unique and locally asymptotically stable. So we only need to prove the global attractivity of E 1 . Define X u 1 , X u 2 , X l 1 X l 2 , and and Y u as those in the proof of Theorem 6.2. Similarly, we have X u 1 ≤ X 10 , X u 2 ≤ X 20 < r/c.
Thus, there exists t 0 > τ M such that X 2 (t) < r/c for all t > t 0 . Thus, we havė Y (t) < −bY 2 (t). In view of Lemma 5.7, we have Y u = 0 and hence lim t→∞ Y (t) = 0.
Thus, for every given ε there exists t 1 > t 0 + τ M such that 0 < Y (t) < ε for all t > t 1 , and hence that for all t > t 1 , . Using a similar argument as the proof of Theorem 6.2, we have It follows from (38) and (40) that lim t→∞ X i (t) = X i0 , i = 1, 2 This completes the proof. Theorems 6.2 and 6.3 means that if the derivative of the delay function is small enough then the local asymptotical stability of the boundary equilibrium E 1 implies that it is also globally asymptotically stable. It is easy to see that {(X 1 , X 2 , 0) | X 1 ≥ 0, X 2 ≥ 0} is a positively invariant set of system (4), that is to say, every trajectory of system (4) with the initial conditions in {(X 1 , X 2 , 0) | X 1 ≥ 0, X 2 ≥ 0} always stays in {(X 1 , X 2 , 0) | X 1 ≥ 0, X 2 ≥ 0}. Furthermore, using a similar argument as the proof of Theorem 6.2, we see that E 1 attracts every solution with initial conditions in {(X 1 , X 2 , 0) | X 1 ≥ 0, X 2 ≥ 0}. Thus, we conjecture that some solutions of system (4) can fluctuate around the positive boundary equilibrium E 1 when E 1 is unstable.
In view of Corollary 3, we see that system (4) with τ (·) ≡ τ undergoes transcritical bifurcation at βr = αce −γτ . Namely, as βr − αce −γτ increases and passes through 0, the two equilibria E * and E 1 of system (4) with τ (·) ≡ τ collide and exchange their stability. However, when βr > αce −γτ , the third component of E * of system (4) with τ (·) ≡ τ is negative and hence this kind of equilibria should be ignored. In view of Corollary 3, it is important to find that the local asymptotical stability of each nonnegative equilibrium (either E * or E 1 ) implies that it is also globally asymptotically stable. 7. Conclusions and simulations. The mathematical model we have proposed consists of three nonlinear ordinary differential equations, corresponding to an immature population, a mature population, and their predator. The predator feed only on the mature prey. We have established the stability conditions for nonnegative equilibria. It is also observed that the system is locally asymptotically stable around the positive interior equilibrium E * (respectively, the positive boundary equilibrium E 1 ) if βr < αce −γτ M (respectively, cX 20 < r). In particular, if the derivative of the delay function is small enough, then the local asymptotical stability of either of the positive interior equilibrium E * and the positive boundary equilibrium E 1 implies that it is also globally asymptotically stable. Let us now give some numerical simulations to illustrate our theoretical results. Consider system (4) with α = 2.8, γ = 0.3, β = 0.4, a = 0.5, and b = 0.1. Namely, we consider the following systemẊ 1 (t) = 2.8X 2 (t) − 0.3X 1 (t) − 2.8e −0.3τ (X) X 2 (t − τ (X)), X 2 (t) = 2.8e −0.3τ (X) X 2 (t − τ (X)) − 0.4X 2 2 (t) − 0.5X 2 (t)Y (t), Y (t) = Y (t) [−r + cX 2 (t) − 0.1Y (t)] , where the delay function is τ (x) = 4 − 2e −0.1x . It is easy to check that τ (x) > 0, τ (x) < 0, τ m = 2, and τ M = 4. We first study system (41) with r = 0.1 and c = 1. Note that τ (x) < 4β/α 2 , βr < αce −γτ M .