# Stability Analysis of a Harvested Prey-Predator Model with Stage Structure and Maturation Delay

Authors: Chao Liu; Qingling Zhang; James Huang
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

A harvested prey-predator model with density-dependent maturation delay and stage structure for prey is proposed, where selective harvest effort on predator population is considered. Conditions which influence positiveness and boundedness of solutions of model system are analytically investigated. Criteria for existence of all equilibria and uniqueness of positive equilibrium are also studied. In order to discuss effects of maturation delay and harvesting on model dynamics, local stability analysis around all equilibria of the proposed model system is discussed due to variation of maturation delay and harvest effort level. Furthermore, global stability of positive equilibrium is investigated by utilizing an iterative technique. Finally, numerical simulations are carried out to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

In the natural world, many species have a life history that takes them through two stages, juvenile stage and adult stage. Individuals in each stage are identical in biological characteristics, and some vital rates (rates of survival, development, and reproduction) of individuals in a population almost always depend on stage structure. Furthermore, many complex biological phenomena arising in prey-predator ecosystem always depend on the past history of system, and it has been recognized that time delay may have complicated impact on dynamics of prey-predator ecosystem [1]. In the past several decades, there has been an increasing interest in prey-predator model system with stage structure and time delay (see [2–26] and the references therein).In the model proposed by Aiello and Freedman [2], stage structure of single population growth with stage structure and time delay representing for maturation of population is considered. Their model predicts a positive steady state as the global attractor, thereby suggesting that stage structure does not generate sustained oscillations frequently observed in single population in the real world. Subsequent work made by other authors [3, 6, 7, 12–14] suggests that time delay to adulthood should be state dependent. Generally, boundedness and persistence of solutions of model system may be affected by introduction of time delay into prey-predator system with stage structure [14, 15, 20–22, 24–26]. Time delay can also cause loss of stability and other complicated dynamical behavior [27]. Especially, there is a well-developed theory of stage-structured models which incorporate time delay into maturity of population [4].It is well known that harvesting has a strong impact on dynamic evolution of a population; there has been considerable interest in the modeling of harvesting of biological resources [1]. In these models, the harvesting effort is considered to be a dynamic variable; several kinds of harvesting policies are utilized to study the dynamical behavior of the model system. In recent years, there has been growing interest in the study of stage-structured prey-predator system with harvesting. Several prey-predator models with stage structure and harvest effort on predator have been investigated in [28–33] and the references therein.Recently, Huo et al. [24] investigated dynamical behavior and stability of the following stage-structured system with time delay:
(1)x˙1(t)=r1x2(t)-dx1(t)-r1e-dτx2(t-τ),x˙2(t)=r1e-dτx2(t-τ)-bx22(t)-a1y(t)x2(t)x2(t)+k1,y˙(t)=y(t)(r2-a2y(t)x2(t)+k2),
where x1(t), x2(t), and y(t) represent the density of immature prey population, mature prey population and predator population, at time t, respectively; r1 is the intrinsic growth rate of mature prey population, and d is the death rate of immature prey population. Constant τ≥0 denotes maturation delay of immature prey population to mature prey population, and the term r1e-dτx2(t-τ) represents the immature prey population who were born at time t-τ and survived at time t. b denotes the intracompetition rate for mature prey population due to overcrowding phenomenon with mature prey population. a1 is the maximum value of the per capita reduction rate of mature prey population due to predator population, and a2 is the maximum value of the per capita reduction rate of predator population due to mature prey population. k1 measures the extent to which the environment provides protection to mature prey population, and k2 measures the extent to which the environment provides protection to predator population. r2 represents the maximal per capita growth rate of predator population. All the parameters mentioned previously are all positive constants. Furthermore, global stability of positive equilibrium of model system (1) is investigated in [26].It is well known that the length of time for prey population to maturity is density dependent; that is, maturation time depends on the total population amount of prey population within prey predator ecosystem, and prey population takes less time to reach maturity with depletion of predator population [23, 34–36]. Density-dependent maturity of population in prey predator ecosystem is discussed in their work, which reveals that density-dependent effects of the predators’ counterparts to prey defenses and the density dependence effect of each type of predator offense are analogous to the corresponding type of prey defense. Dynamical behavior and stability switch is investigated in [23, 34–36]. However, harvest effort on population within prey-predator ecosystem is not considered in [23, 34–36].By assuming maturity delay of prey population is density dependent and predator population is harvested; work done in [24] is extended in this paper, and a harvested prey predator model with density-dependent maturation delay and stage structure for prey population is proposed in the second section of this paper. In the third section of this paper, positiveness and boundedness of solution of the proposed model are studied, and the conditions for existence of equilibria and uniqueness of positive equilibrium are also investigated. Local stability analysis around all equilibria is discussed due to variation of maturation delay as well as harvest effort level. Furthermore, global stability of the positive equilibrium of the proposed model system is studied by utilizing an iterative technique. In the fourth section of this paper, numerical simulations are carried out to show consistency with theoretical analysis. Finally, this paper ends with a conclusion.


## 2. Model Formulation

Based on the previous analysis, the model proposed by Huo et al. in [24] is extended by incorporating harvest effort on predator population and assuming that maturation delay of prey population is density dependent, and the model can be governed by the following differential equations:
(2)x˙1(t)=r1x2(t)-dx1(t)-r1e-dτ(z(t))x2(t-τ(z(t))),x˙2(t)=r1e-dτ(z(t))x2(t-τ(z(t)))-bx22(t)-a1y(t)x2(t)x2(t)+k1,y˙(t)=y(t)(r2-a2y(t)x2(t)+k2)-qEy(t).
The initial conditions for model system (2) take the following form:
(3)x1(0)>0,y(0)>0,x2(θ)=ϕ(θ)>0,θ∈[-τ^,0),
where z(t)=x1(t)+x2(t)+y(t), a scalar E≥0 denotes the harvesting effort to predator population, constant q is the catchability coefficient of predator, and the harvesting term qEy(t) follows the catch per unit effort hypothesis [1]. Furthermore, r1, r2, d, b, a1,a2, k1, and k2 in model system (2) share the same interpretations mentioned in model system (1).In the following section of this paper, model system (2) is derived under the following hypotheses. (H1)
Prey population is divided into two-stage groups, that is, immature and mature. The termr1e-dτ(z(t))x2(t-τ(z(t))) represents the immature prey population born at time t-τ(z(t)) and survive at time t with death rate d, which represents transformation term from immature prey to mature prey.(H2)
Density-dependent time delayτ(z(t)) is taken to be an increasing differentiable bounded function of the total population (immature prey, mature prey, and predator population), which satisfies
(4)d[τ(z(t))]dt≥0,0≤τ0≤τ(z(t))≤τ1,limz(t)→0+τ(z(t))=τ0,limz(t)→+∞τ(z(t))=τ1.(H3)
For the continuity of initial conditions, it is required that(5)x1(0)=∫-τ^0‍r1edsϕ(s)ds,
where ϕ(s) is assumed to be continuous function (for mathematical reason) and nonnegative (for biological reason).(H4)
In order to exclude the possibility of immature prey becoming mature prey except by birth,t-τ(z(t)) is assumed to be a strictly increasing function of t. Otherwise, there are two different times at which the same individual immature prey turns to be mature prey twice at the same instant of time, which is absurd to practical biological interpretations. (For detailed methodology, see [3].)


## 3.  Qualitative Analysis of Model System

In this section, positiveness and boundedness of solution of model system (2) are analytically investigated. Criteria for existence of equilibria and uniqueness of positive equilibrium are also studied. By using differential dynamical system theory and stability theory, local stability analysis around all equilibria of model system is discussed. Furthermore, global stability of the positive equilibrium of the proposed model system is studied by utilizing an iterative technique.


### 3.1. Positiveness and Boundness of Solutions

Theorem 1.
Under hypotheses (H1)–(H4), solutions of model system (2) with given initial conditions are positive for all t>0.Proof.
Assume that there existst0=inf{t>0∣x2(t)=0}. Based on the continuity t0>0, it can be computed by evaluating the model system (2) at time t0:
(6)x˙2(t0)={r1e-dτ(z(t0))ϕ(t0-τ(z(t0))),0≤t0≤τ^,r1e-dτ(z(t0))x2(t0-τ(z(t0))),t0>τ^.
According to (6) and the initial conditions of model system (2), it is easy to show that x˙2(t0)>0. On the other hand, it follows from the definition of t0 that x˙2(t0)=0, which is a contradiction. Consequently, x2(t)>0 for all t>0.
Based on the positiveness ofx2(t) and the third equation of model system (2), it is easy to show that y(t)>0 for y(0)>0,t>0.
Consider the equation(7)u˙(t)=-du(t)-r1e-dτ(z)x2(t-τ(z)),u(0)=x1(0)>0.
It is obvious to show thatu˙(t)<0; that is, u(t) is strictly decreasing. By virtue of positiveness of x2(t), y(t), it derives that
(8)x1(t)>u(t),0<t≤τ^.
By solving (7), it gives that
(9)u(t)=e-dtu(0)-r1e-dt∫0t‍edqe-dτ(z(q))x2(q-τ(z(q)))dq.
According to (5), it derives that
(10)u(t)=e-dt∫-τ^0‍r1edsϕ(s)ds-r1e-dt∫0t‍edqe-dτ(z(q))x2(q-τ(z(q)))dq.
By substitutingp=q-τ(z(q)) in the above equation, it can be obtained that
(11)u(t)=r1e-dt∫-τ^0‍edsϕ(s)ds-r1e-dt∫-τ^t-τ(z(t))‍edpx2(p)1-τ′(z)z˙(p)dp,
which implies that
(12)u(τ^)=r1e-dτ^∫-τ^0‍edsϕ(s)ds-r1e-dτ^∫-τ^τ^-τ(z(τ^))‍edpx2(p)1-τ′(z)z˙(p)dp.
According tox2(t)>0,y(t)>0 for all t>0 and t-τ(z(t)) is an increasing function based on (H4). 1-τ′(z)z˙(t)>0 holds for -τ^≤t≤τ^-τ(z(τ^)), and the following inequality can be obtained:
(13)u(τ^)≥r1e-dτ^∫-τ^0‍edpϕ(p)(1-11-τ′(z)z˙(p))dp>0.
since u(τ^)>0 and u(t) is strictly decreasing, x1(t)>0, -τ^≤t≤0. By repeating this argument to include all positive time, it can be shown that x1(t)>0 for all t>0. Hence, solutions of model system (2) with given initial conditions are positive for all t>0.Theorem 2.
If the hypotheses (H1)–(H4) hold andr2>qE, all solutions of model system (2) are bounded within a region Ω:
(14)Ω={(x1(t),x2(t),y(t))∣60<x1(t)+x2(t)≤R1,0<y(t)≤R2(x1(t),x2(t),y(t))},
where R1=(r1+d)2/4bd,R2=(r2-qE)(R1+k2)/a2.Proof.
Letv(t)=x1(t)+x2(t), and it is easy to show that v(t)>0 based on positiveness of solutions of model system (2). By calculating the derivative of v(t) along the solutions, it gives that
(15)v˙(t)=r1x2(t)-dx1(t)-bx22(t)-a1x2(t)y(t)x2(t)+k1<(r1+d)x2(t)-bx22-dv(t).
By using the standard comparison principle in (15), it derives that
(16)limt→∞sup(v(t))≤(r1+d)24bd=R1.
It follows from the third equation of model system (2) that
(17)y˙(t)=(r2-qE)y(t)-a2y2(t)x2(t)+k2≤(r2-qE)y(t)-a2y2(t)R1+k2.
By utilizing the standard comparison principle in inequality (17) and r2-qE>0, it gives that
(18)limt→∞supy(t)≤(r2-qE)(R1+k2)a2=R2.
Consequently, all solutions of model system (2) are bounded within a region Ω:
(19)Ω={(x1(t),x2(t),y(t))∣60<x1(t)+x2(t)≤R1,0<y(t)≤R2(x1(t),x2(t),y(t))},
where R1=(r1+d)2/4bd, R2=(r2-qE)(R1+k2)/a2.


### 3.2. Existence of Equilibria and Uniqueness of Positive Equilibrium

The existence of biologically reasonable equilibria of model system (2) is investigated in this subsection. Since the biological interpretation of the positive equilibrium implies that immature prey, mature prey, and predator population all exist, uniqueness of positive equilibrium is also studied.By simple computation, there are two equilibria denoted byP0(0,0,0) and P1(0,0,k2(r2-qE)/a2). The biological interpretations of P0, P1 are as follows. For P0(0,0,0), it implies that all population in harvested prey predator ecosystem does not exist. For P1(0,0,k2(r2-qE)/a2), it implies that there is not any predation source for predator population. It follows from the previous biological interpretations that population in such ecosystem cannot be maintained at an ideal level for sustainable development, which are not relevant to major investigation in this paper.Furthermore, there is one or more positive equilibria denoted byP*(x1*,x2*,y*). In order to discuss the existence of P*, it is equivalent to show that the following equations always have at least one positive solution:
(20)r1x2-dx1-r1e-dτ(x1+x2+y)x2=0,r1e-dτ(x1+x2+y)-bx2-a1yx2+k1=0,r2-qE-a2yx2+k2=0.It follows from (20) that
(21)y=(r2-qE)(x2+k2)a2=f(x2),(22)r1x2-dx1=r1x2e-dτ(x1+g(x2)),(23)r1e-dτ(x1+g(x2))=bx2+a1h(x2),
where g(x2)=x2+y=x2+f(x2) and h(x2)=y/(x2+k1)=f(x2)/(x2+k1).Theorem 3 (existence of positive equilibrium).
Supposing that hypotheses (H1)–(H4) hold, ifk1≥k2, r2>qE, and a2r1e-dτ1>a1r2, then there exists at least one positive equilibrium P*.Proof.
LetΓ1 and Γ2 be the solution curves of (22) and (23) for x1≥0, x2≥0, respectively. The analytical properties of curve Γ1 and Γ2 are as follows.
ForΓ1: by simple computing, it can be found that (0,0)∈Γ1.
According to (H2) and positiveness of all solutions of model system (2), it is easy to show that limx1→+∞τ(x1(t)+g(x2))=τ1, and
(24)limx1→+∞x2(x1)=limx1→+∞dx1r1(1-e-dτ(x1+g(x2))=limx1→+∞dx1r1(1-e-dτ1)=+∞.
ForΓ2: by differentiating x2 against x1 along Γ2, it can be obtained that
(25)dx2dx1=-dr1e-dτ(x1+g(x2))×(a1(r2-qE)(k1-k2)a2(x2+k1)2dr1e-dτ(x1+g(x2))τ′(x1+g(x2))a2+r2-qEa26666+b+a1(r2-qE)(k1-k2)a2(x2+k1)2)-1.
It can be shown that(dx2/dx1)<0, provided that k1≥k2, r2>qE, and then Γ2 is strictly decreasing.
Furthermore, according tok1≥k2, r2>qE, and a2r1e-dτ1>a1(r2-qE),
(26)limx1→∞x2(x1)=1b(r1e-dτ1-a1(r2-qE)(x2+k2)a2(x2+k1))>0.
Based on the above analysis,Γ1 and Γ2 intersect at some positive values, which proves the existence of positive equilibrium P*.Theorem 4 (uniqueness of positive equilibrium).
Supposing that hypotheses (H1)–(H4) hold, if the following inequality holds(27)1-dτ′(x1*+g(x2*))(bx2*+a1h(x2*))66×(a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*))66666×(bx2*+a1h(x2*))a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*)))+a1dh′(x2*)66+dτ′(x1*+g(x2*)g′(x2*)(bx2*+a1h(x2*))66×(d+b+dτ′(x1*+g(x2*)))>0,
then there exists a unique positive equilibrium.Proof.
Based on (22) and (23), x2 can be defined as the function of x1:
(28)Γ1:x2=g1(x1),Γ2:x2=g2(x1).
The positive equilibrium P* will be unique, provided that g1′(x1)>g2′(x1) for every such P* otherwise reverse inequality holds.
By differentiating (22) with respect to x1, it can be obtained that
(29)g1′(x1)=d+r1e-dτ(x1+g(x2))(1-dτ′(x1+g(x2)))r1(1+e-dτ(x1+g(x2))dτ′(x1+g(x2))g′(x2)).
By differentiating (23) with respect to x1, it can be obtained that
(30)g2′(x1)666=-b+dr1τ′(x1+g(x2))e-dτ(x1+g(x2))a1h′(x2)+dr1g′(x2)τ′(x1+g(x2))e-dτ(x1+g(x2)).
On the other hand, some expressions about positive equilibriumP*(x1*,x2*,y*) can be obtained based on (22) and (23),
(31)r1x2*-dx1*=r1x2*e-dτ(x1*+g(x2*)),r1e-dτ(x1*+g(x2*))=bx2*+a1h(x2*).
According to (31), g1′(x1*)>g1′(x2*) is equivalent to the following inequality:
(32)1-dτ′(x1*+g(x2*))(bx2*+a1h(x2*))66×(a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*))666666×(bx2*+a1h(x2*))a1h′(x2*))+a1dh′(x2*)66+dτ′(x1*+g(x2*)g′(x2*)66×(bx2*+a1h(x2*))(d+b+dτ′(x1*+g(x2*)))>0.
This completes the proof.


### 3.3. Local Stability Analysis around Equilibria

Local stability of model system (2) around all equilibria of model system (2) is investigated. Furthermore, stability switch due to variation of maturity delay and harvest effort level is also studied in this subsection.The characteristic equation of model system (2) about some equilibrium P~=(x~1,x~2,y~) takes the following form: (33)|λ+(d-A~)e-dτ(z~)r1(1+e-(λ+d)τ(z~))-A~e-dτ(z~)-A~e-dτ(z~)A~e-dτ(z~)λ+2bx~+a1k1y~(x~2+k1)2+A~e-dτ(z~)+r1e-(λ+d)τ(z~)A~e-dτ(z~)+a1x~2k1+x~20-a2y~2(k2+x~2)2λ-r2+qE+2a2y~k2+x~2|=0,where A~=dr1x~2τ′(z~).Theorem 5.
Local stability analysis of model system (2) around P0 and P1 is as follows: (a)
ifr2<qE, then model system is locally stable around P0, and P1 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane;(b)
ifr2>qE, then model system is locally stable around P1, and P0 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane.Proof.
ForP0(0,0,0), (33) reduces to
(34)(λ+de-dτ(0))(λ+r1e-(λ+d)τ(0))(λ-(r2-qE))=0.
By solving (34), it can be found that there are two negative eigenvalues and only one positive eigenvalue, provided r2>qE, which implies that P0 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane. On the other hand, there are three negative eigenvalues, provided that r2<qE, which implies that P0 is a stable point.
ForP1(0,0,(k2(r2-qE)/a2)), (33) reduces to
(35)(λ+de-dτ(y~))(λ+r1e-(λ+d)τ(y~)+a1y~k1)(λ+(r2-qE))=0,
where y~=k2(r2-qE)/a2. It follows from (35) that there are two negative eigenvalues and only one positive eigenvalue, provided r2<qE, which implies that P1 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane. On the other hand, there are three negative eigenvalues, provided that r2>qE, which implies that P1 is a stable point.In order to discuss the local stability of model system (2) around the positive equilibrium P*(x1*,x2*,y*), (33) reduces to(36)|λ+dB*-A*r1(1+B*e-λτ(z*))-A*-A*A*λ+A*+B*e-λτ(z*)+2bx2*+a1k1y*(x2*+k1)2A*+a1x2*x2*+k10-a2y*2(x2*+k2)2λ-r2+qE+2a2y*x2*+k2|=0,where A*=d(r1x2*-dx1*)τ′(z*), B*=(r1x2*-dx1*)/r1x2* and z*=x1*+x2*+y*.It can be computed that(37)M(λ)+N(λ)e-λτ(z*)=0,
where M(λ)=λ3+m1λ2+m2λ+m3 and N(λ)=n1λ2+n2λ+n3,
(38)m1=dB*+2bx2*+a1k1y*(x2*+k1)2-r2+qE+2a2y*x2*+k2,m2=(2a2y*x2*+k2-r2+qE)(A*+2bx2*+a1k1y*(x2*+k1)2)-A*(r1-A*),m3=a2y*2A*2(x2*+k2)2+dB*(A*+2bx2*+a1k1y*(x2*+k1)2)×(2a2y*x2*+k2-r2+qE)-A*(r1+2bx2*+a1k1y*(x2*+k1)2)×(2a2y*x2*+k2-r2+qE),n1=B*,n2=B*(2a2y*x2*+k2-r2+qE-r1A*),n3=B*(2a2y*x2*+k2-r2+qE)(dB*-A*-r1A*).In the following part, dynamical behavior of model system (2) around the positive equilibrium P* is investigated. Furthermore, local stability analysis is discussed due to the variation of maturation delay and harvest effort level. By taking τ′(z*) as a bifurcation parameter, conditions for local stability switch are discussed with the increase of τ′(z*) from zero.Case 1(τ′(z*)=0).
In the case ofτ′(z*)=0, it derives that τ(z*) remains as a constant (zero or a positive constant) for all time t>0 based on (H2). In the following part, τ(z*) is denoted as τ* for simplifying. Furthermore, it can be computed that A*=0, and mi,ni (i=1,2,3) in (37) can be rewritten as follows:
(39)m^1=dB*+2bx2*+a1k1y*(x2*+k1)2-r2+qE+2a2y*x2*+k2,m^2=(2a2y*x2*+k2-r2+qE)(2bx2*+a1k1y*(x2*+k1)2),m^3=dB*(2bx2*+a1k1y*(x2*+k1)2)(2a2y*x2*+k2-r2+qE),n^1=B*,n^2=B*(2a2y*x2*+k2-r2+qE),n^3=dB*2(2a2y*x2*+k2-r2+qE).Theorem 6.
Supposing that hypotheses (H1)–(H4) hold, ifr2-qE>0, then model system (2) is stable around the positive equilibrium P* in the case of τ*=0.Proof.
Whenτ*=0, (37) can be rewritten as follows:
(40)λ3+(m^1+n^1)λ2+(m^2+n^2)λ+m^3+n^3=0.
Based on the above analysis, it can be concluded that the roots of (40) have negative real parts by using the Routh-Hurwitz criteria [1]. Consequently, P* is locally stable in the case of τ*=0.Whenτ*>0, let λ=iω be a root of (37), where ω is positive. Substitute λ=iω into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows:
(41)ω3-m^2ω=(n^1ω2-n^3)sin(ωτ*)+n^2ωcos(ωτ*),m^1ω2-m^3=n^2ωsin(ωτ*)-(n^1ω2-n^3)cos(ωτ*).By squaring and adding (41), it can be obtained that
(42)ω6+B1ω4+B2ω2+B3=0,
where B1=m^12-2m^2-n^12, B2=m^22-2m^1m^3+2n^1n^3-n^22, B3=m^32-n^32, and m^i, n^i (i=1,2,3) have been defined in (40).According to the values ofBi (i=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is B3<0.IfB3<0 holds, then (42) has a positive root ω0, and (37) has a pair of purely imaginary roots of the form ±iω0. Consequently, it can be obtained by eliminating sin(ωτ*) from (41):
(43)cos(ωτ*)666=(n^2-m^1n^1)ω4+(m^1n^3+m^3n^1-m^2n^2)ω2-m^3n^3(n^2ω2)2+(n^3-n^1ω2)2,
The τ~k corresponding to ω0 is as follows:
(44)τ~k=1ω0arccos[((n^2ω2)2+(n^3-n^1ω2)2)-1((n^2-m^1n^1)ω46666666666+(m^1n^3+m^3n^1-m^2n^2)ω2-m^3n^3)666666666×((n^2ω2)2+(n^3-n^1ω2)2)-1]+2kπω0,k=0,1,2,…. By virtue of Butler’s lemma [37], it can be concluded that the positive equilibrium P* remains locally stable for τ*<τ~0, as k=0.Case 2(τ′(z*)>0).
In the case ofτ′(z*)>0, local stability of model system (2) around the positive equilibrium P* can change only if there exists at least one root of (37) such that Reλ=0.
Letλ=iν be one such root, where ν is positive. Substitute λ=iν into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows:
(45)ν3-m2ν=(n1ν2-n3)sin(ντ(z*))+n2νcos(ντ(z*)),m1ν2-m3=n2νsin(ντ(z*))-(n1ν2-n3)cos(ντ(z*)).
By squaring and adding (45), it can be obtained that
(46)ν6+C1ν4+C2ν2+C3=0,
where C1=m12-2m2-n12,C2=m22-2m1m3+2n1n3-n22, C3=m32-n32, and mi, ni (i=1,2,3) have been defined in (37).
According to the values ofCi (i=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is C3<0, which derives that
(47)τ′(z*)66>((x2*+k1)2a1k1r1(r2-qE)x2*(x2*+k2)66666+a2(x2*+k1)2(r1x2*(2bx2*-1)+dx1*))66666×(+r1x2*(2bx2*-d-1))r1x2*[a1k1r1(r2-qE)x2*(x2*+k2)6666666666666+a2(x2*+k1)26666666666666×(dx1*(d+2)+r1x2*(2bx2*-d-1))])-1.
If the above inequality holds, then model system (2) is unstable around the positive equilibrium P* in the case of τ′(z*)>0.


### 3.4. Global Stability Analysis of Positive Equilibrium

In this section, global stability of the positive equilibriumP* is discussed by using an iterative technique in the case of τ′(z*)=0.Lemma 7 (see [29]).
Consider the following equation:(48)x˙=ax(t-τ)-bx(t)-cx2(t),
where a, b, c, and τ are positive constants, and x(t)>0 for t∈[-τ,0]; it follows that (i)
Ifa>b, then limt→+∞x(t)=(a-b)/c;(ii)
Ifa<b, then limt→+∞x(t)=0.Theorem 8.
Supposing that hypotheses (H1)–(H4) andr2-qE>0 hold, if the following inequalities hold
(49)bk2+r1e-dτ*<a2by*,a2k1r1e-dτ*>a1(r2-qE)(R1+k2),
then the positive equilibrium P* is globally asymptotically stable in the case of τ′(z*)=0.Proof.
In the case ofτ′(z*)=0, it derives that τ(z*) remains as a constant (zero or a positive constant) for all time t>0 based on (H2). In the following part, τ(z*) is denoted as τ* for simplifying. Let
(50)U1=limt→+∞supx2(t),V1=limt→+∞infx2(t),U2=limt→+∞supy(t),V2=limt→+∞infy(t).
In the following, we will claim thatU1=V1=x2*, U2=V2=y*.
It follows from Theorem2 that x2(t)≤R1 (R1 has been defined in Theorem 2). From model system (2),
(51)y˙(t)≤(r2-qE)y(t)-a2y2(t)R1+k2.
By standard comparison argument, it derives that(52)U2≤(r2-qE)(R1+k2)a2∶=J1y,
and then for sufficiently small ϵ>0, there exists a T11>0 such that if t>T11, y(t)≤J1y+ϵ. Based on Theorem 1, x2(t)+k1>k1, it can be obtained that for t>T11+τ*,
(53)x˙2(t)≥r1e-dτ*x2(t-τ*)-bx22(t)-a1(J1y+ϵ)x2(t)x2(t)+k1>r1e-dτ*x2(t-τ*)-bx22(t)-a1(J1y+ϵ)x2(t)k1.
Consider the following auxiliary equation:(54)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(J1y+ϵ)v(t)k1.
Under the conditiona2k1r1e-dτ*>a1(r2-qE)(R1+k2), it follows from Lemma 7 that
(55)limt→+∞v(t)=k1r1e-dτ*-a1(J1y+ϵ)bk1:=I1x.
Hence,V1≥I1x. For sufficiently small ϵ>0, there exits T12≥T11+τ* such that if t>T22, then x2(t)≥I1x-ϵ.
We derive from the model system (2) that for t>T12,
(56)y˙(t)≥(r2-qE)y(t)-a2y2(t)k2+I1x-ϵ.
A standard comparison argument shows that(57)limt→+∞y(t)=(r2-qE)(k2+I1x-ϵ)a2:=I1y.
Hence, for sufficiently smallϵ>0, there is a T21≥T12 satisfying if t>T21, then y(t)≥I1y-ϵ. Consequently, for t>T21+τ*,
(58)x˙2(t)≤r1e-dτ*x2(t-τ*)-bx22(t)-a1(I1y-ϵ)x2(t)R1+k1.
Consider the following auxiliary equation:(59)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(I1y-ϵ)v(t)R1+k1.
It follows from Lemma7 that
(60)limt→+∞v(t)=(r1e-dτ*-b)(R1+k1)a1(I1y-ϵ)∶=J1x.
Hence,U1≤J1x. For sufficiently small ϵ>0, there exists a T22≥T21+τ* satisfying that if t>T22, then x2(t)≤J1x+ϵ. For t>T22, it gives that
(61)y˙(t)≤(r2-qE)y(t)-a2y2(t)k2+J1x+ϵ.
By standard comparison argument, it derives that(62)limt→+∞y(t)=(r2-qE)(J1x+k2+ϵ)a2∶=J2y.
Hence, for sufficiently smallϵ>0, there exists T31≥T22 satisfying that if t>T31, y(t)≤J2y+ϵ, the for t>T31+τ*(63)x˙2(t)≥r1e-dτ*x2(t-τ*)-bx22(t)-a1(J2y+ϵ)x2(t)k1.
Consider the following auxiliary equation:(64)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(J2y+ϵ)v(t)k1.
By using Lemma7, it can be obtained that
(65)limt→+∞v(t)=k1r1e-dτ*-a1(J2y+ϵ)bk1:=I2x.
Since it is true for any sufficiently smallϵ>0, V1≥I2x. Therefore, there exists T32≥T31+τ* such that if t>T32, then x2(t)≥I2x-ϵ.
It follows from model system (2) that for t>T32,
(66)y˙(t)≥(r2-qE)y(t)-a2y2(t)I2x-ϵ+k2.
By using standard comparison argument, it derives that(67)limt→+∞y(t)=(r2-qE)(I2x-ϵ+k2)a2:=I2y.
Since this is true for any sufficiently smallϵ>0, V2≥I2y. Consequently, there exists T41≥T32 satisfying if t>T41, then y(t)≥I2y-ϵ.
It follows from model system (2) that for t>T41+τ*,
(68)x˙2(t)≤r1e-dτ*x2(t-τ*)-bx22(t)-a1(I2y-ϵ)x2(t)R1+k1.
Consider the following auxiliary equation,(69)v˙(t)≤r1e-dτ*v(t-τ*)-bv2(t)-a1(I2y-ϵ)v(t)R1+k1.
By using Lemma7, it derives that
(70)limt→+∞x2(t)=r1e-dτ*(R1+k1)-a1(I2y-ϵ)b(R1+k1):=J2x.
Continuing the above process, four sequences{Inx}, {Iny}, {Jnx}, {Jny}, n=1,2,…, are obtained which take the following form
(71)Jnx=r1e-dτ*(R1+k1)+a1ϵ-a1Inyb(R1+k1),Jny=(ϵ+k2)(r2-qE)+(r2-qE)Jn-1xa2,Inx=k1r1e-dτ*-a1ϵ-a1Jnybk1,Iny=(r2-qE)(k2-ϵ)+(r2-qE)Inxa2.
It is easy to show that(72)Inx≤V1≤U1≤Jnx,Iny≤V2≤U2≤Jny.
By virtue of (71), it derives that
(73)Jny6=(r2-qE)[b(ϵ+k2)+r1e-dτ*]a2b+a1ϵ(r2-qE)a2b(R1+k1)66-a1(r2-qE)2[k1(r1e-dτ*+bk2)+ϵ(bk1-a1)-a1Jn-1y]a22b2k1(R1+k1).
Furthermore,(74)Jny-Jn-1y666=(r2-qE)[b(ϵ+k2)+r1e-dτ*]a2b+a1ϵ(r2-qE)a2b(R1+k1)666-a1(r2-qE)2[k1(r1e-dτ*+bk2)+ϵ(bk1-a1)-a1Jn-1y]a22b2k1(R1+k1)666-Jn-1y.
If the following inequalities hold(75)bk2+r1e-dτ*<a2by*,a2k1r1e-dτ*>a1(r2-qE)(R1+k2),
then Jny-Jn-1y≤0, which implies that {Jny∣Jny≥y*,n=1,2,…} is monotonically decreasing. Hence, it can be shown that limitation of sequence {Jny} exists. Taking n→+∞, it follows from (73) that
(76)limn→+∞Jny=y*.
By using (71) and (76), it can be shown that
(77)limn→+∞Jnx=x2*,limn→+∞Iny=y*,limn→+∞Inx=x2*.
According to the definition ofU1, U2, V1, and V2, it derives that
(78)U1=V1=x2*,U2=V2=y*.
Hence,(79)limt→+∞x2(t)=x2*,limt→+∞y(t)=y*.
Based on (5), it derives that
(80)x1(t)=∫t-τ^t‍r1e-d(t-s)ϕ(s)ds.
By using L’Hospital’s rule, it derives that(81)limt→+∞x1(t)=limt→+∞r1[edtϕ(t)-ed(t-τ^)ϕ(t-τ^)]dedt=limt→+∞r1d[ϕ(t)-e-dτ^ϕ(t-τ^)]=limt→+∞r1d(1-e-dτ^)x2*.
According to (20), it is easy to show that
(82)limt→+∞x1(t)=limt→+∞r1d(1-e-dτ^)x2*=x1*.
This completes the proof.


## 3.1. Positiveness and Boundness of Solutions

Theorem 1.
Under hypotheses (H1)–(H4), solutions of model system (2) with given initial conditions are positive for all t>0.Proof.
Assume that there existst0=inf{t>0∣x2(t)=0}. Based on the continuity t0>0, it can be computed by evaluating the model system (2) at time t0:
(6)x˙2(t0)={r1e-dτ(z(t0))ϕ(t0-τ(z(t0))),0≤t0≤τ^,r1e-dτ(z(t0))x2(t0-τ(z(t0))),t0>τ^.
According to (6) and the initial conditions of model system (2), it is easy to show that x˙2(t0)>0. On the other hand, it follows from the definition of t0 that x˙2(t0)=0, which is a contradiction. Consequently, x2(t)>0 for all t>0.
Based on the positiveness ofx2(t) and the third equation of model system (2), it is easy to show that y(t)>0 for y(0)>0,t>0.
Consider the equation(7)u˙(t)=-du(t)-r1e-dτ(z)x2(t-τ(z)),u(0)=x1(0)>0.
It is obvious to show thatu˙(t)<0; that is, u(t) is strictly decreasing. By virtue of positiveness of x2(t), y(t), it derives that
(8)x1(t)>u(t),0<t≤τ^.
By solving (7), it gives that
(9)u(t)=e-dtu(0)-r1e-dt∫0t‍edqe-dτ(z(q))x2(q-τ(z(q)))dq.
According to (5), it derives that
(10)u(t)=e-dt∫-τ^0‍r1edsϕ(s)ds-r1e-dt∫0t‍edqe-dτ(z(q))x2(q-τ(z(q)))dq.
By substitutingp=q-τ(z(q)) in the above equation, it can be obtained that
(11)u(t)=r1e-dt∫-τ^0‍edsϕ(s)ds-r1e-dt∫-τ^t-τ(z(t))‍edpx2(p)1-τ′(z)z˙(p)dp,
which implies that
(12)u(τ^)=r1e-dτ^∫-τ^0‍edsϕ(s)ds-r1e-dτ^∫-τ^τ^-τ(z(τ^))‍edpx2(p)1-τ′(z)z˙(p)dp.
According tox2(t)>0,y(t)>0 for all t>0 and t-τ(z(t)) is an increasing function based on (H4). 1-τ′(z)z˙(t)>0 holds for -τ^≤t≤τ^-τ(z(τ^)), and the following inequality can be obtained:
(13)u(τ^)≥r1e-dτ^∫-τ^0‍edpϕ(p)(1-11-τ′(z)z˙(p))dp>0.
since u(τ^)>0 and u(t) is strictly decreasing, x1(t)>0, -τ^≤t≤0. By repeating this argument to include all positive time, it can be shown that x1(t)>0 for all t>0. Hence, solutions of model system (2) with given initial conditions are positive for all t>0.Theorem 2.
If the hypotheses (H1)–(H4) hold andr2>qE, all solutions of model system (2) are bounded within a region Ω:
(14)Ω={(x1(t),x2(t),y(t))∣60<x1(t)+x2(t)≤R1,0<y(t)≤R2(x1(t),x2(t),y(t))},
where R1=(r1+d)2/4bd,R2=(r2-qE)(R1+k2)/a2.Proof.
Letv(t)=x1(t)+x2(t), and it is easy to show that v(t)>0 based on positiveness of solutions of model system (2). By calculating the derivative of v(t) along the solutions, it gives that
(15)v˙(t)=r1x2(t)-dx1(t)-bx22(t)-a1x2(t)y(t)x2(t)+k1<(r1+d)x2(t)-bx22-dv(t).
By using the standard comparison principle in (15), it derives that
(16)limt→∞sup(v(t))≤(r1+d)24bd=R1.
It follows from the third equation of model system (2) that
(17)y˙(t)=(r2-qE)y(t)-a2y2(t)x2(t)+k2≤(r2-qE)y(t)-a2y2(t)R1+k2.
By utilizing the standard comparison principle in inequality (17) and r2-qE>0, it gives that
(18)limt→∞supy(t)≤(r2-qE)(R1+k2)a2=R2.
Consequently, all solutions of model system (2) are bounded within a region Ω:
(19)Ω={(x1(t),x2(t),y(t))∣60<x1(t)+x2(t)≤R1,0<y(t)≤R2(x1(t),x2(t),y(t))},
where R1=(r1+d)2/4bd, R2=(r2-qE)(R1+k2)/a2.


## 3.2. Existence of Equilibria and Uniqueness of Positive Equilibrium

The existence of biologically reasonable equilibria of model system (2) is investigated in this subsection. Since the biological interpretation of the positive equilibrium implies that immature prey, mature prey, and predator population all exist, uniqueness of positive equilibrium is also studied.By simple computation, there are two equilibria denoted byP0(0,0,0) and P1(0,0,k2(r2-qE)/a2). The biological interpretations of P0, P1 are as follows. For P0(0,0,0), it implies that all population in harvested prey predator ecosystem does not exist. For P1(0,0,k2(r2-qE)/a2), it implies that there is not any predation source for predator population. It follows from the previous biological interpretations that population in such ecosystem cannot be maintained at an ideal level for sustainable development, which are not relevant to major investigation in this paper.Furthermore, there is one or more positive equilibria denoted byP*(x1*,x2*,y*). In order to discuss the existence of P*, it is equivalent to show that the following equations always have at least one positive solution:
(20)r1x2-dx1-r1e-dτ(x1+x2+y)x2=0,r1e-dτ(x1+x2+y)-bx2-a1yx2+k1=0,r2-qE-a2yx2+k2=0.It follows from (20) that
(21)y=(r2-qE)(x2+k2)a2=f(x2),(22)r1x2-dx1=r1x2e-dτ(x1+g(x2)),(23)r1e-dτ(x1+g(x2))=bx2+a1h(x2),
where g(x2)=x2+y=x2+f(x2) and h(x2)=y/(x2+k1)=f(x2)/(x2+k1).Theorem 3 (existence of positive equilibrium).
Supposing that hypotheses (H1)–(H4) hold, ifk1≥k2, r2>qE, and a2r1e-dτ1>a1r2, then there exists at least one positive equilibrium P*.Proof.
LetΓ1 and Γ2 be the solution curves of (22) and (23) for x1≥0, x2≥0, respectively. The analytical properties of curve Γ1 and Γ2 are as follows.
ForΓ1: by simple computing, it can be found that (0,0)∈Γ1.
According to (H2) and positiveness of all solutions of model system (2), it is easy to show that limx1→+∞τ(x1(t)+g(x2))=τ1, and
(24)limx1→+∞x2(x1)=limx1→+∞dx1r1(1-e-dτ(x1+g(x2))=limx1→+∞dx1r1(1-e-dτ1)=+∞.
ForΓ2: by differentiating x2 against x1 along Γ2, it can be obtained that
(25)dx2dx1=-dr1e-dτ(x1+g(x2))×(a1(r2-qE)(k1-k2)a2(x2+k1)2dr1e-dτ(x1+g(x2))τ′(x1+g(x2))a2+r2-qEa26666+b+a1(r2-qE)(k1-k2)a2(x2+k1)2)-1.
It can be shown that(dx2/dx1)<0, provided that k1≥k2, r2>qE, and then Γ2 is strictly decreasing.
Furthermore, according tok1≥k2, r2>qE, and a2r1e-dτ1>a1(r2-qE),
(26)limx1→∞x2(x1)=1b(r1e-dτ1-a1(r2-qE)(x2+k2)a2(x2+k1))>0.
Based on the above analysis,Γ1 and Γ2 intersect at some positive values, which proves the existence of positive equilibrium P*.Theorem 4 (uniqueness of positive equilibrium).
Supposing that hypotheses (H1)–(H4) hold, if the following inequality holds(27)1-dτ′(x1*+g(x2*))(bx2*+a1h(x2*))66×(a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*))66666×(bx2*+a1h(x2*))a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*)))+a1dh′(x2*)66+dτ′(x1*+g(x2*)g′(x2*)(bx2*+a1h(x2*))66×(d+b+dτ′(x1*+g(x2*)))>0,
then there exists a unique positive equilibrium.Proof.
Based on (22) and (23), x2 can be defined as the function of x1:
(28)Γ1:x2=g1(x1),Γ2:x2=g2(x1).
The positive equilibrium P* will be unique, provided that g1′(x1)>g2′(x1) for every such P* otherwise reverse inequality holds.
By differentiating (22) with respect to x1, it can be obtained that
(29)g1′(x1)=d+r1e-dτ(x1+g(x2))(1-dτ′(x1+g(x2)))r1(1+e-dτ(x1+g(x2))dτ′(x1+g(x2))g′(x2)).
By differentiating (23) with respect to x1, it can be obtained that
(30)g2′(x1)666=-b+dr1τ′(x1+g(x2))e-dτ(x1+g(x2))a1h′(x2)+dr1g′(x2)τ′(x1+g(x2))e-dτ(x1+g(x2)).
On the other hand, some expressions about positive equilibriumP*(x1*,x2*,y*) can be obtained based on (22) and (23),
(31)r1x2*-dx1*=r1x2*e-dτ(x1*+g(x2*)),r1e-dτ(x1*+g(x2*))=bx2*+a1h(x2*).
According to (31), g1′(x1*)>g1′(x2*) is equivalent to the following inequality:
(32)1-dτ′(x1*+g(x2*))(bx2*+a1h(x2*))66×(a1h′(x2*)+dg′(x2*)τ′(x1*+g(x2*))666666×(bx2*+a1h(x2*))a1h′(x2*))+a1dh′(x2*)66+dτ′(x1*+g(x2*)g′(x2*)66×(bx2*+a1h(x2*))(d+b+dτ′(x1*+g(x2*)))>0.
This completes the proof.


## 3.3. Local Stability Analysis around Equilibria

Local stability of model system (2) around all equilibria of model system (2) is investigated. Furthermore, stability switch due to variation of maturity delay and harvest effort level is also studied in this subsection.The characteristic equation of model system (2) about some equilibrium P~=(x~1,x~2,y~) takes the following form: (33)|λ+(d-A~)e-dτ(z~)r1(1+e-(λ+d)τ(z~))-A~e-dτ(z~)-A~e-dτ(z~)A~e-dτ(z~)λ+2bx~+a1k1y~(x~2+k1)2+A~e-dτ(z~)+r1e-(λ+d)τ(z~)A~e-dτ(z~)+a1x~2k1+x~20-a2y~2(k2+x~2)2λ-r2+qE+2a2y~k2+x~2|=0,where A~=dr1x~2τ′(z~).Theorem 5.
Local stability analysis of model system (2) around P0 and P1 is as follows: (a)
ifr2<qE, then model system is locally stable around P0, and P1 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane;(b)
ifr2>qE, then model system is locally stable around P1, and P0 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane.Proof.
ForP0(0,0,0), (33) reduces to
(34)(λ+de-dτ(0))(λ+r1e-(λ+d)τ(0))(λ-(r2-qE))=0.
By solving (34), it can be found that there are two negative eigenvalues and only one positive eigenvalue, provided r2>qE, which implies that P0 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane. On the other hand, there are three negative eigenvalues, provided that r2<qE, which implies that P0 is a stable point.
ForP1(0,0,(k2(r2-qE)/a2)), (33) reduces to
(35)(λ+de-dτ(y~))(λ+r1e-(λ+d)τ(y~)+a1y~k1)(λ+(r2-qE))=0,
where y~=k2(r2-qE)/a2. It follows from (35) that there are two negative eigenvalues and only one positive eigenvalue, provided r2<qE, which implies that P1 is a saddle point which is unstable in the y-direction and stable in the x1-x2 plane. On the other hand, there are three negative eigenvalues, provided that r2>qE, which implies that P1 is a stable point.In order to discuss the local stability of model system (2) around the positive equilibrium P*(x1*,x2*,y*), (33) reduces to(36)|λ+dB*-A*r1(1+B*e-λτ(z*))-A*-A*A*λ+A*+B*e-λτ(z*)+2bx2*+a1k1y*(x2*+k1)2A*+a1x2*x2*+k10-a2y*2(x2*+k2)2λ-r2+qE+2a2y*x2*+k2|=0,where A*=d(r1x2*-dx1*)τ′(z*), B*=(r1x2*-dx1*)/r1x2* and z*=x1*+x2*+y*.It can be computed that(37)M(λ)+N(λ)e-λτ(z*)=0,
where M(λ)=λ3+m1λ2+m2λ+m3 and N(λ)=n1λ2+n2λ+n3,
(38)m1=dB*+2bx2*+a1k1y*(x2*+k1)2-r2+qE+2a2y*x2*+k2,m2=(2a2y*x2*+k2-r2+qE)(A*+2bx2*+a1k1y*(x2*+k1)2)-A*(r1-A*),m3=a2y*2A*2(x2*+k2)2+dB*(A*+2bx2*+a1k1y*(x2*+k1)2)×(2a2y*x2*+k2-r2+qE)-A*(r1+2bx2*+a1k1y*(x2*+k1)2)×(2a2y*x2*+k2-r2+qE),n1=B*,n2=B*(2a2y*x2*+k2-r2+qE-r1A*),n3=B*(2a2y*x2*+k2-r2+qE)(dB*-A*-r1A*).In the following part, dynamical behavior of model system (2) around the positive equilibrium P* is investigated. Furthermore, local stability analysis is discussed due to the variation of maturation delay and harvest effort level. By taking τ′(z*) as a bifurcation parameter, conditions for local stability switch are discussed with the increase of τ′(z*) from zero.Case 1(τ′(z*)=0).
In the case ofτ′(z*)=0, it derives that τ(z*) remains as a constant (zero or a positive constant) for all time t>0 based on (H2). In the following part, τ(z*) is denoted as τ* for simplifying. Furthermore, it can be computed that A*=0, and mi,ni (i=1,2,3) in (37) can be rewritten as follows:
(39)m^1=dB*+2bx2*+a1k1y*(x2*+k1)2-r2+qE+2a2y*x2*+k2,m^2=(2a2y*x2*+k2-r2+qE)(2bx2*+a1k1y*(x2*+k1)2),m^3=dB*(2bx2*+a1k1y*(x2*+k1)2)(2a2y*x2*+k2-r2+qE),n^1=B*,n^2=B*(2a2y*x2*+k2-r2+qE),n^3=dB*2(2a2y*x2*+k2-r2+qE).Theorem 6.
Supposing that hypotheses (H1)–(H4) hold, ifr2-qE>0, then model system (2) is stable around the positive equilibrium P* in the case of τ*=0.Proof.
Whenτ*=0, (37) can be rewritten as follows:
(40)λ3+(m^1+n^1)λ2+(m^2+n^2)λ+m^3+n^3=0.
Based on the above analysis, it can be concluded that the roots of (40) have negative real parts by using the Routh-Hurwitz criteria [1]. Consequently, P* is locally stable in the case of τ*=0.Whenτ*>0, let λ=iω be a root of (37), where ω is positive. Substitute λ=iω into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows:
(41)ω3-m^2ω=(n^1ω2-n^3)sin(ωτ*)+n^2ωcos(ωτ*),m^1ω2-m^3=n^2ωsin(ωτ*)-(n^1ω2-n^3)cos(ωτ*).By squaring and adding (41), it can be obtained that
(42)ω6+B1ω4+B2ω2+B3=0,
where B1=m^12-2m^2-n^12, B2=m^22-2m^1m^3+2n^1n^3-n^22, B3=m^32-n^32, and m^i, n^i (i=1,2,3) have been defined in (40).According to the values ofBi (i=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is B3<0.IfB3<0 holds, then (42) has a positive root ω0, and (37) has a pair of purely imaginary roots of the form ±iω0. Consequently, it can be obtained by eliminating sin(ωτ*) from (41):
(43)cos(ωτ*)666=(n^2-m^1n^1)ω4+(m^1n^3+m^3n^1-m^2n^2)ω2-m^3n^3(n^2ω2)2+(n^3-n^1ω2)2,
The τ~k corresponding to ω0 is as follows:
(44)τ~k=1ω0arccos[((n^2ω2)2+(n^3-n^1ω2)2)-1((n^2-m^1n^1)ω46666666666+(m^1n^3+m^3n^1-m^2n^2)ω2-m^3n^3)666666666×((n^2ω2)2+(n^3-n^1ω2)2)-1]+2kπω0,k=0,1,2,…. By virtue of Butler’s lemma [37], it can be concluded that the positive equilibrium P* remains locally stable for τ*<τ~0, as k=0.Case 2(τ′(z*)>0).
In the case ofτ′(z*)>0, local stability of model system (2) around the positive equilibrium P* can change only if there exists at least one root of (37) such that Reλ=0.
Letλ=iν be one such root, where ν is positive. Substitute λ=iν into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows:
(45)ν3-m2ν=(n1ν2-n3)sin(ντ(z*))+n2νcos(ντ(z*)),m1ν2-m3=n2νsin(ντ(z*))-(n1ν2-n3)cos(ντ(z*)).
By squaring and adding (45), it can be obtained that
(46)ν6+C1ν4+C2ν2+C3=0,
where C1=m12-2m2-n12,C2=m22-2m1m3+2n1n3-n22, C3=m32-n32, and mi, ni (i=1,2,3) have been defined in (37).
According to the values ofCi (i=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is C3<0, which derives that
(47)τ′(z*)66>((x2*+k1)2a1k1r1(r2-qE)x2*(x2*+k2)66666+a2(x2*+k1)2(r1x2*(2bx2*-1)+dx1*))66666×(+r1x2*(2bx2*-d-1))r1x2*[a1k1r1(r2-qE)x2*(x2*+k2)6666666666666+a2(x2*+k1)26666666666666×(dx1*(d+2)+r1x2*(2bx2*-d-1))])-1.
If the above inequality holds, then model system (2) is unstable around the positive equilibrium P* in the case of τ′(z*)>0.


## 3.4. Global Stability Analysis of Positive Equilibrium

In this section, global stability of the positive equilibriumP* is discussed by using an iterative technique in the case of τ′(z*)=0.Lemma 7 (see [29]).
Consider the following equation:(48)x˙=ax(t-τ)-bx(t)-cx2(t),
where a, b, c, and τ are positive constants, and x(t)>0 for t∈[-τ,0]; it follows that (i)
Ifa>b, then limt→+∞x(t)=(a-b)/c;(ii)
Ifa<b, then limt→+∞x(t)=0.Theorem 8.
Supposing that hypotheses (H1)–(H4) andr2-qE>0 hold, if the following inequalities hold
(49)bk2+r1e-dτ*<a2by*,a2k1r1e-dτ*>a1(r2-qE)(R1+k2),
then the positive equilibrium P* is globally asymptotically stable in the case of τ′(z*)=0.Proof.
In the case ofτ′(z*)=0, it derives that τ(z*) remains as a constant (zero or a positive constant) for all time t>0 based on (H2). In the following part, τ(z*) is denoted as τ* for simplifying. Let
(50)U1=limt→+∞supx2(t),V1=limt→+∞infx2(t),U2=limt→+∞supy(t),V2=limt→+∞infy(t).
In the following, we will claim thatU1=V1=x2*, U2=V2=y*.
It follows from Theorem2 that x2(t)≤R1 (R1 has been defined in Theorem 2). From model system (2),
(51)y˙(t)≤(r2-qE)y(t)-a2y2(t)R1+k2.
By standard comparison argument, it derives that(52)U2≤(r2-qE)(R1+k2)a2∶=J1y,
and then for sufficiently small ϵ>0, there exists a T11>0 such that if t>T11, y(t)≤J1y+ϵ. Based on Theorem 1, x2(t)+k1>k1, it can be obtained that for t>T11+τ*,
(53)x˙2(t)≥r1e-dτ*x2(t-τ*)-bx22(t)-a1(J1y+ϵ)x2(t)x2(t)+k1>r1e-dτ*x2(t-τ*)-bx22(t)-a1(J1y+ϵ)x2(t)k1.
Consider the following auxiliary equation:(54)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(J1y+ϵ)v(t)k1.
Under the conditiona2k1r1e-dτ*>a1(r2-qE)(R1+k2), it follows from Lemma 7 that
(55)limt→+∞v(t)=k1r1e-dτ*-a1(J1y+ϵ)bk1:=I1x.
Hence,V1≥I1x. For sufficiently small ϵ>0, there exits T12≥T11+τ* such that if t>T22, then x2(t)≥I1x-ϵ.
We derive from the model system (2) that for t>T12,
(56)y˙(t)≥(r2-qE)y(t)-a2y2(t)k2+I1x-ϵ.
A standard comparison argument shows that(57)limt→+∞y(t)=(r2-qE)(k2+I1x-ϵ)a2:=I1y.
Hence, for sufficiently smallϵ>0, there is a T21≥T12 satisfying if t>T21, then y(t)≥I1y-ϵ. Consequently, for t>T21+τ*,
(58)x˙2(t)≤r1e-dτ*x2(t-τ*)-bx22(t)-a1(I1y-ϵ)x2(t)R1+k1.
Consider the following auxiliary equation:(59)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(I1y-ϵ)v(t)R1+k1.
It follows from Lemma7 that
(60)limt→+∞v(t)=(r1e-dτ*-b)(R1+k1)a1(I1y-ϵ)∶=J1x.
Hence,U1≤J1x. For sufficiently small ϵ>0, there exists a T22≥T21+τ* satisfying that if t>T22, then x2(t)≤J1x+ϵ. For t>T22, it gives that
(61)y˙(t)≤(r2-qE)y(t)-a2y2(t)k2+J1x+ϵ.
By standard comparison argument, it derives that(62)limt→+∞y(t)=(r2-qE)(J1x+k2+ϵ)a2∶=J2y.
Hence, for sufficiently smallϵ>0, there exists T31≥T22 satisfying that if t>T31, y(t)≤J2y+ϵ, the for t>T31+τ*(63)x˙2(t)≥r1e-dτ*x2(t-τ*)-bx22(t)-a1(J2y+ϵ)x2(t)k1.
Consider the following auxiliary equation:(64)v˙(t)=r1e-dτ*v(t-τ*)-bv2(t)-a1(J2y+ϵ)v(t)k1.
By using Lemma7, it can be obtained that
(65)limt→+∞v(t)=k1r1e-dτ*-a1(J2y+ϵ)bk1:=I2x.
Since it is true for any sufficiently smallϵ>0, V1≥I2x. Therefore, there exists T32≥T31+τ* such that if t>T32, then x2(t)≥I2x-ϵ.
It follows from model system (2) that for t>T32,
(66)y˙(t)≥(r2-qE)y(t)-a2y2(t)I2x-ϵ+k2.
By using standard comparison argument, it derives that(67)limt→+∞y(t)=(r2-qE)(I2x-ϵ+k2)a2:=I2y.
Since this is true for any sufficiently smallϵ>0, V2≥I2y. Consequently, there exists T41≥T32 satisfying if t>T41, then y(t)≥I2y-ϵ.
It follows from model system (2) that for t>T41+τ*,
(68)x˙2(t)≤r1e-dτ*x2(t-τ*)-bx22(t)-a1(I2y-ϵ)x2(t)R1+k1.
Consider the following auxiliary equation,(69)v˙(t)≤r1e-dτ*v(t-τ*)-bv2(t)-a1(I2y-ϵ)v(t)R1+k1.
By using Lemma7, it derives that
(70)limt→+∞x2(t)=r1e-dτ*(R1+k1)-a1(I2y-ϵ)b(R1+k1):=J2x.
Continuing the above process, four sequences{Inx}, {Iny}, {Jnx}, {Jny}, n=1,2,…, are obtained which take the following form
(71)Jnx=r1e-dτ*(R1+k1)+a1ϵ-a1Inyb(R1+k1),Jny=(ϵ+k2)(r2-qE)+(r2-qE)Jn-1xa2,Inx=k1r1e-dτ*-a1ϵ-a1Jnybk1,Iny=(r2-qE)(k2-ϵ)+(r2-qE)Inxa2.
It is easy to show that(72)Inx≤V1≤U1≤Jnx,Iny≤V2≤U2≤Jny.
By virtue of (71), it derives that
(73)Jny6=(r2-qE)[b(ϵ+k2)+r1e-dτ*]a2b+a1ϵ(r2-qE)a2b(R1+k1)66-a1(r2-qE)2[k1(r1e-dτ*+bk2)+ϵ(bk1-a1)-a1Jn-1y]a22b2k1(R1+k1).
Furthermore,(74)Jny-Jn-1y666=(r2-qE)[b(ϵ+k2)+r1e-dτ*]a2b+a1ϵ(r2-qE)a2b(R1+k1)666-a1(r2-qE)2[k1(r1e-dτ*+bk2)+ϵ(bk1-a1)-a1Jn-1y]a22b2k1(R1+k1)666-Jn-1y.
If the following inequalities hold(75)bk2+r1e-dτ*<a2by*,a2k1r1e-dτ*>a1(r2-qE)(R1+k2),
then Jny-Jn-1y≤0, which implies that {Jny∣Jny≥y*,n=1,2,…} is monotonically decreasing. Hence, it can be shown that limitation of sequence {Jny} exists. Taking n→+∞, it follows from (73) that
(76)limn→+∞Jny=y*.
By using (71) and (76), it can be shown that
(77)limn→+∞Jnx=x2*,limn→+∞Iny=y*,limn→+∞Inx=x2*.
According to the definition ofU1, U2, V1, and V2, it derives that
(78)U1=V1=x2*,U2=V2=y*.
Hence,(79)limt→+∞x2(t)=x2*,limt→+∞y(t)=y*.
Based on (5), it derives that
(80)x1(t)=∫t-τ^t‍r1e-d(t-s)ϕ(s)ds.
By using L’Hospital’s rule, it derives that(81)limt→+∞x1(t)=limt→+∞r1[edtϕ(t)-ed(t-τ^)ϕ(t-τ^)]dedt=limt→+∞r1d[ϕ(t)-e-dτ^ϕ(t-τ^)]=limt→+∞r1d(1-e-dτ^)x2*.
According to (20), it is easy to show that
(82)limt→+∞x1(t)=limt→+∞r1d(1-e-dτ^)x2*=x1*.
This completes the proof.


## 4. Numerical Simulation

With the help of MATLAB, numerical simulations are provided to understand the theoretical results which have been established in the previous sections of this paper. In order to facilitate the numerical simulation, it is assumed thatτ(z(t)) takes the following form [23]:
(83)τ(z(t))=τ0+τm-τme-z(t),
where τm∈(τ0,τ1) satisfying τ0+τm=τ1. Based on Theorem 1, it follows from simple computation that
(84)limt→0τ(z(t))=τ0,limt→+∞τ(z(t))=τ1,
which implies that (H2) holds.Values of parameters are taken from [24] which are used in Example 1 of [24] and set in appropriate units. r1=12, d=0.2, b=1.2, a1=0.5, a2=2, k1=2, k2=1, r2=2, q=0.25, and E=4. According to the given values of parameters, it follows from Theorems 3 and 4 that there exists a unique positive equilibrium P*(141.5454,5.6252,3.3126). Furthermore, it can be verified that P* is globally attractive based on Theorem 8. Responses of model system (2) are indicated in Figure 1, and the phase portrait of model system (2) with different initial values is plotted in Figure 2.Figure 1
Dynamical responses of model system (2).Figure 2
Phase portrait of model system (2) with different initial values.


## 5. Conclusion

In this paper, a harvested prey predator model is proposed to investigate the effects of density-dependent maturation delay and harvest effort on the dynamics. Conditions which influence positiveness and boundedness of solutions of model system are obtained in Theorems1 and 2, respectively. Existence of all equilibria of model system and uniqueness of the positive equilibrium are studied in Theorems 3 and 4, respectively. Biological interpretations of the positive equilibrium mean immature prey, mature prey, predator and harvest effort on predator population all exist in the harvested ecosystem. Consequently, we mainly concentrate on dynamical analysis around positive equilibrium in this paper. Local stability analysis in Theorem 6 reveals that local stability of the positive equilibrium loses due to variation of maturation delay and harvest effort level. Furthermore, global stability of the positive equilibrium is discussed by utilizing an iterative technique in Theorem 8, which is utilized to investigate the coexistence and interaction mechanism of harvested prey-predator ecosystem.Compared with the work done in [24] and the related work in [26], maturation delay for prey population in this paper relates to the density of all population within the harvested ecosystem, which accurately reflects the practical phenomena in the real world [23, 34–36]. Furthermore, it should be noted that dynamics of prey predator model with density-dependent delay for predator population is investigated in [23], while dynamics of harvest effort on population within ecosystem is not considered. Compared with the work done in [23], harvest effort on predator population is introduced, and the effect of harvesting on model dynamics is also investigated in this paper. With the rapid development of commercial harvesting on prey predator ecosystem in the real world, the introduction of harvest effort and related qualitative analysis makes the work done in this paper have some new and positive feature.



---

*Source: 329592-2013-05-12.xml*