# The Dynamic Complexity of a Holling Type-IV Predator-Prey System with Stage Structure and Double Delays

Authors: Yakui Xue; Xiafeng Duan
Journal: Discrete Dynamics in Nature and Society
 (2011)

---

## Abstract

We invest a predator-prey model of Holling type-IV functional response with stage structure and double delays due to maturation time
for both prey and predator. The dynamical behavior of the system is investigated from the point of view of stability switches aspects. We assume that the
immature and mature individuals of each species are divided by a fixed age,
and the mature predator only attacks the mature prey. Based on some comparison arguments, sharp threshold conditions which are both necessary and
sufficient for the global stability of the equilibrium point of predator extinction
are obtained. The most important outcome of this paper is that the variation of
predator stage structure can affect the existence of the interior equilibrium point
and drive the predator into extinction by changing the maturation (through-stage) time delay. Our linear stability work and numerical results show that
if the resource is dynamic, as in nature, there is a window in maturation time
delay parameters that generate sustainable oscillatory dynamics.


---

## Body

## 1. Introduction

Predator-prey models are arguably the most fundamental building blocks of the any bio- and ecosystems as all biomasses are grown out of their resource masses. Species compete, evolve and disperse often simply for the purpose of seeking resources to sustain their struggle for their very existence. Their extinctions are often the results of their failure in obtaining the minimum level of resources needed for their subsistence. Depending on their specific settings of applications, predator-prey models can take the forms of resource-consumer, plant-herbivore, parasite-host, tumor cells (virus)-immune system, susceptible-infectious interactions, and so forth. They deal with the general loss-win interactions and hence may have applications outside of ecosystems. When seemingly competitive interactions are carefully examined, they are often in fact some forms of predator-prey interaction in disguise. The dynamic relationship between predators and their prey has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology [1] due to its universal existence and importance [2]. These problems may appear to be simple mathematically at first sight, but they are, in fact very challenging and complicated. There are many different kinds of predator-prey models in the literature [3–5]; for more details, we can refer to [2, 6]. In general, a predator-prey system may have the form (1.1)dx(t)dt=rx(1-xk)-φ(x)y,dy(t)dt=y(μφ(x)-D),
where φ(x) is the functional response function, which reflects the capture ability of the predator to prey. For more biological meaning, the reader may consult Freedman [6], May [7], and Murray [8].Some experiments and observations indicate that a nonmonotonic response occurs at this level: when the nutrient concentration reaches a high level, an inhibitory effect on the specific growth rate may occur. To model such an inhibitory effect, Andrews [9] proposed the response function φ(x)=mx/(a+bx+x2), called the Monod-Haldane function, and also called a Holling type-IV function.In the past several decades, the predator-prey systems play an important role in the modeling of multispecies population dynamics [10, 11]. Many models of population growth were studied with time delays [12–17]. Some other age- and stage-structured models of various types (discrete and distributed time delays, stochastic, etc.) have been utilized [18–23]. In the pioneering work [23], a stage-structured model of population growth consisting of immature and mature individuals was proposed, where the stage-structure was modeled by the introduction of a constant time delay, reflecting a delayed birth of immature and a reduced survival of immature to their maturity. The model takes the form (1.2)dxi(t)dt=αxm(t)-γxi(t)-αe-γτxm(t-τ),dxm(t)dt=αe-γτxm(t-τ)-βxm2(t),
where xi(t) and xm(t) represent the immature and mature populations densities, respectively, to model stage-structured population growth. There, α>0 represents the birth rate, γ>0 is the immature death rate, β>0 is the mature death and overcrowding rate, and τ is the time to maturity. The term αe-γτxm(t-τ) represents the immature who were born at time t-τ and survive at time t with the immature death rate γ and thus represents the transformation of immature to mature.Motivated by the above important works [23–25], in the present paper, we consider the following stage-structured predator-prey system with Holling type-IV functional response, which takes the form: (1.3)dxj(t)dt=bx-d1xj-be-d1τ1x(t-τ1),dx(t)dt=be-d1τ1x(t-τ1)-ax2-lxyx2+dx+e,dyj(t)dt=mxyx2+dx+e-me-d2τ2x(t-τ2)y(t-τ2)x2(t-τ2)+dx(t-τ2)+e-d2yj,dy(t)dt=me-d2τ2x(t-τ2)y(t-τ2)x2(t-τ2)+dx(t-τ2)+e-d3y,
where τ=max{τ1,τ2}, x(θ),y(θ)>0 is continuous on -τ≤θ≤0, and xj(0),x(0), yj(0)>0, y(0)>0, x and y represent mature prey and the predator densities, respectively. And xj and yj represent the immature or juvenile prey and the predator densities. The constant b is the birth rate of the mature prey. We assume that immature prey suffer a mortality rate of d1 and take τ1 units of time to mature, thus e-d1τ1 is the surviving rate of each immature prey to reach maturity. The constant d2 is the death rate of the juvenile predator and take τ2 units of time to mature, thus e-d2τ2 is the surviving rate of each immature predator to reach maturity. d3 is the death rate of the mature predator. The mature predator consumes the mature prey with functional response of Holling type-IV lx/(x2+dx+e). It is assumed in model (1.3) that the predator population only feeds on the mature prey and immature individual predators do not feed on prey and do not have the ability to reproduce. Obviously, all the constants are positive for their biological sense.For the continuity of the solutions to system (1.3), in this paper, we require(1.4)xj(0)=b∫-τ10(ed1sx(s))ds.
By the first equation of system (1.3), the initial conditions (1.4), and the arguments similar to Lemma  3.1 in [26, page 672], we have(1.5)xj(s)=b∫-τ10(ed1tx(t+s))dt,
that is, xj(t) is completely determined by x(t). Thus, the following system can be separated from system (1.3)(1.6)dx(t)dt=be-d1τ1x(t-τ1)-ax2-lxyx2+dx+e,dyj(t)dt=mxyx2+dx+e-me-d2τ2x(t-τ2)y(t-τ2)x2(t-τ2)+dx(t-τ2)+e-d2yj,dy(t)dt=me-d2τ2x(t-τ2)y(t-τ2)x2(t-τ2)+dx(t-τ2)+e-d3y,
where x(θ),yj(θ),y(θ)≥0 are continuous on -τ≤θ≤0, and x(0),yj(0),y(0)>0.Notice that, mathematically, no information on the past history ofyj is needed for system (1.6). The dynamics of model (1.6) are determined by the first equation and the third equation. Therefore, in the rest of this paper, we will study the following: (1.7)dx(t)dt=be-d1τ1x(t-τ1)-ax2-lxyx2+dx+e,dy(t)dt=me-d2τ2x(t-τ2)y(t-τ2)x2(t-τ2)+dx(t-τ2)+e-d3y.In our study, we assume that the initial conditions of system (1.7) take the form(1.8)x(θ)=ϕ(θ)≥0,y(θ)=ψ(θ)≥0,θ∈[-τ,0],ϕ(0)>0,ψ(0)>0,
the Banach space of continuous functions mapping the interval [-τ,0] into R+02, where R+02={(x,y):x≥0,y≥0}.Remark 1.1.
Becausexj(t) is completely determined by x(t) and yj(t) is completely determined by x(t) and y(t), we can get all the dynamical behaviors at the equilibria of system (1.3). Hence, we only study the system (1.7) in the following sections.In the present paper, we present a qualitative analysis for the predator-prey system (1.7) by incorporating stage structures for both prey and predator. The main goal of this paper is to study the combined effects of the stage structure on prey and predator on the dynamics of the system. The rest of the paper is organized as follows. In the next section, we present some important lemmas. In Section 3, we get all the equilibria and their feasibility. In Section 4, both necessary and sufficient for the global stability of the boundary equilibrium is established. In Section 5, the stability switches of the coexistence equilibrium of system (1.7) are gotten. Finally, we numerically illustrate our results and obtain very rich dynamics of our model. The paper ends with a discussion.


## 2. Preliminary Analysis

To prove the main results, we need the following lemmas. In the biology significance, we only study inR+2={(x,y)∈R2∣x≥0,y≥0}.Using the similar arguments to Lemma  1 in [27], we directly have.Lemma 2.1.
System (1.7) with initial conditions x(t),y(t)>0(-τ≤t≤0) and x(0),y(0)>0 has strictly positive solutions for all t>0.Lemma 2.2 (see [24, 27]).
For equation(2.1)x′(t)=bx(t-τ)-a1x(t)-a2x2(t),
where a1≥0, a2,b,τ>0, x(0)>0 and x(t)>0  for all -τ≤t≤0, one has(i)
ifb>a1, then limt→+∞x(t)=(b-a1)/a2,(ii)
ifb<a1, then limt→+∞x(t)=0.Lemma 2.3 (see [21]).
Ifa<b, then the solution of the equation
(2.2)u′(t)=au(t-τ)-bu(t),
where a,b,τ>0, and u(t)>0 for -τ≤t≤0, satisfies
(2.3)lim⁡t→∞u(t)=0.Lemma 2.4.
Positive solutions of System (1.7) with initial conditions x(t),y(t)>0(-τ≤t≤0) and x(0),y(0)>0 are ultimately bounded.Proof.
Let(x(t),y(t)) be a positive solution of system (1.7) with initial conditions x(t),y(t)>0(-τ≤t≤0) and x(0),y(0)>0. It follows from the first equation of system (1.7) that
(2.4)ẋ(t)≤be-d1τ1x(t-τ1)-ax2.
By Lemma 2.2, a comparison argument shows that
(2.5)limsup⁡t→+∞x(t)≤be-d1τ1a≤M.
It implies that x(t) is ultimately bounded. No less of generality, we suppose that there exists T1>0 and M>a/b such that x(t)<M for all t>T1. Define W(t)=(1/l)x(t)+(1/m)y(t), we get
(2.6)Ẇ(t)≤be-d1τ1lx(t-τ1)-ax2(t)l.
By Lemma 2.2 and the comparison theorem, we get
(2.7)Ẇ(t)≤be-d1τ1a,
which implies that there exists a constant C>be-d1τ1/a>0, such that all trajectories initiating in R+2 enter the region Ω:={(x,y)∈R+2∣(x/l)+(y/m)≤C+ε for any ε>0}, proving Lemma 2.4.


## 3. Equilibria and Their Feasibility

Apart from the zero solution, system (1.7) always has Eb=(be-d1τ1/a,0) as an boundary equilibrium. The components of any interior equilibrium must satisfy (3.1)y=(be-d1τ1-ax)(x2+dx+e)l,me-d2τ2xx2+dx+e=d3.
Whether an interior equilibrium E*=(x*,y*) is feasible or not depends on the values of the parameters. Figure 1 shows the x- and y-isoclines of the system obtained from (1.7) by removing the delays from the arguments. It is clear that a unique interior equilibrium will exist if and only if (3.2)A2+AB+C<0
holds. Here, A=be-d1τ1/a, B=d-(me-d2τ2/d3), C=e,(3.3)x*=-B-B2-4C2,y*=(be-d1τ1-ax*)(x*2+dx*+e)l.
From (3.2) we can easily see that B<0. Hence, the positive equilibrium E* exists for all prey's maturation times τ1 and predator's maturation times τ2 in the interval I=[0,τ1*)×[0,τ2*), where τ1* and τ2* are the maximum values which satisfy the following equation: (3.4)d-me-d2τ2d3<0,be-d1τ1>ax*,x*=[(me-d2τ2/d3)-d-((me-d2τ2/d3)-d)2-4e]2.Thex- and y-isoclines, for various of τ1 and τ2, of the system (1.7) by removing the delays from the arguments.
(a)(b)Figure1 also makes clear how the interior equilibrium E*, if feasible, depends on the prey's maturation time delay τ1, the predator's maturation time delay τ2, and the other parameters. In Figure 1(a), we can easily see that for higher τ2, there is no interior equilibrium. Also, increasing τ1 lowers the x-isocline in Figure 1(b), causing the coincidence of E* with ((be-d1τ1/a)/a,0) at a finite value of τ1. For enough higher τ1 there is only boundary equilibrium (b/a,0) and no interior equilibrium. In all our analysis, it will be important to keep track of how the interior equilibrium depends on the parameters.


## 4. Global Stability of the Equilibrium(be-d1τ1/a,0)

The following result gives conditions which are both necessary and sufficient for the global stability of the equilibrium(x,y)=(be-d1τ1/a,0) of system (1.7).Theorem 4.1.
limt→+∞(x(t),y(t))=(be-d1τ1/a,0) holds true if and only if me-d2τ2/d≤d3 holds true.Proof.
For the sufficiency of the theorem, by positivity of solutions,(4.1)ẋ(t)≤be-d1τ1x(t-τ1)-ax2.
This implies that limt→+∞x(t)=be-d1τ1/a≤b/a and therefore there exists Tε>0 and a positive constant N>b/a such that x(t)<N for all t≥Tε. Then, for t≥Tε+τ2,
(4.2)ẏ(t)≤me-d2τ2y(t-τ2)NN2+dN+e-d3y(t)≤me-d2τ2y(t-τ2)d-d3y(t).
By comparison, y(t) is bounded above by the solution u(t) of
(4.3)u′(t)=me-d2τ2u(t-τ2)d-d3u(t),t>Tε+τ2
satisfying u(t)=y(t) for t∈[Tε,Tε+τ2]. Since me-d2τ2/d≤d3, Lemma 2.3 yields that u(t)→0 which proves y(t)→0. Hence by the first equation of system (1.7), we get limt→+∞x(t)=be-d1τ1/a. This proves me-d2τ2/d≤d3 is the sufficient condition for limt→+∞(x(t),y(t))=(be-d1τ1/a,0).
Now, we provelimt→+∞(x(t),y(t))=(be-d1τ1/a,0)⇒me-d2τ2/d≤d3. Assume the contrary, that is, me-d2τ2/d>d3, then system (1.7) has a positive equilibrium (x*,y*), contradicting limt→+∞(x(t),y(t))=(be-d1τ1/a,0) for all solution (x(t),y(t)). Hence, there must be me-d2τ2/d≤d3, and this proves Theorem 4.1.Theorem 4.2.
The system (1.7) is permanent if and only if it satisfies (3.2).To prove Theorem4.2, we engage the persistence theory by Hale and Waltman [28] for infinite dimensional systems; we also refer to Thieme [29]. Now, we present the persistence theory [28] as follows.Consider a metric spaceX with metric d. T is a continuous semiflow on X, that is, a continuous mapping T:[0,∞)×X→X with the following properties:(4.4)Tt∘Ts=Tt+s,t,s≥0,T0(x)=x,x∈X.
Here 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 (4.5)d(x,y)=inf⁡y∈Y⁡d(x,y).
Recall that the positive orbit γ+(x) through x is defined as γ+(x)=⋃t≥0{T(t)x}, and its ω-limit set is ω(x)=⋂τ≥0CL⋃t≥τ{T(t)x}, where CL means closure. Define Ws(A), the stable set of a compact invariant set A as(4.6)Ws(A)={x:x∈X,ω(x)≠ϕ,ω(x)⊂A};
define Aα̃ the particular invariant sets of interest as
(4.7)Aα̃=⋃x∈Aαω(x).(H0)
AssumeX is the closure of open set X0; ∂X0 is nonempty and is the boundary of X0. Moreover, the C0-semigroup T(t) on X satisfies
(4.8)T(t):X0⟶X0,T(t):∂X0⟶∂X0.Lemma 4.3 (see [28, Theorem  4.1, page 392]).
SupposeT(t) satisfies (H0) and (i)
there is at0≥0 such that T(t) is compact for t>t0,(ii)
T(t) is point dissipative in X,(iii)
Aα̃ is isolated and has an acyclic covering M.Then, T(t) is uniformly persistent if and only if for each Mi∈M, Ws(Mi)∩X0=∅.Proof of Theorem4.2.
Claim 1.
The condition (3.2) leads to the permanence of system (1.7).
LetC+([-τ,0],R+2) denote the space of continuous functions mapping [-τ,0] into R+2. We choose
(4.9)C1={(φ0,φ1)∈C+([-τ,0],R+2):φ0(θ)≡0,φ1(θ)>0,θ∈[-τ,0]},C2={(φ0,φ1)∈C+([-τ,0],R+2):φ0(θ)>0,φ1(θ)≡0,θ∈[-τ,0]}.
Denote C=C1∪C2, X=C+([-τ,0],R+2), and X0=IntC+([-τ,0],R+2), then C=∂X0. It is easy to see that system (1.7) possesses two constant solutions in C=∂X0:E0̃∈C1, E1̃∈C2 with
(4.10)E0̃={(φ0,φ1)∈C+([-τ,0],R+2):φ0(θ)≡φ1(θ)≡0,θ∈[-τ,0]},E1̃={(φ0,φ1)∈C+([-τ,0],R+2):φ0(θ)≡be-d1τ1a,φ1(θ)≡0,θ∈[-τ,0]}.
We verify below that the conditions of Lemma 4.3 are satisfied. By the definition of X0 and ∂X0 and system (1.7), it is easy to see that conditions (i) and (ii) of Lemma 4.3 are satisfied and that X0 and ∂X0 are invariant. Hence, (H0) is also satisfied.
Consider condition (iii) of Lemma4.3. We have
(4.11)ẋ(t)|(φ0,φ1)∈C1≡0,
thus x(t)|(φ0,φ1)∈C1≡0 for all t≥0. Hence, we have
(4.12)ẏ(t)|(φ0,φ1)∈C1=-d3y≤0,
from which follows that all points in C1 approach E0̃, that is, C1=Ws(E0̃). Similarly, we can prove that all points in C2 approach E1̃, that is, C2=Ws(E1̃). Hence, Aα̃=E0̃∪E1̃ and clearly it is isolated. Noting that C1∩C2=∅, it follows from these structural features that the flow in Aα̃ is acyclic, satisfying condition (iii) of Lemma 4.3.
Now, we show thatWs(Eĩ)∩X0=∅, i=0,1. By Lemma 4.3, we have x(t),y(t)>0 for all t>0. Assume Ws(E0̃)∩X0≠∅, that is, there exists a positive solution (x(t),y(t)) with limt→+∞(x(t),y(t))=(0,0), then using the first equation of (1.7), we get
(4.13)d(lnx(t))dt=be-d1τ1x(t-τ1)x(t)-ax(t)-ly(x2+dx+e)>be-d1τ12,
for all sufficiently large t. Hence, we have
(4.14)lim⁡t→+∞x(t)=+∞,
contradicting limt→+∞x(t)=0; this proves Ws(E0̃)∩X0=∅.
Now, we verifyWs(E1̃)∩X0=∅; assume the contrary, that is, Ws(E0̃)∩X0≠∅. Then, there exists a positive solution (x(t),y(t)) to system (1.7) with
(4.15)lim⁡t→+∞(x(t),y(t))=(be-d1τ1a,0),
and for sufficiently small positive constant ε with
(4.16)ε<2A+B+B2-4C2,
there exists a positive constant T=T(ε) such that
(4.17)x(t)>be-d1τ1a-ε>0,y(t)<ε,∀t≥T.
By the second equation of (1.7), we have
(4.18)y′(t)>me-d2τ2(be-d1τ1a-ε)y(t-τ2)(be-d1τ1a-ε)2+d(be-d1τ1a-ε)+e-d3y(t),t≥T+τ.
Consider the equation
(4.19)v′(t)=me-d2τ2(be-d1τ1a-ε)v(t-τ2)(be-d1τ1a-ε)2+d(be-d1τ1a-ε)+e-d3v(t),t≥T+τ,v(t)=y(t),t∈[T,T+τ].
By (4.18) and the comparison theorem, we have y(t)≥v(t) for all t>T. On the other hand, using Theorem  4.9.1 of [16, page 159], we have limt→+∞y(t)>ε, contradicting y(t)<ε as t≥T. Thus we have Ws(Eĩ)∩X0=∅, i=0,1. Now, we get that system (1.7) satisfies all conditions of Lemma 4.3, thus (x(t),y(t)) is uniformly persistent, that is, there exists positive constants ϵ and T=T(ϵ) such that (x(t),y(t))≥ϵ for all t≥T; noting Lemma 2.4 shows that (x,y) are ultimately bounded, and this proves the permanence of system (1.7). This completes the proof.


## 5. Linearised Analysis

We shall concentrate on the dynamics analysis in the two cases whend2>0 and d2=0. These cases correspond, respectively, to the presence or absence of mortality among the immature predators.We will begin by examining the linear stability of the equilibrium(x*,y*), assuming of course that (3.2) holds, so that the equilibrium is feasible. If d2>0 (i.e., there is mortality among the immature predators) then, as one increases the delay τ2, the equilibrium loses feasibility at a finite value of τ2. If d2=0 then, as prey's maturation time delay τ1 increases, the equilibrium loses feasibility at a finite value of τ1.Let us linearise (1.7) at the interior equilibrium (x*,y*). Setting x=x*+u, y=y*+v, where u and v are small, and linearising, gives (5.1)du(t)dt=(-2ax*-lpx*′)u(t)-lpy*′v(t)+be-d1τ1u(t-τ1),dv(t)dt=-d3v(t)+me-d2τ2px*′u(t-τ2)+d3v(t-τ2),
where p(x,y)=xy/(x2+dx+e), px′=((e-x2)y)/((x2+dx+e)2), py′=x/(x2+dx+e).The characteristic equation atE* is as follows: (5.2)G(λ,τ1,τ2)=λ2-be-λτ1(λe-d1τ1+d3e-d1τ1)-d3e-λτ2(λ+2ax*)+bd3e-λτ1-λτ2-d1τ1+λ(d3+2ax*+lpx*′)+d3(2ax*+lpx*′)=0.
Notice that G(0,τ1,τ2)=ld3px*′>0 if x*<e. We always assume x*<e below. So λ=0 is not a solution of the characteristic equation (5.2). Thus, if there is any stability switch of the trivial solution of the linearized system (5.1), there must exist a pair of pure conjugate imaginary roots of the characteristic equation (5.2). When τ1=τ2=0, the original model (1.7) is an ODE model. The characteristic equation of its linearized equation is given by(5.3)G(λ,0,0)=λ2+[2ax*-b]λ+ld3px*′=0.
Clearly, all roots of (5.2) with τ1=τ2=0 have negative real parts provided (5.4)2ax*-b>0.
Hence, we get the following conclusion.Theorem 5.1 (Kar and Pahari [30, Theorem  2.1]).
Assume (3.2) and (5.4) hold. Then the positive equilibrium E*=(x*,y*) of (1.7) is globally asymptotically stable if τ1=τ2=0.In the remaining part of this section, we assume thatτ1 is a small delay, that is, let e-λτ1-d1τ1=1-λτ1-d1τ1. So, (5.2) takes the general form (5.5)P(λ,τ1,τ2)+Q(λ,τ1,τ2)e-λτ2=0,
where(5.6)P(λ,τ1,τ2)=(1+bτ1)λ2+A(τ1,τ2)λ+B(τ1,τ2),Q(λ,τ1,τ2)=C(τ1,τ2)λ+D(τ1,τ2).
Here,(5.7)A(τ1,τ2)=d3-b+lpx*′+2ax*+bd1τ1+bd3τ1,B(τ1,τ2)=bd3(d1τ1-1)+d3(lpx*′+2ax*),C(τ1,τ2)=-d3(1+bτ1),D(τ1,τ2)=-2ad3x*+bd3(1-d1τ1).Let us first consider the caseτ1>0 and τ2=0, then the characteristic equation (5.5) is (5.8)(1+bτ1)λ2+[A(τ1,0)+C(τ1,0)]λ+[B(τ1,0)+D(τ1,0)]=0.
Clearly, all roots of (5.5) with τ1>0 and τ2=0 have negative real parts provided(H1)
A(τ1,0)+C(τ1,0)>0 and B(τ1,0)+D(τ1,0)>0.We summarize the above analysis in the following theorem for model (1.7).Theorem 5.2.
Assume that (3.2), (5.4), and (H1) hold. Then, the positive equilibrium E* of (1.7) is locally asymptotically stable if τ1*>τ1>0 and τ2=0. (Figure 2).a=1.4, b=0.4, d=0.5, e=0.05, and l=m=d1=d3=0.5. E* is asymptotically stable for τ1=0.08 and τ2=0. (a) shows the trajectories graphs of the system (1.7) with initial conditions near E*. (b) shows the phase portrait of system (1.7) corresponding to (a).
(a)(b)For now, fixτ1>0, assume τ2>0, and regard τ2 as a bifurcation parameter to obtain finer results on the stability of E*. Note that (5.2) takes the form of a second-degree exponential polynomial in λ, with all the coefficients of P and Q depending on τ1 and τ2. Thus, we use the method introduced by Beretta and Kuang [31], which gives the existence of purely imaginary roots of a characteristic equation with delay-dependent coefficients (see also [32, 33]). In order to apply the criterion in [31], we need to verify the following properties for ω>0 and τ2∈I with I defined in (5.6). For simplicity, we drop the dependence of τ1 and always assume that τ2∈I:(i)
P(0,τ2)+Q(0,τ2)≠0,(ii)
P(iω,τ2)+Q(iω,τ2)≠0,(iii)
limsup|λ|→∞,Reλ≥0|Q(λ,τ2)/P(λ,τ2)|<1,(iv)
F(ω)=|P(iω,τ2)|2-|Q(iω,τ2)|2 has a finite number of zeros,(v)
each positive rootω(τ2) of F(ω(τ2))=0 is continuous and differentiable in τ2 whenever it exists.Here,P(λ,τ2) and Q(λ,τ2) are defined by (5.5) and (5.9)I={τ2≥0:∃ω(τ2)>0suchthatF(ω(τ2))=0}.LetF be defined as in (iv). By (5.5), we obtain (5.10)F(ω)=(1+bτ1)2ω4+a1(τ2)ω2+a2(τ2)
with (5.11)a1(τ2)=A2(τ2)-2B(τ2)(1+bτ1)-C2(τ2),a2(τ2)=B2(τ2)-D2(τ2).Before proceeding further, let us analyze the structure ofI described by (5.9). Set (5.12)z±(τ2)=-a1(τ2)±a12(τ2)-4a2(τ2)(1+bτ1)22(1+bτ1)2.We need the following hypotheses.(H2)
A(τ2)+C(τ2)≠0 and B(τ2)+D(τ2)≠0.(H3)
0<a2(τ2)<a12(τ2)/4(1+bτ1)2 and a1(τ2)<0.(H4)
a2(τ2)≤0.Theorem 5.3.
Assume (3.2) holds.(i)
If (H3) is satisfied, thenF(ω(τ2))=0 has two different positive roots z±(τ2) denoted by ω±, respectively.(ii)
If (H4) is met, thenF(ω(τ2))=0 has a unique positive root ω+=z+(τ2).DefineI1={τ2≥0: (3.2), (H2) and (H3) hold},I2={τ2≥0: (3.2), (H2) and (H4) hold}.Then,(5.13)I=I1⋃I2.
Now, we verify the properties (i)–(v) for τ2∈I. Firstly, (i) and (ii) are satisfied. Indeed, (5.14)P(0,τ2)+Q(0,τ2)=B(τ2)+D(τ2)≠0,P(iω,τ2)+Q(iω,τ2)=[B(τ2)+D(τ2)-ω2]+iω[A(τ2)+C(τ2)]≠0.
Then, from (5.5), we know (5.15)lim⁡|λ|→∞|Q(λ,τ2)P(λ,τ2)|=0.
Therefore, (iii) follows. Finally, (5.10) implies (iv) and (v).Letλ=iω(ω>0) be a root of (5.2). Substituting it into (5.2) and separating the real and imaginary parts yield (5.16)(1+bτ1)ω2-B(τ2)=D(τ2)cos⁡ωτ2+C(τ2)ωsinωτ2,A(τ2)ω=D(τ2)sinωτ2-C(τ2)ωcos⁡ωτ2.
Then, we have(5.17)sinωτ2=((1+bτ1)ω2-B(τ2))Cω+ωADω2C2+D2,(5.18)cos⁡ωτ2=((1+bτ1)ω2-B(τ2))D-ω2ACω2C2+D2.
Here, the dependence of the coefficients on τ2 is implicitly assumed.By the definitions ofP(λ,τ1,τ2) and Q(λ,τ1,τ2) in (5.6) and applying the property (ii), (5.18) can be written as (5.19)sinωτ2=Im⁡(P(iω,τ1,τ2)Q(iω,τ1,τ2)),cos⁡ωτ2=-Re(P(iω,τ1,τ2)Q(iω,τ1,τ2)),
which lead us to (5.20)F(ω)=0.
Noting that Theorem 5.3 has given the explicit expressions of ω(τ2) that satisfies F(ω(τ2))=0 on τ2∈I, define θ(τ2)∈[0,2π] for τ2∈I(5.21)sinθ(τ2)=((1+bτ1)ω2-B(τ2))Cω+ωADω2C2+D2,(5.22)cos⁡θ(τ2)=((1+bτ1)ω2-B(τ2))D-ω2ACω2C2+D2.
Such θ(τ2) defined above is well and uniquely defined for all τ2∈I (see [31]).Therefore,iω* with ω*=ω(τ2*)>0 is a purely imaginary root of (5.2) if and only if τ2* is a zero of the Sn, where (5.23)Sn(τ2)=τ2-θ(τ2)+2πnw(τ2),τ2∈I,
with n∈N0={0,1,2,…}.Hence, once we know suchτ2*, this will give us a pair of delay values (τ1,τ2*) at which the stability switch may be possible when increasing the value of τ2 while keeping τ1 fixed (Figure 3). The following result is due to Beretta and Kuang [31].A solution of model (1.7) with initial conditions (0.16, 0.05), where a=1.4, b=1.1, d=0.5, e=0.05, l=m=d1=d3=0.5, d2=0 (i.e., no mortality among the immature predators), τ1=1.5, and τ2 varies from 0.1 to 20.
(a)(b)(c)(d)Lemma 5.4.
The characteristic equation (5.2) has a pair of simple pure imaginary roots λ=±iω(τ2*) at τ2*∈I, provided Sn(τ2*)=0 for n∈N0. Moreover, if ω(τ2*)=ω+(τ2*), then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if δ+(τ2*)>0 and crosses the imaginary axis from right to left if δ+(τ2*)<0, where
(5.24)δ+(τ2*)=sign⁡{dRe(λ)dτ2|λ=iω+(τ2*)}=sign⁡{dSn(τ2)dτ2|τ2=τ2*},
and if ω(τ2*)=ω-(τ2*), then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if δ-(τ2*)>0 and crosses the imaginary axis from right to left if δ-(τ2*)<0, where
(5.25)δ-(τ2*)=sign⁡{dRe(λ)dτ2|λ=iω-(τ2*)}=-sign⁡{dSn(τ2)dτ2|τ2=τ2*}.IfI=I2, then only ω+ is feasible. In this case, we can easily observe that Sn(0)≤0 and Sn(τ2)>Sn+1(τ2) for τ2∈I and n∈N0, then without loss of generality, we may suppose (5.26)dSn(τ2nj)dτ2≠0withSn(τ2nj)=0.
Then, for one thing, stability switches occur at the zeros of S0(τ2), denoted by τ20j, if (H1) holds. For another, applying Theorem 5.2 and Hopf bifurcation theorem for functional differential equations (see Hale's book [34]), we can conclude the existence of a Hopf bifurcation. Before stating the main theorem, denote(5.27)τ2m=min⁡{τ2∈I:S0(τ2)=0},τ2M=max⁡{τ2∈I:S0(τ2)=0}.Theorem 5.5.
Assume (3.2), (H1), (H2), and (H4) hold.(i)
IfS0(τ2) has no positive zeros in I, then E* is locally asymptotically stable for all τ2≥0.(ii)
IfSn(τ2) has at least one positive zero in I and (5.26) is met, then E* is locally asymptotically stable for τ2∈[0,τ2m)∪(τ2M,∞) and unstable for τ2∈(τ2m,τ2M), that is, stability switches of stability-instability-stability occur. Hopf bifurcation takes place when τ2=τ2nj, n∈N0.IfI=I1, then both ω+ and ω- are feasible for τ2∈I. We have the following two sequences of functions on I:
(5.28)Sn+(τ2)=τ2-θ+(τ2)+2nπω+(τ2),Sn-(τ2)=τ2-θ-(τ2)+2nπω-(τ2),
where θ±(τ2) is the solution of (5.22) when ω=ω±. Similarly, it is obtained for τ2∈I that Sn±(0)≤0 and Sn±(τ2)>Sn+1±(τ2) with n∈N0. Furthermore, if S0+(τ2)>S0-(τ2), then Sn+(τ2)>Sn-(τ2) for n∈N0. Other than the case above, stability switch may depend on all real roots of Sn+(τ2)=0 and Sn-(τ2)=0. In addition, Hopf bifurcations can also appear at each root of Sn±(τ2)=0. Naturally, we can also obtain results similar to Theorem 5.2. Here, we omit the corresponding statements.IfI=I1∪I2 with I1≠∅ and I2≠∅, then we can discuss different cases on their own sets, respectively. The remaining discussion is the same as above.


## 6. Discussion

In this paper, we study the predator-prey model (1.3) of Holling type-IV functional response with stage structure on prey and predator which we consider a fairly realistic one in this category.By Theorem4.1, We have that limt→∞(x(t),y(t))=(be-d1τ1/a,0) if and only if the following condition holds true: me-d2τ2/d≤d3.  It is implied that the boundary equilibrium Eb of system (1.7) is globally stable.Our main purpose of this paper is to analyze a two-species predator-prey autonomous model with stage structure for both prey and predator, in which there are two discrete delays due to the maturity for both immature prey and immature predator, respectively. Under certain initial conditions, the boundary and the existence of the coexistence equilibrium of system (1.7) were investigated, and also the stability switches of the coexistence equilibrium of system (1.7) were shown by analyzing the corresponding characteristic equation as predator's maturation time delay (through-stage time delay). τ2 is increased from zero. Additionally, in the last section, we have acquired very rich dynamical behaviours of the nontrivial equilibrium point E* when varying the value of τ1(τ2) while keeping τ2(τ1) fixed. Particularly, when we let τ2=0 (see Figure 2), our simulations show that the system does not have sustained  oscillations when the other delay parameter τ1 is not large. In the special case when keeping τ1>0 fixed and τ2 varied, the oscillatory dynamics will  persist and gain complexity when we increase the delay τ2 (see Figure 3). Such distinct dynamical outcomes highlight the importance of incorporating the through-stage death rate in stage structured population models.Observable delay effects are often gradual (distributed) and smooth in most dynamical systems, it is thus natural to utilize distributed delay parameters rather than discrete delays when modeling these systems. In other words, discrete delay is often a simplification of the complicated dynamical process that is almost always best represented by smooth (continuous and distributed) delay. However, mathematically, a single discrete delay alone can often generate rich dynamics that enable and facilitate nontrivial and interesting biological observations as evidenced by this work. Nevertheless, we plan to pursue additional studies on the predator-prey system and sustained oscillations through models with distributed time delays.



---

*Source: 509871-2011-05-10.xml*