Dynamics of Stage-Structured Predator–Prey Model with Beddington–DeAngelis Functional Response and Harvesting

In this paper, we investigate the stability of equilibrium in the stage-structured and density-dependent predator–prey system with Beddington–DeAngelis functional response. First, by checking the sign of the real part for eigenvalue, local stability of origin equilibrium and boundary equilibrium are studied. Second, we explore the local stability of the positive equilibrium for τ=0 and τ≠0 (time delay τ is the time taken from immaturity to maturity predator), which shows that local stability of the positive equilibrium is dependent on parameter τ. Third, we qualitatively analyze global asymptotical stability of the positive equilibrium. Based on stability theory of periodic solutions, global asymptotical stability of the positive equilibrium is obtained when τ=0; by constructing Lyapunov functions, we conclude that the positive equilibrium is also globally asymptotically stable when τ≠0. Finally, examples with numerical simulations are given to illustrate the obtained results.


Introduction
The dynamical behavior of the predator-prey system is one of the main research topics in mathematical ecology and theoretical biology [1][2][3][4][5][6][7][8][9].Functional response is the core component of the community and food web model, and their mathematical form strongly affects the dynamics and stability of the ecosystem [10][11][12][13][14]. Beddington [15] and DeAngelis [16] originally proposed the predator-prey model as follows: m 1 +m 2 x(t)+m 3 y(t) , dy(t) dt = y(t) −d + f x(t) m 1 +m 2 x(t)+m 3 y(t) , (1) where x(t) is the prey density, y(t) is the predator density, a is the intrinsic growth rate of the prey, b represents the intensity of intraspecific competition of the prey, and d represents the predator's death rate.The Beddington-DeAngelis functional response is similar to the well-known Holling II type [17,18] functional response, but it has an extra term m 3 y in the denominator modeling mutual interference among predators.
The functional response of consumers is a function of resource density [17,19,20].However, it has been shown that other species and other predators can alter the predation process directly or indirectly [21].In addition, Kratina [21] showed that predator dependence is very important not only when the predator density on per capita predation rate is very high, but also when the predator density is low.Therefore, we need to consider the realistic level of predator density when we study the predator-prey system.
Moreover, in nature, many species undergo two stages [22,23]: immature and mature, and species at these two stages may have different behaviors.The model of single-species stage-structured dynamics [22] was described as where x 1 (t) and x 2 (t) represent the immature and mature populations densities, respectively, and τ represents a constant time to maturity.Therefore, in order to be in accord with the natural phenomenon, the system with the stage structure recently has been extensively studied [24][25][26][27][28][29][30].Liu and Beretta [27] studied time delay τ in the response term x(t) m 1 +m 2 x(t)+m 3 y(t) of (1) in the predator equation, that is, 1+k 1 x(t)+k 2 y(t) , dy(t) dt = −dy(t) + nge −d j τ x(t−τ)y(t−τ) 1+k 1 x(t−τ)+k 2 y(t−τ) , dy i (t) dt = ngx(t)y(t) 1+k 1 x(t)+k 2 y(t) − nge −d j τ x(t−τ)y(t−τ) 1+k 1 x(t−τ)+k 2 y(t−τ) − d j y j (t), where g > 0 (units: 1/time) and k 1 > 0 (units: 1/prey) stand the effects of capture rate and handling time, respectively, on the feeding rate; n is the birth rate of the predator; and k 2 ≥ 0 (units: 1/predator) stands the magnitude of interference among predators.Liu and Beretta [27] pointed out the difference between Beddington-DeAngelis functional response and Holling type II, and the effect of k 2 y (describing mutual interference by predators) on the dynamic of the system (3).On the basis of the system (1) and the system (2), She and Li [29] investigated a predator-prey system with density-dependence for predator and stage structure for prey where p stands for predator density-dependent mortality rate, the predator consumes prey with functional response of Beddington-DeAngelis type cx(t)y(t) m 1 +m 2 x(t)+m 3 y(t) and contributes to its growth with rate f x(t)y(t) m 1 +m 2 x(t)+m 3 y(t) .Note that compared with the system (1), the system (4) contains not only bx 2 (t) (which stands for intraspecific competition of prey species), but also py 2 (t) (which stands for intraspecific competition of predator species).That is, they consider both the prey density dependence and the predator density dependence in the predator-prey model (4).She and Li [29] studied the dynamics of the system (4) and pointed out the impact of the predator density-dependent mortality rate p on the global attraction and permanence of the system (4).
The development of biological resources, the management of renewable resources, and the harvest of populations are universal human purposes for realizing the economic benefits of fishery, forestry and wildlife management [31,32].Many researchers [28,30,33,34] have extensively studied the predator-prey model with harvesting and the role of harvesting in renewable resource management.Brauer [35] introduced the predator-prey system with constant-rate prey harvesting where prey is harvested at a constant time rate F. May [36] put forward two types of harvesting regimes: constant-yield harvesting (representing harvested biomass independent of the size of the population) and constant-effort harvesting (representing harvested biomass proportional to the size of the population).
Based on the system (3), we construct a density-dependent and constant-effort harvesting predator-prey model: where y i (t) denotes the immature or juvenile predator density, juveniles suffer a mortality rate d i and take τ units of time to be mature and e −d i τ is the surviving rate of each immature predator to reach maturity.E 1 and E 2 denote the harvesting effort of the mature population of prey and predator, respectively.Further, all the parameters a, d i , b, c, d, f , m 1 , m 2 , m 3 , p, E 1 , E 2 , and τ are positive.
In the system (6), as y i (t) does not intervene in the dynamics of x(t) and y(t), system ( 6) is equal to the following system: The initial conditions of the system (7) is where This paper mainly investigates the local and global stability of positive equilibrium in the system (7) on parameter τ, which is organized as follows.In Section 2, we study local stability of origin equilibrium and boundary equilibrium.In Section 3, we derive local stability of the positive equilibrium for τ = 0 and τ = 0, respectively.In Section 4, we obtain the global asymptotical stability of the positive equilibrium for τ = 0 and τ = 0, respectively.Last, we conclude the paper.

Local Stability of Origin Equilibrium and Boundary Equilibrium
For any value of all parameters, system (7) has the equilibria E 0 (0, 0) and E 1 a−E 1 b , 0 , denoted as the origin and the boundary equilibrium, respectively.In the following, we determine the local stability of two equilibria by the sign of eigenvalue for the corresponding characteristic matrix.
First, for origin equilibrium E 0 , the corresponding characteristic matrix is and the eigenvalues are Clearly, E 0 (0, 0) is hyperbolic saddle and is unstable.
Next, for boundary equilibrium E 1 , the corresponding characteristic matrix is We can obtain one eigenvalue λ 1 = −(a − E 1 ) < 0, the second eigenvalue is determined by the following equation: Let n be the independent variable and g be the dependent variable; then, the straight line g = n + (d + E 2 ) and the curve g = (a−E 1 ) f e −d i τ bm 1 +m 2 (a−E 1 ) e −nτ must intersect at a unique point (λ 2 , g).Further, we can get the following: Therefore, we have the following conclusion about local stability of boundary equilibrium.
b , 0 is locally asymptotically stable; when τ < τ * , a−E 1 b , 0 is unstable; when τ = τ * , a−E 1 b , 0 is linearly neutrally stable, where That is, τ * is a threshold for the stability of the boundary equilibrium.

Local Stability of the Positive Equilibrium
We denote (x * , y * ) as equilibrium other than the origin equilibrium and boundary equilibrium for the system (7) and satisfying the algebraic equations By considering the intersection of curves y = F(x) and x = H(y) in the first quadrant, we can obtain if the condition holds, (x * , y * ) is a positive equilibrium of system (7).By the condition (10), we can directly obtain the following result.
The parameter τ is the time from immature to mature predator.Time delay τ = 0, it shows the predator population is divided into immature and mature; time delay τ = 0, it shows the model does not consider the immature predators.In the following, we discuss the local stability of positive equilibrium (x * , y * ) for τ = 0 and τ = 0, respectively.
Theorem 2. Let (x * , y * ) is a positive equilibrium of system (7), if the inequality (10) and the following inequalities hold, the positive equilibrium (x * , y * ) of system ( 7) is locally asymptotically stable for τ = 0.

Local Stability of the Positive Equilibrium for τ = 0
In this subsection, we further discuss the sign of the real part of the characteristic root for the characteristic Equation (11) with the change of parameter τ, and determine stability switch.
First, assume that characteristic Equation ( 11) has a characteristic root with zero real part, and let it be λ = iω(ω ∈ R and ω = 0).Putting it into the characteristic Equation ( 11), we have Here, a 1 (τ), a 2 (τ), a 3 (τ) and a 4 (τ) are abbreviated as a 1 , a 2 , a 3 and a 4 , respectively.From Equation ( 15), we can get Regarding ω 2 as an invariant, by the quadratic root formula, we have We discuss the case of roots for the Equation ( 16): ω 2 + and ω 2 − are negative, which obviously contradicts with ω ∈ R. Therefore, Equation ( 16) does not have real roots, that is, characteristic Equation (11) does not have purely imaginary roots.Moreover, if the conditions (10) and ( 13) hold, all roots of the characteristic Equation ( 12) have negative real part for τ = 0. Therefore, by Rouche's theorem, it follows that the roots of the characteristic Equation (11) also have negative real part. ( Equation ( 16) has a positive root ω 2 0 .Furthermore, putting ω 2 0 into Equations ( 15), we have Equation ( 16) has two positive roots ω 2 ± .Putting ω 2 ± into Equations ( 15), we obtain From the above analysis, we have the following conclusion.
Second, we have that when τ > τ 0 , τ > τ + j , and τ < τ − j , there are roots with positive real part in the characteristic Equation (11).Let these characteristic roots with positive real part be According to the characteristic Equation ( 11), we get Obviously, sign d dτ Reλ(τ) = sign Re( dλ(τ) dτ ) −1 .Therefore, in order to judge the sign of d dτ Reλ(τ), we only need to judge the sign of Re dλ(τ) dτ −1 . According to (23), we can get Last, we can easily prove that λ 0 and τ ± j satisfy the following transversality conditions: It follows that λ 0 and τ ± j are bifurcation values.Thus, we have the following theorem about the distribution of the characteristic roots of Equation (11).
(iii) If conditions ( 10) and ( 21) hold, there is a positive integer k such that the positive equilibrium (x * , y * ) of system (7) has k stability switches from stability to instability to stability.That is, when ), all roots of Equation ( 11) have negative real parts, and the positive equilibrium (x * , y * ) of system ( 7) is locally asymptotically stable.When Equation ( 11) has at least one root with positive real part, the positive equilibrium (x * , y * ) of system ( 7) is unstable.Furthermore, it shows when time delay τ passes the critical value τ + j , j = 0, 1, 2, ..., k − 1, the positive equilibrium (x * , y * ) of the system (7) loses its stability and Hopf bifurcation occurs.

Global Asymptotic Stability of Positive Equilibrium
In order to maintain healthy and sustainable ecosystem, the relevant departments need to plan harvesting.From the perspective of ecological managers, the positive equilibrium is best to be globally asymptotically stable.Therefore, we will investigate the global asymptotic stability of the positive equilibrium (x * , y * ) for τ = 0 and τ = 0, respectively.

Global Asymptotic Stability of Positive Equilibrium for τ = 0
In this subsection, when τ = 0, our aim is to show that the positive equilibrium (x * , y * ) is globally asymptotically stable.
When τ = 0, then system (7) becomes The variational matrix of the system ( 25) is given by and the trace of variational matrix is Let Γ(t) = (x(t), y(t)) be any nontrivial periodic orbit of system (25) with period T > 0. Next, following the proof of Theorem 2.1 in [38], we prove T 0 tr(J)dt < 0 and obtain the positive equilibrium (x * , y * ) of system ( 7) is globally asymptotically stable. As Obviously, if cm 2 ≤ f m 3 , we have tr(J) < 0 for all t ≥ 0. Therefore, T 0 tr(J)dt < 0. By Lemma 3.1 (P224) in [39], we can obtain that any non-trivial periodic orbit Γ(t) = (x(t), y(t)) with period T(T > 0) is orbitally stable.However, it is impossible because the equilibrium (x * , y * ) is inside the periodic orbit Γ(t), Γ(t) is orbitally stable, and (x * , y * ) is locally asymptotically stable, there must exist an unstable periodic orbit between (x * , y * ) and Γ(t).This leads to a contradiction, and the assumption of nontrivial periodic orbit Γ(t) is not true.Therefore, the system (25) does not have periodic orbits in R 2 + .In the system (25), (x * , y * ) is locally asymptotically stable, equilibria (0, 0) and a−E 1 b , 0 are all hyperbolic saddles, and there is not other attractors, then the other trajectories tend to (x * , y * ).Therefore, (x * , y * ) is globally asymptotically stable.
Example 2 (Example 1 continued).For system (14), the positive equilibrium (1.2040, 0.2505) is globally asymptotically stable, which can be seen from Figure 1.Note that in Figure 1, four trajectories start from initial points (0.5, 0.5), (0.1, 0.1), (1.5, 0.1), and (1.5, 0.5), respectively, and all converge (1.2040, 0.2505) as t tends to +∞.Remark 3. Generally speaking, the exploited population should be mature population, which is more in line with the economic and biological views of renewable resource management.If then the positive equilibrium (x * , y * ) of system ( 7) is globally asymptotically stable for τ = 0. From the perspective of eco-managers, in order to plan harvest strategies and maintain the sustainable development of the ecosystem, a positive equilibrium that is asymptotically stable may be required [28].

Global Asymptotic Stability of Positive Equilibrium for τ = 0
In this subsection, we will study the global asymptotical stability of positive equilibrium (x * , y * ) for system (7) by constructing Lyapunov function.
Let x(t) = x * + u(t), y(t) = y * + v(t), linearized system of system ( 7) is where Rewrite the system (28) as follows: In the following, we will construct a positive definite Lyapunov function V(t) and prove that V (t) < 0. Thus, if the global asymptotic stability of origin equilibrium (0, 0) for system (30) is obtained, we can have the global asymptotic stability of positive equilibrium (x * , y * ) for system (7).
First, let V 11 (t) = v 2 (t), then By the inequality 2ab ≤ a 2 + b 2 , we have Second, let Therefore, According to the definition of A 4 and A 5 , we have holds, then V (t) < 0. Through the above analysis, V(t) is a positive definite function and V (t) is negative definite when the conditions ( 10) and (31) hold.Thus, (0, 0) of the system (30) is globally asymptotically stable.That is, the positive equilibrium (x * , y * ) for system (7) is globally asymptotically stable.We can get the following theorem.

Conclusions
In this paper, we investigate dynamics of a stage-structured density-dependent predator-prey system with Beddington-DeAngelis functional response and harvesting.
First, according to the sign for real part of eigenvalue, the local stability of the origin equilibrium (0, 0) and the boundary equilibrium a−E 1 b , 0 is determined.In addition, we discuss local stability of positive equilibrium (x * , y * ) for the case of τ = 0 and τ = 0, respectively.When τ = 0, the positive equilibrium (x * , y * ) of system (7) is locally asymptotically stable.However, when τ = 0, the local stability of positive equilibrium (x * , y * ) of the system (7) includes the following three cases: (i) If (10), (13), and (18) hold, the positive equilibrium (x * , y * ) of system ( 7) is locally asymptotically stable for any τ ≥ 0; (ii) If ( 10) and ( 19) hold, when τ increases and passes τ 0 , the stability of positive equilibrium (x * , y * ) for system (7) switches from local asymptotical stability to unstability and bifurcates into two periodic solutions with small amplitudes; (iii) If ( 10) and ( 21) hold, there is a positive integer k such that the positive equilibrium (x * , y * ) of system (7) has k switches from stability to instability to stability, which shows when the delay τ passes the critical value τ + j , j = 0, 1, 2, ..., k − 1, the positive equilibrium (x * , y * ) of system (7) loses its stability and Hopf bifurcation occurs.
Last, we consider the global asymptotical stability of the positive equilibrium (x * , y * ) for system (7).When τ = 0, stability theory of the periodic solution is used to prove that the positive equilibrium (x * , y * ) is globally asymptotically stable by contradiction; when τ = 0, by constructing Lyapunov functions we prove the positive equilibrium (x * , y * ) for system (7) is globally asymptotically stable.Moreover, examples are given to illustrate the obtained results.
The predator-prey model studied in this paper, where not only the prey density dependence but also the predator density dependence are considered such that the system conforms to real biological environment.In Remarks 1, 2, and 4, we explain the effect of time delay on the stability of boundary equilibrium, the existence of the positive equilibrium, and the global asymptotic stability of the positive equilibrium, respectively.Moreover, we study the effect of harvesting on the global asymptotic stability of the posi-tive equilibrium in Remark 3. Simultaneously, regarding the predator density-dependent mortality rate p, by Theorem 1 and the condition 10, the parameter p does not affect the stability of boundary equilibrium, the existence of the positive equilibrium, respectively.As x * , y * , a i (i = 1, 2, • • •, 4), A j (j = 1, 2, • • •, 5) all contain parameter p, the predator density-dependent mortality rate p not only affects the local asymptotic stability of the positive equilibrium by Theorems 2 and 4, but also the global asymptotic stability of the positive equilibrium by Theorems 5 and 6.However, we cannot describe exactly how the predator density-dependent mortality rate affects the stability of the positive equilibrium.In our future work, we will further study the influence of parameter p on the predator-prey system.