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

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.


Introduction
Predator-prey models are arguably the most fundamental building blocks of the any bioand 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 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 x 2 , 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 where x i t and x m 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 −γτ x m 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: where τ max{τ 1 , τ 2 }, x θ , y θ > 0 is continuous on −τ ≤ θ ≤ 0, and x j 0 , x 0 , y j 0 > 0, y 0 > 0, x and y represent mature prey and the predator densities, respectively.And x j and y j 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 d 1 and take τ 1 units of time to mature, thus e −d 1 τ 1 is the surviving rate of each immature prey to reach maturity.The constant d 2 is the death rate of the juvenile predator and take τ 2 units of time to mature, thus e −d 2 τ 2 is the surviving rate of each immature predator to reach maturity.d 3 is the death rate of the mature predator.The mature predator consumes the mature prey with functional response of Holling type-IV lx/ x 2 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 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 that is, x j t is completely determined by x t .Thus, the following system can be separated from system 1.3 where x θ , y j θ , y θ ≥ 0 are continuous on −τ ≤ θ ≤ 0, and x 0 , y j 0 , y 0 > 0.
Notice that, mathematically, no information on the past history of y j 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 In our study, we assume that the initial conditions of system 1.7 take the form the Banach space of continuous functions mapping the interval −τ, 0 into R 2 0 , where R 2 0 { x, y : x ≥ 0, y ≥ 0}.
Remark 1.1.Because x j t is completely determined by x t and y j 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.

Preliminary Analysis
To prove the main results, we need the following lemmas.In the biology significance, we only study in Using the similar arguments to Lemma 1 in 27 , we directly have.
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 By Lemma 2.2, a comparison argument shows that lim sup It implies that x t is ultimately bounded.No less of generality, we suppose that there exists By Lemma 2.2 and the comparison theorem, we get which implies that there exists a constant C > be −d 1 τ 1 /a > 0, such that all trajectories initiating in R 2 enter the region

Equilibria and Their Feasibility
Apart from the zero solution, system 1.7 always has E b be −d 1 τ 1 /a, 0 as an boundary equilibrium.The components of any interior equilibrium must satisfy

3.1
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.3
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 where τ * 1 and τ * 2 are the maximum values which satisfy the following equation:

3.4
Figure 1 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 −d 1 τ 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.

Global Stability of the Equilibrium be
The following result gives conditions which are both necessary and sufficient for the global stability of the equilibrium x, y be −d 1 τ 1 /a, 0 of system 1.7 .Proof.For the sufficiency of the theorem, by positivity of solutions, 4.1 This implies that lim t → ∞ x t be −d 1 τ 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 By comparison, y t is bounded above by the solution u t of 3 yields that u t → 0 which proves y t → 0. Hence by the first equation of system 1.7 , we get lim t → ∞ x t be −d 1 τ 1 /a.This proves me Assume the contrary, that is, me −d 2 τ 2 /d > d 3 , then system 1.7 has a positive equilibrium x * , y * , contradicting lim t → ∞ x t , y t be −d 1 τ 1 /a, 0 for all solution x t , y t .Hence, there must be me −d 2 τ 2 /d ≤ d 3 , 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 Theorem 4.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 space X with metric d.T is a continuous semiflow on X, that is, a continuous mapping T : 0, ∞ × X → X with the following properties: Here T t denotes the mapping from X to X given by T t x T t, x .The distance d x, y of a point x ∈ X from a subset Y of X is defined by Recall that the positive orbit γ x through x is defined as γ x t≥0 {T t x}, and its ωlimit set is ω x τ≥0 CL t≥τ {T t x}, where CL means closure.Define W s A , the stable set of a compact invariant set A as define A α the particular invariant sets of interest as H 0 Assume X is the closure of open set X 0 ; ∂X 0 is nonempty and is the boundary of X 0 .Moreover, the C 0 -semigroup T t on X satisfies iii A α is isolated and has an acyclic covering M.

Then, T t is uniformly persistent if and only if for each
Proof of Theorem 4.2.
Claim 1.The condition 3.2 leads to the permanence of system 1.7 .Let C −τ, 0 , R 2 denote the space of continuous functions mapping −τ, 0 into R 2 .We choose It is easy to see that system 1.7 possesses two constant solutions in C ∂X 0 :

4.10
We verify below that the conditions of Lemma 4.3 are satisfied.By the definition of X 0 and ∂X 0 and system 1.7 , it is easy to see that conditions i and ii of Lemma 4.3 are satisfied and that X 0 and ∂X 0 are invariant.Hence, H 0 is also satisfied.Consider condition iii of Lemma 4.3.We have ẋ t | ϕ 0 ,ϕ 1 ∈C 1 ≡ 0, 4.11 thus x t | ϕ 0 ,ϕ 1 ∈C 1 ≡ 0 for all t ≥ 0. Hence, we have ẏ t ϕ 0 ,ϕ 1 ∈C 1 −d 3 y ≤ 0, 4.12 from which follows that all points in C 1 approach E 0 , that is, C 1 W s E 0 .Similarly, we can prove that all points in C 2 approach E 1 , that is, C 2 W s E 1 .Hence, A α E 0 ∪ E 1 and clearly it is isolated.Noting that C 1 ∩ C 2 ∅, it follows from these structural features that the flow in A α is acyclic, satisfying condition iii of Lemma 4.3.Now, we show that W s E i ∩ X 0 ∅, i 0, 1.By Lemma 4.3, we have x t , y t > 0 for all t > 0. Assume W s E 0 ∩ X 0 / ∅, that is, there exists a positive solution x t , y t with lim t → ∞ x t , y t 0, 0 , then using the first equation of 1.7 , we get for all sufficiently large t.Hence, we have lim Now, we verify W s E 1 ∩ X 0 ∅; assume the contrary, that is, W s E 0 ∩ X 0 / ∅.Then, there exists a positive solution x t , y t to system 1.7 with lim t → ∞ x t , y t be −d 1 τ 1 a , 0 , 4.15 and for sufficiently small positive constant ε with there exists a positive constant T T ε such that By the second equation of 1.7 , we have Consider the equation

4.19
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 lim t → ∞ y t > ε, contradicting y t < ε as t ≥ T .Thus we have W s E i ∩ X 0 ∅, 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.

Linearised Analysis
We shall concentrate on the dynamics analysis in the two cases when d 2 > 0 and d 2 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 d 2 > 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 d 2 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
where p x, y xy/ x 2 dx e , p x e − x 2 y / x 2 dx e 2 , p y x/ x 2 dx e .
The characteristic equation at E * is as follows:

5.2
Notice that G 0, τ 1 , τ 2 ld 3 p x * > 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 Clearly, all roots of 5.2 with τ 1 τ 2 0 have negative real parts provided 2ax * − b > 0.

5.4
Hence, we get the following conclusion.In the remaining part of this section, we assume that τ 1 is a small delay, that is, let e −λτ 1 −d 1 τ 1 1 − λτ 1 − d 1 τ 1 .So, 5.2 takes the general form where

5.7
Let us first consider the case τ 1 > 0 and τ 2 0, then the characteristic equation 5.5 is Clearly, all roots of 5.5 with τ 1 > 0 and τ 2 0 have negative real parts provided We summarize the above analysis in the following theorem for model 1.7 .
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: 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 .
Here, P λ, τ 2 and Q λ, τ 2 are defined by 5.5 and Let F be defined as in iv .By 5.5 , we obtain with

5.11
Before proceeding further, let us analyze the structure of I described by 5.9 .Set

5.12
We need the following hypotheses.

5.13
Now, we verify the properties i -v for τ 2 ∈ I. Firstly, i and ii are satisfied.Indeed,

5.15
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
Then, we have Here, the dependence of the coefficients on τ 2 is implicitly assumed.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 .If I I 1 , then both ω and ω − are feasible for τ 2 ∈ I.We have the following two sequences of functions on I: where θ ± τ 2 is the solution of 5.22 when ω ω ± .Similarly, it is obtained for τ 2 ∈ I that S ± n 0 ≤ 0 and n τ 2 for n ∈ N 0 .Other than the case above, stability switch may depend on all real roots of S n τ 2 0 and S − n τ 2 0. In addition, Hopf bifurcations can also appear at each root of S ± n τ 2 0. Naturally, we can also obtain results similar to Theorem 5.2.Here, we omit the corresponding statements.
If I I 1 ∪ I 2 with I 1 / ∅ and I 2 / ∅, then we can discuss different cases on their own sets, respectively.The remaining discussion is the same as above.

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 Theorem 4.1, We have that lim t → ∞ x t , y t be −d 1 τ 1 /a, 0 if and only if the following condition holds true: me −d 2 τ 2 /d ≤ d 3 .It is implied that the boundary equilibrium E b 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.

Figure 1 :
Figure 1: The x-and y-isoclines, for various of τ 1 and τ 2 , of the system 1.7 by removing the delays from the arguments.

Figure 2 :
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 .

8
Lemma 4.3 see 28, Theorem 4.1, page 392 .Suppose T t satisfies H 0 and i there is a t 0 By the definitions of P λ, τ 1 , τ 2 and Q λ, τ 1 , τ 2 in 5.6 and applying the property ii , 5.18 can be written as 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 imaginary axis from right to left if δ τ * If I I 2 , then only ω is feasible.In this case, we can easily observe that S n 0 ≤ 0 and S n τ 2 > S n 1 τ 2 for τ 2 ∈ I and n ∈ N 0 , then without loss of generality, we may suppose 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 i If S 0 τ 2 has no positive zeros in I, then E * is locally asymptotically stable for all τ 2 ≥ 0.ii If S n τ 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 τ j 2n , n ∈ N 0 .