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

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.

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) where 𝑥₁(𝑡), 𝑥₂(𝑡), and 𝑦(𝑡) represent the density of immature prey population, mature prey population and predator population, at time 𝑡, respectively; 𝑟₁ is the intrinsic growth rate of mature prey population, and 𝑑 is the death rate of immature prey population. Constant 𝜏≥0 denotes maturation delay of immature prey population to mature prey population, and the term 𝑟₁𝑒^−𝑑𝜏𝑥₂(𝑡−𝜏) represents the immature prey population who were born at time 𝑡−𝜏 and survived at time 𝑡. 𝑏 denotes the intracompetition rate for mature prey population due to overcrowding phenomenon with mature prey population. 𝑎₁ is the maximum value of the per capita reduction rate of mature prey population due to predator population, and 𝑎₂ is the maximum value of the per capita reduction rate of predator population due to mature prey population. 𝑘₁ measures the extent to which the environment provides protection to mature prey population, and 𝑘₂ measures the extent to which the environment provides protection to predator population. 𝑟₂ 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) The initial conditions for model system (2) take the following form: 𝑥₁(0)(0)𝑥>0,𝑦>0,₂(𝜃)=𝜙(𝜃)>0,𝜃∈[−𝜏,0), (3) where 𝑧(𝑡)=𝑥₁(𝑡)+𝑥₂(𝑡)+𝑦(𝑡), a scalar 𝐸≥0 denotes the harvesting effort to predator population, constant 𝑞 is the catchability coefficient of predator, and the harvesting term 𝑞𝐸𝑦(𝑡) follows the catch per unit effort hypothesis [1]. Furthermore, 𝑟₁, 𝑟₂, 𝑑, 𝑏, 𝑎₁,𝑎₂, 𝑘₁, and 𝑘₂ 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 term 𝑟₁𝑒^{−𝑑𝜏(𝑧(𝑡))}𝑥₂(𝑡−𝜏(𝑧(𝑡))) represents the immature prey population born at time 𝑡−𝜏(𝑧(𝑡)) and survive at time 𝑡 with death rate 𝑑, which represents transformation term from immature prey to mature prey.(H2) Density-dependent time delay 𝜏(𝑧(𝑡)) is taken to be an increasing differentiable bounded function of the total population (immature prey, mature prey, and predator population), which satisfies 𝑑[𝜏(𝑧(𝑡))]𝑑𝑡≥0,0≤𝜏₀≤𝜏(𝑧(𝑡))≤𝜏₁,lim_{𝑧(𝑡)→0+}𝜏(𝑧(𝑡))=𝜏₀,lim_{𝑧(𝑡)→+∞}𝜏(𝑧(𝑡))=𝜏₁. (4)(H3) For the continuity of initial conditions, it is required that 𝑥₁(0)=⁰₋𝜏𝑟₁𝑒^𝑑𝑠𝜙(𝑠)𝑑𝑠, (5) where 𝜙(𝑠) 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, 𝑡−𝜏(𝑧(𝑡)) is assumed to be a strictly increasing function of 𝑡. 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 𝑡>0.

Proof. Assume that there exists 𝑡₀=inf{𝑡>0∣𝑥₂(𝑡)=0}. Based on the continuity 𝑡₀>0, it can be computed by evaluating the model system (2) at time 𝑡₀: 𝑥̇₂𝑡₀=𝑟₁𝑒^{−𝑑𝜏(𝑧(𝑡0))}𝜙𝑡₀−𝜏𝑧𝑡₀,0≤𝑡₀≤𝑟𝜏,₁𝑒^{−𝑑𝜏(𝑧(𝑡0))}𝑥₂𝑡₀−𝜏𝑧𝑡₀,𝑡₀>𝜏. (6)
According to (6) and the initial conditions of model system (2), it is easy to show that 𝑥̇₂(𝑡₀)>0. On the other hand, it follows from the definition of 𝑡₀ that 𝑥̇₂(𝑡₀)=0, which is a contradiction. Consequently, 𝑥₂(𝑡)>0 for all 𝑡>0.
Based on the positiveness of 𝑥₂(𝑡) and the third equation of model system (2), it is easy to show that 𝑦(𝑡)>0 for 𝑦(0)>0,𝑡>0.
Consider the equation 𝑢̇(𝑡)=−𝑑𝑢(𝑡)−𝑟₁𝑒^{−𝑑𝜏(𝑧)}𝑥₂((𝑧,𝑢(0)𝑡−𝜏))=𝑥₁(0)>0. (7)
It is obvious to show that 𝑢̇(𝑡)<0; that is, 𝑢(𝑡) is strictly decreasing. By virtue of positiveness of 𝑥₂(𝑡), 𝑦(𝑡), it derives that 𝑥₁(𝑡)>𝑢(𝑡),0<𝑡≤𝜏. (8)
By solving (7), it gives that 𝑢(𝑡)=𝑒^−𝑑𝑡𝑢(0)−𝑟₁𝑒^{−𝑑𝑡𝑡}₀𝑒^𝑑𝑞𝑒^{−𝑑𝜏(𝑧(𝑞))}𝑥₂𝑞−𝜏𝑧𝑞𝑑𝑞. (9)
According to (5), it derives that 𝑢(𝑡)=𝑒^−𝑑𝑡0₋𝜏𝑟₁𝑒^𝑑𝑠𝜙(𝑠)𝑑𝑠−𝑟₁𝑒^{−𝑑𝑡𝑡}₀𝑒^𝑑𝑞𝑒^{−𝑑𝜏(𝑧(𝑞))}𝑥₂𝑞−𝜏𝑧𝑞𝑑𝑞. (10)
By substituting 𝑝=𝑞−𝜏(𝑧(𝑞)) in the above equation, it can be obtained that 𝑢(𝑡)=𝑟₁𝑒^−𝑑𝑡0₋𝜏𝑒^𝑑𝑠𝜙(𝑠)𝑑𝑠−𝑟₁𝑒^−𝑑𝑡₋^{𝑡−𝜏(𝑧(𝑡))}𝜏𝑒^𝑑𝑝𝑥₂𝑝1−𝜏^′(𝑧)𝑧̇𝑝𝑑𝑝, (11) which implies that 𝑢(𝜏)=𝑟₁𝑒^−𝑑𝜏⁰₋𝜏𝑒^𝑑𝑠𝜙(𝑠)𝑑𝑠−𝑟₁𝑒^−𝑑𝜏𝜏_{−𝜏(𝑧(}𝜏₋₎₎𝜏𝑒^𝑑𝑝𝑥₂𝑝1−𝜏^′(𝑧)𝑧̇𝑝𝑑𝑝. (12)
According to 𝑥₂(𝑡)>0,𝑦(𝑡)>0 for all 𝑡>0 and 𝑡−𝜏(𝑧(𝑡)) is an increasing function based on (H4). 1−𝜏^′(𝑧)𝑧(̇𝑡)>0 holds for −𝜏≤𝑡≤𝜏−𝜏(𝑧(𝜏)), and the following inequality can be obtained: 𝑢(𝜏)≥𝑟₁𝑒^−𝑑𝜏⁰₋𝜏𝑒^𝑑𝑝𝜙𝑝11−1−𝜏^′(𝑧)𝑧̇𝑝𝑑𝑝>0. (13) since 𝑢(𝜏)>0 and 𝑢(𝑡) is strictly decreasing, 𝑥₁(𝑡)>0, −𝜏≤𝑡≤0. By repeating this argument to include all positive time, it can be shown that 𝑥₁(𝑡)>0 for all 𝑡>0. Hence, solutions of model system (2) with given initial conditions are positive for all 𝑡>0.

Theorem 2. If the hypotheses (H1)–(H4) hold and 𝑟₂>𝑞𝐸, all solutions of model system (2) are bounded within a region Ω: Ω=𝑥₁(𝑡),𝑥₂(𝑡),𝑦(𝑡)∣0<𝑥₁(𝑡)+𝑥₂(𝑡)≤𝑅₁,0<𝑦(𝑡)≤𝑅₂, (14) where 𝑅₁=(𝑟₁+𝑑)²/4𝑏𝑑,𝑅₂=(𝑟₂−𝑞𝐸)(𝑅₁+𝑘₂)/𝑎₂.

Proof. Let 𝑣(𝑡)=𝑥₁(𝑡)+𝑥₂(𝑡), and it is easy to show that 𝑣(𝑡)>0 based on positiveness of solutions of model system (2). By calculating the derivative of 𝑣(𝑡) along the solutions, it gives that 𝑣̇(𝑡)=𝑟₁𝑥₂(𝑡)−𝑑𝑥₁(𝑡)−𝑏𝑥²₂(𝑡)−𝑎₁𝑥₂(𝑡)𝑦(𝑡)𝑥₂(𝑡)+𝑘₁<𝑟₁+𝑑𝑥₂(𝑡)−𝑏𝑥²₂−𝑑𝑣(𝑡). (15)
By using the standard comparison principle in (15), it derives that lim_𝑡→∞(𝑣(𝑡≤𝑟sup))¹+𝑑²4𝑏𝑑=𝑅₁. (16)
It follows from the third equation of model system (2) that 𝑦̇(𝑡)=𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝑥₂(𝑡)+𝑘₂≤𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝑅₁+𝑘₂. (17)
By utilizing the standard comparison principle in inequality (17) and 𝑟₂−𝑞𝐸>0, it gives that lim_𝑡→∞(𝑡)≤𝑟sup𝑦²−𝑞𝐸𝑅₁+𝑘₂𝑎₂=𝑅₂. (18)
Consequently, all solutions of model system (2) are bounded within a region Ω: Ω=𝑥₁(𝑡),𝑥₂(𝑡),𝑦(𝑡)∣0<𝑥₁(𝑡)+𝑥₂(𝑡)≤𝑅₁,0<𝑦(𝑡)≤𝑅₂, (19) where 𝑅₁=(𝑟₁+𝑑)²/4𝑏𝑑, 𝑅₂=(𝑟₂−𝑞𝐸)(𝑅₁+𝑘₂)/𝑎₂.

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 by 𝑃₀(0,0,0) and 𝑃₁(0,0,𝑘₂(𝑟₂−𝑞𝐸)/𝑎₂). The biological interpretations of 𝑃₀, 𝑃₁ are as follows. For 𝑃₀(0,0,0), it implies that all population in harvested prey predator ecosystem does not exist. For 𝑃₁(0,0,𝑘₂(𝑟₂−𝑞𝐸)/𝑎₂), 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 by 𝑃^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗). In order to discuss the existence of 𝑃^∗, it is equivalent to show that the following equations always have at least one positive solution: 𝑟₁𝑥₂−𝑑𝑥₁−𝑟₁𝑒^{−𝑑𝜏(𝑥1+𝑥2+𝑦)}𝑥₂𝑟=0,₁𝑒^{−𝑑𝜏(𝑥1+𝑥2+𝑦)}−𝑏𝑥₂−𝑎₁𝑦𝑥₂+𝑘₁𝑟=0,₂𝑎−𝑞𝐸−²𝑦𝑥₂+𝑘₂=0. (20)

It follows from (20) that 𝑟𝑦=²−𝑞𝐸𝑥₂+𝑘₂𝑎₂=𝑓𝑥₂,𝑟₁𝑥₂−𝑑𝑥₁=𝑟₁𝑥₂𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))},𝑟₁𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}=𝑏𝑥₂+𝑎₁ℎ𝑥₂, (21)(22)(23) where 𝑔(𝑥₂)=𝑥₂+𝑦=𝑥₂+𝑓(𝑥₂) and ℎ(𝑥₂)=𝑦/(𝑥₂+𝑘₁)=𝑓(𝑥₂)/(𝑥₂+𝑘₁).

Theorem 3 (existence of positive equilibrium). Supposing that hypotheses (H1)–(H4) hold, if 𝑘₁≥𝑘₂, 𝑟₂>𝑞𝐸, and 𝑎₂𝑟₁𝑒^−𝑑𝜏1>𝑎₁𝑟₂, then there exists at least one positive equilibrium 𝑃^∗.

Proof. Let Γ₁ and Γ₂ be the solution curves of (22) and (23) for 𝑥₁≥0, 𝑥₂≥0, respectively. The analytical properties of curve Γ₁ and Γ₂ are as follows.
For Γ₁: by simple computing, it can be found that (0,0)∈Γ₁.
According to (H2) and positiveness of all solutions of model system (2), it is easy to show that lim_𝑥1→+∞𝜏(𝑥₁(𝑡)+𝑔(𝑥₂))=𝜏₁, and lim_𝑥1→+∞𝑥₂𝑥₁=lim_𝑥1→+∞𝑑𝑥₁𝑟₁1−𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2)}=lim_𝑥1→+∞𝑑𝑥₁𝑟₁1−𝑒^−𝑑𝜏1=+∞. (24)
For Γ₂: by differentiating 𝑥₂ against 𝑥₁ along Γ₂, it can be obtained that 𝑑𝑥₂𝑑𝑥₁=−𝑑𝑟₁𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}×𝑑𝑟₁𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}𝜏^′𝑥₁+𝑔𝑥₂𝑎₂+𝑟₂−𝑞𝐸𝑎₂𝑎+𝑏+¹𝑟₂−𝑞𝐸𝑘₁−𝑘₂𝑎₂𝑥₂+𝑘₁²⁻¹. (25)
It can be shown that (𝑑𝑥₂/𝑑𝑥₁)<0, provided that 𝑘₁≥𝑘₂, 𝑟₂>𝑞𝐸, and then Γ₂ is strictly decreasing.
Furthermore, according to 𝑘₁≥𝑘₂, 𝑟₂>𝑞𝐸, and 𝑎₂𝑟₁𝑒^−𝑑𝜏1>𝑎₁(𝑟₂−𝑞𝐸), lim_𝑥1→∞𝑥₂𝑥₁=1𝑏𝑟₁𝑒^−𝑑𝜏1−𝑎₁𝑟₂−𝑞𝐸𝑥₂+𝑘₂𝑎₂𝑥₂+𝑘₁>0. (26)
Based on the above analysis, Γ₁ and Γ₂ intersect at some positive values, which proves the existence of positive equilibrium 𝑃^∗.

Theorem 4 (uniqueness of positive equilibrium). Supposing that hypotheses (H1)–(H4) hold, if the following inequality holds 1−𝑑𝜏^′𝑥^∗₁+𝑔𝑥^∗₂𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂×𝑎₁ℎ^′𝑥^∗₂+𝑑𝑔^′𝑥^∗₂𝜏^′𝑥^∗₁+𝑔𝑥^∗₂×𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂+𝑎₁𝑑ℎ^′𝑥^∗₂+𝑑𝜏^′(𝑥^∗₁+𝑔𝑥^∗₂𝑔^′𝑥^∗₂𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂×𝑑+𝑏+𝑑𝜏^′𝑥^∗₁+𝑔𝑥^∗₂>0, (27) then there exists a unique positive equilibrium.

Proof. Based on (22) and (23), 𝑥₂ can be defined as the function of 𝑥₁: Γ₁:𝑥₂=𝑔₁𝑥₁,Γ₂:𝑥₂=𝑔₂𝑥₁. (28) The positive equilibrium 𝑃^∗ will be unique, provided that 𝑔^′₁(𝑥₁)>𝑔^′₂(𝑥₁) for every such 𝑃^∗ otherwise reverse inequality holds.
By differentiating (22) with respect to 𝑥₁, it can be obtained that 𝑔^′₁𝑥₁=𝑑+𝑟₁𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}1−𝑑𝜏^′𝑥₁+𝑔𝑥₂𝑟₁1+𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}𝑑𝜏^′𝑥₁+𝑔𝑥₂𝑔^′𝑥₂. (29)
By differentiating (23) with respect to 𝑥₁, it can be obtained that 𝑔^′₂𝑥₁=−𝑏+𝑑𝑟₁𝜏^′𝑥₁+𝑔𝑥₂𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}𝑎₁ℎ^′𝑥₂+𝑑𝑟₁𝑔^′𝑥₂𝜏^′𝑥₁+𝑔𝑥₂𝑒^{−𝑑𝜏(𝑥1+𝑔(𝑥2))}. (30)
On the other hand, some expressions about positive equilibrium 𝑃^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗) can be obtained based on (22) and (23), 𝑟₁𝑥^∗₂−𝑑𝑥^∗₁=𝑟₁𝑥^∗₂𝑒^{−𝑑𝜏(𝑥∗1+𝑔(𝑥∗2))},𝑟₁𝑒^{−𝑑𝜏(𝑥∗1+𝑔(𝑥∗2))}=𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂. (31)
According to (31), 𝑔^′₁(𝑥^∗₁)>𝑔^′₁(𝑥^∗₂) is equivalent to the following inequality: 1−𝑑𝜏^′𝑥^∗₁+𝑔𝑥^∗₂𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂×𝑎₁ℎ^′𝑥^∗₂+𝑑𝑔^′𝑥^∗₂𝜏^′𝑥^∗₁+𝑔𝑥^∗₂×𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂+𝑎₁𝑑ℎ^′𝑥^∗₂+𝑑𝜏^′(𝑥^∗₁+𝑔𝑥^∗₂𝑔^′𝑥^∗₂×𝑏𝑥^∗₂+𝑎₁ℎ𝑥^∗₂𝑑+𝑏+𝑑𝜏^′𝑥^∗₁+𝑔𝑥^∗₂>0. (32)
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 𝑃=(𝑥₁,𝑥₂,𝑦) takes the following form: 𝜆+𝑑−𝐴𝑒^−𝑑𝜏(𝑧₎𝑟₁1+𝑒^{−(𝜆+𝑑)𝜏(}𝑧₎−𝐴𝑒^−𝑑𝜏(𝑧₎−𝐴𝑒^−𝑑𝜏(𝑧₎𝐴𝑒^−𝑑𝜏(𝑧₎𝑎𝜆+2𝑏𝑥+¹𝑘₁𝑦𝑥₂+𝑘₁²+𝐴𝑒^−𝑑𝜏(𝑧₎+𝑟₁𝑒^{−(𝜆+𝑑)𝜏(}𝑧₎𝐴𝑒^−𝑑𝜏(𝑧₎+𝑎₁𝑥₂𝑘₁+𝑥₂𝑎0−²𝑦²𝑘₂+𝑥₂²𝜆−𝑟₂+𝑞𝐸+2𝑎₂𝑦𝑘₂+𝑥₂=0, (33)where 𝐴=𝑑𝑟₁𝑥₂𝜏^′(𝑧).

Theorem 5. Local stability analysis of model system (2) around 𝑃₀ and 𝑃₁ is as follows: (a)if 𝑟₂<𝑞𝐸, then model system is locally stable around 𝑃₀, and 𝑃₁ is a saddle point which is unstable in the 𝑦-direction and stable in the 𝑥₁-𝑥₂ plane; (b)if 𝑟₂>𝑞𝐸, then model system is locally stable around 𝑃₁, and 𝑃₀ is a saddle point which is unstable in the 𝑦-direction and stable in the 𝑥₁-𝑥₂ plane.

Proof. For 𝑃₀(0,0,0), (33) reduces to 𝜆+𝑑𝑒^{−𝑑𝜏(0)}𝜆+𝑟₁𝑒^{−(𝜆+𝑑)𝜏(0)}𝜆−𝑟₂−𝑞𝐸=0. (34)
By solving (34), it can be found that there are two negative eigenvalues and only one positive eigenvalue, provided 𝑟₂>𝑞𝐸, which implies that 𝑃₀ is a saddle point which is unstable in the 𝑦-direction and stable in the 𝑥₁-𝑥₂ plane. On the other hand, there are three negative eigenvalues, provided that 𝑟₂<𝑞𝐸, which implies that 𝑃₀ is a stable point.
For 𝑃₁(0,0,(𝑘₂(𝑟₂−𝑞𝐸)/𝑎₂)), (33) reduces to 𝜆+𝑑𝑒^−𝑑𝜏(𝑦₎𝜆+𝑟₁𝑒^{−(𝜆+𝑑)𝜏(}𝑦₎+𝑎₁𝑦𝑘₁𝜆+𝑟₂−𝑞𝐸=0, (35) where 𝑦=𝑘₂(𝑟₂−𝑞𝐸)/𝑎₂. It follows from (35) that there are two negative eigenvalues and only one positive eigenvalue, provided 𝑟₂<𝑞𝐸, which implies that 𝑃₁ is a saddle point which is unstable in the 𝑦-direction and stable in the 𝑥₁-𝑥₂ plane. On the other hand, there are three negative eigenvalues, provided that 𝑟₂>𝑞𝐸, which implies that 𝑃₁ is a stable point.

In order to discuss the local stability of model system (2) around the positive equilibrium 𝑃^∗(𝑥^∗₁,𝑥^∗₂,𝑦^∗), (33) reduces to𝜆+𝑑𝐵^∗−𝐴^∗𝑟₁1+𝐵^∗𝑒^{−𝜆𝜏(𝑧∗)}−𝐴^∗−𝐴^∗𝐴^∗𝜆+𝐴^∗+𝐵^∗𝑒^{−𝜆𝜏(𝑧∗)}+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²𝐴^∗+𝑎₁𝑥^∗₂𝑥^∗₂+𝑘₁𝑎0−²𝑦^∗2𝑥^∗₂+𝑘₂²𝜆−𝑟₂+𝑞𝐸+2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂=0, (36)where 𝐴^∗=𝑑(𝑟₁𝑥^∗₂−𝑑𝑥^∗₁)𝜏^′(𝑧^∗), 𝐵^∗=(𝑟₁𝑥^∗₂−𝑑𝑥^∗₁)/𝑟₁𝑥^∗₂ and 𝑧^∗=𝑥^∗₁+𝑥^∗₂+𝑦^∗.

It can be computed that 𝑀(𝜆)+𝑁(𝜆)𝑒^{−𝜆𝜏(𝑧∗)}=0, (37) where 𝑀(𝜆)=𝜆³+𝑚₁𝜆²+𝑚₂𝜆+𝑚₃ and 𝑁(𝜆)=𝑛₁𝜆²+𝑛₂𝜆+𝑛₃, 𝑚₁=𝑑𝐵^∗+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²−𝑟₂+𝑞𝐸+2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂,𝑚₂=2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸𝐴^∗+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²−𝐴^∗𝑟₁−𝐴^∗,𝑚₃=𝑎₂𝑦^∗2𝐴^∗2𝑥^∗₂+𝑘₂²+𝑑𝐵^∗𝐴^∗+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²×2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸−𝐴^∗𝑟₁+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²×2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂,𝑛+𝑞𝐸₁=𝐵^∗,𝑛₂=𝐵^∗2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸−𝑟₁𝐴^∗,𝑛₃=𝐵^∗2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸𝑑𝐵^∗−𝐴^∗−𝑟₁𝐴^∗. (38)

In the following part, dynamical behavior of model system (2) around the positive equilibrium 𝑃^∗ is investigated. Furthermore, local stability analysis is discussed due to the variation of maturation delay and harvest effort level. By taking 𝜏^′(𝑧^∗) as a bifurcation parameter, conditions for local stability switch are discussed with the increase of 𝜏^′(𝑧^∗) from zero.

Case 1 (𝜏^′(𝑧^∗)=0). In the case of 𝜏^′(𝑧^∗)=0, it derives that 𝜏(𝑧^∗) remains as a constant (zero or a positive constant) for all time 𝑡>0 based on (H2). In the following part, 𝜏(𝑧^∗) is denoted as 𝜏^∗ for simplifying. Furthermore, it can be computed that 𝐴^∗=0, and 𝑚_𝑖,𝑛_𝑖 (𝑖=1,2,3) in (37) can be rewritten as follows: 𝑚₁=𝑑𝐵^∗+2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²−𝑟₂+𝑞𝐸+2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂,𝑚₂=2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²,𝑚₃=𝑑𝐵^∗2𝑏𝑥^∗₂+𝑎₁𝑘₁𝑦^∗𝑥^∗₂+𝑘₁²2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂,𝑛+𝑞𝐸₁=𝐵^∗,𝑛₂=𝐵^∗2𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂,𝑛+𝑞𝐸₃=𝑑𝐵^∗22𝑎₂𝑦^∗𝑥^∗₂+𝑘₂−𝑟₂+𝑞𝐸. (39)

Theorem 6. Supposing that hypotheses (H1)–(H4) hold, if 𝑟₂−𝑞𝐸>0, then model system (2) is stable around the positive equilibrium 𝑃^∗ in the case of 𝜏^∗=0.

Proof. When 𝜏^∗=0, (37) can be rewritten as follows: 𝜆³+𝑚₁+𝑛₁𝜆²+𝑚₂+𝑛₂𝜆+𝑚₃+𝑛₃=0. (40)
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, 𝑃^∗ is locally stable in the case of 𝜏^∗=0.

When 𝜏^∗>0, let 𝜆=𝑖𝜔 be a root of (37), where 𝜔 is positive. Substitute 𝜆=𝑖𝜔 into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows: 𝜔³−𝑚₂𝜔=𝑛₁𝜔²−𝑛₃sin𝜔𝜏^∗+𝑛₂𝜔cos𝜔𝜏^∗,𝑚₁𝜔²−𝑚₃=𝑛₂𝜔sin𝜔𝜏^∗−𝑛₁𝜔²−𝑛₃cos𝜔𝜏^∗. (41)

By squaring and adding (41), it can be obtained that 𝜔⁶+𝐵₁𝜔⁴+𝐵₂𝜔²+𝐵₃=0, (42) where 𝐵₁=𝑚²₁−2𝑚₂−𝑛²₁, 𝐵₂=𝑚²₂−2𝑚₁𝑚₃+2𝑛₁𝑛₃−𝑛²₂, 𝐵₃=𝑚²₃−𝑛²₃, and 𝑚_𝑖, 𝑛_𝑖 (𝑖=1,2,3) have been defined in (40).

According to the values of 𝐵_𝑖 (𝑖=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is 𝐵₃<0.

If 𝐵₃<0 holds, then (42) has a positive root 𝜔₀, and (37) has a pair of purely imaginary roots of the form ±𝑖𝜔₀. Consequently, it can be obtained by eliminating sin(𝜔𝜏^∗) from (41): cos𝜔𝜏^∗=𝑛₂−𝑚₁𝑛₁𝜔⁴+𝑚₁𝑛₃+𝑚₃𝑛₁−𝑚₂𝑛₂𝜔²−𝑚₃𝑛₃𝑛₂𝜔²²+𝑛₃−𝑛₁𝜔²², (43) The 𝜏_𝑘 corresponding to 𝜔₀ is as follows: 𝜏_𝑘=1𝜔₀arccos𝑛₂−𝑚₁𝑛₁𝜔⁴+𝑚₁𝑛₃+𝑚₃𝑛₁−𝑚₂𝑛₂𝜔²−𝑚₃𝑛₃×𝑛₂𝜔²²+𝑛₃−𝑛₁𝜔²²⁻¹+2𝑘𝜋𝜔₀, (44)𝑘=0,1,2,…. By virtue of Butler’s lemma [37], it can be concluded that the positive equilibrium 𝑃^∗ remains locally stable for 𝜏^∗<𝜏₀, as 𝑘=0.

Case 2 (𝜏^′(𝑧^∗)>0). In the case of 𝜏^′(𝑧^∗)>0, local stability of model system (2) around the positive equilibrium 𝑃^∗ can change only if there exists at least one root of (37) such that Re𝜆=0.
Let 𝜆=𝑖𝜈 be one such root, where 𝜈 is positive. Substitute 𝜆=𝑖𝜈 into (37) and separate the real and imaginary parts, and then two transcendental equations can be obtained as follows: 𝜈³−𝑚₂𝜈=𝑛₁𝜈²−𝑛₃sin𝜈𝜏𝑧^∗+𝑛₂𝜈cos𝜈𝜏𝑧^∗,𝑚₁𝜈²−𝑚₃=𝑛₂𝜈sin𝜈𝜏𝑧^∗−𝑛₁𝜈²−𝑛₃cos𝜈𝜏𝑧^∗. (45)
By squaring and adding (45), it can be obtained that 𝜈⁶+𝐶₁𝜈⁴+𝐶₂𝜈²+𝐶₃=0, (46) where 𝐶₁=𝑚²₁−2𝑚₂−𝑛²₁,𝐶₂=𝑚²₂−2𝑚₁𝑚₃+2𝑛₁𝑛₃−𝑛²₂, 𝐶₃=𝑚²₃−𝑛²₃, and 𝑚_𝑖, 𝑛_𝑖 (𝑖=1,2,3) have been defined in (37).
According to the values of 𝐶_𝑖 (𝑖=1,2,3) and the Routh-Hurwitz criteria [1], a simple assumption of the existence of a positive root for (42) is 𝐶₃<0, which derives that 𝜏^′𝑧^∗>𝑎₁𝑘₁𝑟₁𝑟₂−𝑞𝐸𝑥^∗₂𝑥^∗₂+𝑘₂+𝑎₂𝑥^∗₂+𝑘₁²𝑟₁𝑥^∗₂2𝑏𝑥^∗₂−1+𝑑𝑥^∗₁×𝑟₁𝑥^∗₂𝑎₁𝑘₁𝑟₁𝑟₂−𝑞𝐸𝑥^∗₂𝑥^∗₂+𝑘₂+𝑎₂𝑥^∗₂+𝑘₁²×𝑑𝑥^∗₁(𝑑+2)+𝑟₁𝑥^∗₂2𝑏𝑥^∗₂−𝑑−1⁻¹. (47)
If the above inequality holds, then model system (2) is unstable around the positive equilibrium 𝑃^∗ in the case of 𝜏^′(𝑧^∗)>0.

3.4. Global Stability Analysis of Positive Equilibrium

In this section, global stability of the positive equilibrium 𝑃^∗ is discussed by using an iterative technique in the case of 𝜏^′(𝑧^∗)=0.

Lemma 7 (see [29]). Consider the following equation: 𝑥̇=𝑎𝑥(𝑡−𝜏)−𝑏𝑥(𝑡)−𝑐𝑥²(𝑡), (48) where 𝑎, 𝑏, 𝑐, and 𝜏 are positive constants, and 𝑥(𝑡)>0 for 𝑡∈[−𝜏,0]; it follows that (i)If 𝑎>𝑏, then lim_𝑡→+∞𝑥(𝑡)=(𝑎−𝑏)/𝑐; (ii)If 𝑎<𝑏, then lim_𝑡→+∞𝑥(𝑡)=0.

Theorem 8. Supposing that hypotheses (H1)–(H4) and 𝑟₂−𝑞𝐸>0 hold, if the following inequalities hold 𝑏𝑘₂+𝑟₁𝑒^{−𝑑𝜏∗}<𝑎₂𝑏𝑦^∗,𝑎₂𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}>𝑎₁𝑟₂−𝑞𝐸𝑅₁+𝑘₂, (49) then the positive equilibrium 𝑃^∗ is globally asymptotically stable in the case of 𝜏^′(𝑧^∗)=0.

Proof. In the case of 𝜏^′(𝑧^∗)=0, it derives that 𝜏(𝑧^∗) remains as a constant (zero or a positive constant) for all time 𝑡>0 based on (H2). In the following part, 𝜏(𝑧^∗) is denoted as 𝜏^∗ for simplifying. Let 𝑈₁=lim_𝑡→+∞sup𝑥₂(𝑡),𝑉₁=lim_𝑡→+∞inf𝑥₂(𝑡),𝑈₂=lim_𝑡→+∞sup𝑦(𝑡),𝑉₂=lim_𝑡→+∞inf𝑦(𝑡). (50)
In the following, we will claim that 𝑈₁=𝑉₁=𝑥^∗₂, 𝑈₂=𝑉₂=𝑦^∗.
It follows from Theorem 2 that 𝑥₂(𝑡)≤𝑅₁ (𝑅₁ has been defined in Theorem 2). From model system (2), 𝑦̇(𝑡)≤𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝑅₁+𝑘₂. (51)
By standard comparison argument, it derives that 𝑈₂≤𝑟₂−𝑞𝐸𝑅₁+𝑘₂𝑎₂∶=𝐽^𝑦₁, (52) and then for sufficiently small 𝜖>0, there exists a 𝑇₁₁>0 such that if 𝑡>𝑇₁₁, 𝑦(𝑡)≤𝐽^𝑦₁+𝜖. Based on Theorem 1, 𝑥₂(𝑡)+𝑘₁>𝑘₁, it can be obtained that for 𝑡>𝑇₁₁+𝜏^∗, 𝑥̇₂(𝑡)≥𝑟₁𝑒^{−𝑑𝜏∗}𝑥₂𝑡−𝜏^∗−𝑏𝑥²₂(𝑡)−𝑎₁𝐽^𝑦₁+𝜖𝑥₂(𝑡)𝑥₂(𝑡)+𝑘₁>𝑟₁𝑒^{−𝑑𝜏∗}𝑥₂𝑡−𝜏^∗−𝑏𝑥²₂(𝑡)−𝑎₁𝐽^𝑦₁+𝜖𝑥₂(𝑡)𝑘₁. (53)
Consider the following auxiliary equation: 𝑣̇(𝑡)=𝑟₁𝑒^{−𝑑𝜏∗}𝑣𝑡−𝜏^∗−𝑏𝑣²(𝑡)−𝑎₁𝐽^𝑦₁+𝜖𝑣(𝑡)𝑘₁. (54)
Under the condition 𝑎₂𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}>𝑎₁(𝑟₂−𝑞𝐸)(𝑅₁+𝑘₂), it follows from Lemma 7 that lim_𝑡→+∞𝑣(𝑡)=𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}−𝑎₁𝐽^𝑦₁+𝜖𝑏𝑘₁∶=𝐼^𝑥₁. (55)
Hence, 𝑉₁≥𝐼^𝑥₁. For sufficiently small 𝜖>0, there exits 𝑇₁₂≥𝑇₁₁+𝜏^∗ such that if 𝑡>𝑇₂₂, then 𝑥₂(𝑡)≥𝐼^𝑥₁−𝜖.
We derive from the model system (2) that for 𝑡>𝑇₁₂, 𝑦̇(𝑡)≥𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝑘₂+𝐼^𝑥₁.−𝜖 (56)
A standard comparison argument shows that lim_𝑡→+∞𝑦(𝑡)=𝑟₂−𝑞𝐸𝑘₂+𝐼^𝑥₁−𝜖𝑎₂∶=𝐼^𝑦₁. (57)
Hence, for sufficiently small 𝜖>0, there is a 𝑇₂₁≥𝑇₁₂ satisfying if 𝑡>𝑇₂₁, then 𝑦(𝑡)≥𝐼^𝑦₁−𝜖. Consequently, for 𝑡>𝑇₂₁+𝜏^∗, 𝑥̇₂(𝑡)≤𝑟₁𝑒^{−𝑑𝜏∗}𝑥₂𝑡−𝜏^∗−𝑏𝑥²₂(𝑡)−𝑎₁𝐼^𝑦₁−𝜖𝑥₂(𝑡)𝑅₁+𝑘₁. (58)
Consider the following auxiliary equation: 𝑣̇(𝑡)=𝑟₁𝑒^{−𝑑𝜏∗}𝑣𝑡−𝜏^∗−𝑏𝑣²(𝑡)−𝑎₁𝐼^𝑦₁−𝜖𝑣(𝑡)𝑅₁+𝑘₁. (59)
It follows from Lemma 7 that lim_𝑡→+∞𝑣(𝑡)=𝑟₁𝑒^{−𝑑𝜏∗}−𝑏𝑅₁+𝑘₁𝑎₁𝐼^𝑦₁−𝜖∶=𝐽^𝑥₁. (60)
Hence, 𝑈₁≤𝐽^𝑥₁. For sufficiently small 𝜖>0, there exists a 𝑇₂₂≥𝑇₂₁+𝜏^∗ satisfying that if 𝑡>𝑇₂₂, then 𝑥₂(𝑡)≤𝐽^𝑥₁+𝜖. For 𝑡>𝑇₂₂, it gives that 𝑦̇(𝑡)≤𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝑘₂+𝐽^𝑥₁.+𝜖 (61)
By standard comparison argument, it derives that lim_𝑡→+∞𝑦(𝑡)=𝑟₂−𝑞𝐸𝐽^𝑥₁+𝑘₂+𝜖𝑎₂∶=𝐽^𝑦₂. (62)
Hence, for sufficiently small 𝜖>0, there exists 𝑇₃₁≥𝑇₂₂ satisfying that if 𝑡>𝑇₃₁, 𝑦(𝑡)≤𝐽^𝑦₂+𝜖, the for 𝑡>𝑇₃₁+𝜏^∗𝑥̇₂(𝑡)≥𝑟₁𝑒^{−𝑑𝜏∗}𝑥₂𝑡−𝜏^∗−𝑏𝑥²₂(𝑡)−𝑎₁𝐽^𝑦₂+𝜖𝑥₂(𝑡)𝑘₁. (63)
Consider the following auxiliary equation: 𝑣̇(𝑡)=𝑟₁𝑒^{−𝑑𝜏∗}𝑣𝑡−𝜏^∗−𝑏𝑣²(𝑡)−𝑎₁𝐽^𝑦₂+𝜖𝑣(𝑡)𝑘₁. (64)
By using Lemma 7, it can be obtained that lim_𝑡→+∞𝑣(𝑡)=𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}−𝑎₁𝐽^𝑦₂+𝜖𝑏𝑘₁∶=𝐼^𝑥₂. (65)
Since it is true for any sufficiently small 𝜖>0, 𝑉₁≥𝐼^𝑥₂. Therefore, there exists 𝑇₃₂≥𝑇₃₁+𝜏^∗ such that if 𝑡>𝑇₃₂, then 𝑥₂(𝑡)≥𝐼^𝑥₂−𝜖.
It follows from model system (2) that for 𝑡>𝑇₃₂, 𝑦̇(𝑡)≥𝑟₂𝑦(𝑡)−𝑎−𝑞𝐸²𝑦²(𝑡)𝐼^𝑥₂−𝜖+𝑘₂. (66)
By using standard comparison argument, it derives that lim_𝑡→+∞𝑦(𝑡)=𝑟₂−𝑞𝐸𝐼^𝑥₂−𝜖+𝑘₂𝑎₂∶=𝐼^𝑦₂. (67)
Since this is true for any sufficiently small 𝜖>0, 𝑉₂≥𝐼^𝑦₂. Consequently, there exists 𝑇₄₁≥𝑇₃₂ satisfying if 𝑡>𝑇₄₁, then 𝑦(𝑡)≥𝐼^𝑦₂−𝜖.
It follows from model system (2) that for 𝑡>𝑇₄₁+𝜏^∗, 𝑥̇₂(𝑡)≤𝑟₁𝑒^{−𝑑𝜏∗}𝑥₂𝑡−𝜏^∗−𝑏𝑥²₂(𝑡)−𝑎₁𝐼^𝑦₂−𝜖𝑥₂(𝑡)𝑅₁+𝑘₁. (68)
Consider the following auxiliary equation, 𝑣̇(𝑡)≤𝑟₁𝑒^{−𝑑𝜏∗}𝑣𝑡−𝜏^∗−𝑏𝑣²(𝑡)−𝑎₁𝐼^𝑦₂−𝜖𝑣(𝑡)𝑅₁+𝑘₁. (69)
By using Lemma 7, it derives that lim_𝑡→+∞𝑥₂(𝑡)=𝑟₁𝑒^{−𝑑𝜏∗}𝑅₁+𝑘₁−𝑎₁𝐼^𝑦₂−𝜖𝑏𝑅₁+𝑘₁∶=𝐽^𝑥₂. (70)
Continuing the above process, four sequences {𝐼^𝑥_𝑛}, {𝐼^𝑦_𝑛}, {𝐽^𝑥_𝑛}, {𝐽^𝑦_𝑛}, 𝑛=1,2,…, are obtained which take the following form 𝐽^𝑥_𝑛=𝑟₁𝑒^{−𝑑𝜏∗}𝑅₁+𝑘₁+𝑎₁𝜖−𝑎₁𝐼^𝑦_𝑛𝑏𝑅₁+𝑘₁,𝐽^𝑦_𝑛=𝜖+𝑘₂𝑟₂−q𝐸+𝑟₂−𝑞𝐸𝐽^𝑥_𝑛−1𝑎₂,𝐼^𝑥_𝑛=𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}−𝑎₁𝜖−𝑎₁𝐽^𝑦_𝑛𝑏𝑘₁,𝐼^𝑦_𝑛=𝑟₂−𝑞𝐸𝑘₂−𝜖+𝑟₂−𝑞𝐸𝐼^𝑥_𝑛𝑎₂. (71)
It is easy to show that 𝐼^𝑥_𝑛≤𝑉₁≤𝑈₁≤𝐽^𝑥_𝑛,𝐼^𝑦_𝑛≤𝑉₂≤𝑈₂≤𝐽^𝑦_𝑛. (72)
By virtue of (71), it derives that 𝐽^𝑦_𝑛=𝑟₂−𝑞𝐸𝑏𝜖+𝑘₂+𝑟₁𝑒^{−𝑑𝜏∗}𝑎₂𝑏+𝑎₁𝜖𝑟₂−𝑞𝐸𝑎₂𝑏𝑅₁+𝑘₁−𝑎₁𝑟₂−𝑞𝐸²𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}+𝑏𝑘₂+𝜖𝑏𝑘₁−𝑎₁−𝑎₁𝐽^𝑦_𝑛−1𝑎²₂𝑏²𝑘₁𝑅₁+𝑘₁. (73)
Furthermore, 𝐽^𝑦_𝑛−𝐽^𝑦_𝑛−1=𝑟₂−𝑞𝐸𝑏𝜖+𝑘₂+𝑟₁𝑒^{−𝑑𝜏∗}𝑎₂𝑏+𝑎₁𝜖𝑟₂−𝑞𝐸𝑎₂𝑏𝑅₁+𝑘₁−𝑎₁𝑟₂−𝑞𝐸²𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}+𝑏𝑘₂+𝜖𝑏𝑘₁−𝑎₁−𝑎₁𝐽^𝑦_𝑛−1𝑎²₂𝑏²𝑘₁𝑅₁+𝑘₁−𝐽^𝑦_𝑛−1. (74)
If the following inequalities hold 𝑏𝑘₂+𝑟₁𝑒^{−𝑑𝜏∗}<𝑎₂𝑏𝑦^∗,𝑎₂𝑘₁𝑟₁𝑒^{−𝑑𝜏∗}>𝑎₁𝑟₂−𝑞𝐸𝑅₁+𝑘₂, (75) then 𝐽^𝑦_𝑛−𝐽^𝑦_𝑛−1≤0, which implies that {𝐽^𝑦_𝑛∣𝐽^𝑦_𝑛≥𝑦^∗,𝑛=1,2,…} is monotonically decreasing. Hence, it can be shown that limitation of sequence {𝐽^𝑦_𝑛} exists. Taking 𝑛→+∞, it follows from (73) that lim_𝑛→+∞𝐽^𝑦_𝑛=𝑦^∗. (76)
By using (71) and (76), it can be shown that lim_𝑛→+∞𝐽^𝑥_𝑛=𝑥^∗₂,lim_𝑛→+∞𝐼^𝑦_𝑛=𝑦^∗,lim_𝑛→+∞𝐼^𝑥_𝑛=𝑥^∗₂. (77)
According to the definition of 𝑈₁, 𝑈₂, 𝑉₁, and 𝑉₂, it derives that 𝑈₁=𝑉₁=𝑥^∗₂,𝑈₂=𝑉₂=𝑦^∗. (78)
Hence, lim_𝑡→+∞𝑥₂(𝑡)=𝑥^∗₂,lim_𝑡→+∞𝑦(𝑡)=𝑦^∗. (79)
Based on (5), it derives that 𝑥₁(𝑡)=^𝑡_𝑡−𝜏𝑟₁𝑒^{−𝑑(𝑡−𝑠)}𝜙(𝑠)𝑑𝑠. (80)
By using L’Hospital’s rule, it derives that lim_𝑡→+∞𝑥₁(𝑡)=lim_𝑡→+∞𝑟₁𝑒^𝑑𝑡𝜙(𝑡)−𝑒^𝑑(𝑡−𝜏₎𝜙(𝑡−𝜏)𝑑𝑒^𝑑𝑡=lim_𝑡→+∞𝑟₁𝑑𝜙(𝑡)−𝑒^−𝑑𝜏𝜙(𝑡−𝜏)=lim_𝑡→+∞𝑟₁𝑑1−𝑒^−𝑑𝜏𝑥^∗₂. (81)
According to (20), it is easy to show that lim_𝑡→+∞𝑥₁(𝑡)=lim_𝑡→+∞𝑟₁𝑑1−𝑒^−𝑑𝜏𝑥^∗₂=𝑥^∗₁. (82)
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 𝜏(𝑧(𝑡)) takes the following form [23]: 𝜏(𝑧(𝑡))=𝜏₀+𝜏_𝑚−𝜏_𝑚𝑒^{−𝑧(𝑡)}, (83) where 𝜏_𝑚∈(𝜏₀,𝜏₁) satisfying 𝜏₀+𝜏_𝑚=𝜏₁. Based on Theorem 1, it follows from simple computation that lim_𝑡→0𝜏(𝑧(𝑡))=𝜏₀,lim_𝑡→+∞𝜏(𝑧(𝑡))=𝜏₁, (84) 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. 𝑟₁=12, 𝑑=0.2, 𝑏=1.2, 𝑎₁=0.5, 𝑎₂=2, 𝑘₁=2, 𝑘₂=1, 𝑟₂=2, 𝑞=0.25, and 𝐸=4. According to the given values of parameters, it follows from Theorems 3 and 4 that there exists a unique positive equilibrium 𝑃^∗(141.5454,5.6252,3.3126). Furthermore, it can be verified that 𝑃^∗ 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 Theorems 1 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.

Acknowledgments

The authors gratefully acknowledge editors and anonymous reviewers' comments and patient work. This work is supported by National Natural Science Foundation of China, grant No. 61104003, grant No. 61273008 and grant No. 61104093. Research Foundation for Doctoral Program of Higher Education of Education Ministry, grant No. 20110042120016. Hebei Province Natural Science Foundation, grant No. F2011501023. Fundamental Research Funds for the Central Universities, grant No. N120423009. Research Foundation for Science and Technology Pillar Program of Northeastern University at Qinhuangdao, grant No. XNK201301. This work is supported by State Key Laboratory of Integrated Automation of Process Industry, Northeastern University, supported by Hong Kong Admission Scheme for Mainland Talents and Professionals, Hong Kong Special Administrative Region.