# Dynamics and Optimal Taxation Control in a Bioeconomic Model with Stage Structure and Gestation Delay

Authors: Yi Zhang; Qingling Zhang; Fenglan Bai
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

A prey-predator model with gestation delay, stage structure for predator, and selective harvesting effort on mature predator is proposed, where taxation is considered as a control instrument to protect the population resource in prey-predator biosystem from overexploitation. It shows that interior equilibrium is locally asymptotically stable when the gestation delay is zero, and there is no periodic orbit within the interior of the first quadrant of state space around the interior equilibrium. An optimal harvesting policy can be obtained by virtue of Pontryagin's Maximum Principle without considering gestation delay; on the other hand, the interior equilibrium of model system loses as gestation delay increases through critical certain threshold, a phenomenon of Hopf bifurcation occurs, and a stable limit cycle corresponding to the periodic solution of model system is also observed. Finally, numerical simulations are carried out to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

Recently, the dynamics of a class of stage-structured prey-predator models with gestation delay have been studied by several authors [1–9]. Especially, there is a well-developed theory of stage structured models which incorporates time delay into maturity of population [1]. The prey-predator models with stage structure and gestation delay play an important role in the modelling of multispecies population dynamics. Generally, individuals in each stage are identical in biological characteristics, and the reproduction of mature predator population after predating the prey is not instantaneous but mediates by some discrete time lag required for gestation of mature predator. Xu and Ma [9] proposed the following model,
(1.1)x˙(t)=x(t)(r-αx(t)-a1y2(t)),y˙1(t)=a2x(t-τ)y2(t-τ)-r1y1(t)-dy1(t),y˙2(t)=dy1(t)-r2y2(t)-βy22(t),
where x(t) represents the density of prey population at time t, y1(t) and y2(t) represent the density of immature and mature predator population at time t, respectively; r is the intrinsic growth rate of prey population, r1 is death rate of immature predator and r2 is the death rate of mature predator, α and β represents the intraspecific competition rate of prey and mature predator population, respectively; a2/a1 is the rate of converting prey population into new immature predator population. The constant τ≥0 denotes the gestation delay of mature predator, and τ≥0 is based on the assumption that the reproduction of predator population after predating the prey population is not instantaneous but mediates by some discrete time lag required for gestation of mature predator population. d>0 denotes the proportional transforming rate from immature predator population to mature predator population. All the parameters mentioned above are positive constants.It is well known that exploitation of several biological resources has been increased by the growing human needs for more food and energy, which attracts a global concern to protect the limited biological resources. Consequently, regulation of exploitation of biological resources has become a problem of major concern in view of dwindling resource stocks and deteriorating environment. It should be noted that some techniques and issues associated with bioeconomic exploitation have been discussed in details by Clark [10]. Due to its economic flexibility, taxation is usually considered as possible governing instruments in regulation for harvesting to keep the damage to the ecosystem minimal.Recently, there has been considerable interest in the modeling of harvesting of biological resources. 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. Furthermore, optimal harvesting policies with taxation are also discussed. Ganguly and Chaudhuri [11], Krishna et al. [12], and Dubey et al. [13, 14] investigated the optimal harvesting of a class of models of a single fishery species with taxation as a control. Chaudhuri et al. [15–17], Pradhan and Chaudhuri [18], and Kar et al. [19–25] studied the optimal taxation policies for harvesting of the prey-predator system. However, from the above literature survey, it may be pointed out that no attempt has been made to study the optimal taxation policy of a stage-structured prey-predator system. Furthermore, taxation instrument is discussed to control overharvesting from prey-predator system with gestation delay in [26]. However, stage structure of predator population is not considered, and the periodic orbit within the interior of the first quadrant of state space around interior equilibrium is also not investigated in [26].The stability analysis of interior equilibrium is performed in the third section. It reveals that when gestation delay is zero, the interior equilibrium is locally asymptotically stable. It is also found that equilibrium switch occurs due to variation of gestation delay. Furthermore, an optimal harvesting policy for mature predator is also discussed in the absence of gestation delay. We aim to find an optimal harvesting policy which guarantees an ever-lasting exploitation of the biological resource and maximizes the benefits resulting from the harvesting. Numerical simulations are provided to support the analytical findings in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

Based on the above analysis, the work done by Xu and Ma in [9] is extended by incorporating harvest effort on mature predator, and taxation is chosen to control the conservation of biological resource. In this paper, a prey-predator model with gestation delay and stage structure for predator is established. It is assumed that mature predator is subject to a dynamic harvesting. To conserve the population in the prey-predator ecosystem, the regulatory agency imposes a taxation σ>0 per unit biomass of mature predator (σ<0 denotes the subsidies given to the harvesting effort). Based on the above aspects, the model can be governed by the following differential equations:
(2.1)x˙(t)=x(t)(r-αx(t)-a1y2(t)),y˙1(t)=a2x(t-τ)y2(t-τ)-r1y1(t)-dy1(t),y˙2(t)=dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t),E˙(t)=α0E(t)[(p-σ)qy2(t)-c],
where initial conditions are as follows:
(2.2)x(t)=ψ1(t)>0,y2(t)=ψ2(t)>0,t∈[-τ,0),y1(0)=ψ3(0)>0,E(0)=ψ4(0)>0.The harvesting termE(t) is assumed to be proportional to both stock level and effort, which follows the catch per unit effort hypothesis [10]. The constant q is the catchability coefficient, p is the fixed price per unit of predator species, c is the fixed cost of harvesting per unit of effort, and α0 is called stiffness parameter measuring the strength of reaction of harvesting effort. The parameters mentioned above are all positive constants.


## 3. Qualitative Analysis of Model System

From the view of ecological management, we only concentrate on the interior equilibrium of the model system in this paper, since the biological meaning of the interior equilibrium implies that juvenile preys, mature preys, predators, and harvesting effort on predators all exist, which are relevant to our study.It can be obtained that the only interior equilibrium of the model system (2.1) is P*(x*,y1*,y2*,E*), where x*=(r-a1y2*)/α, y1*=a2y2*(r-a1y2*)/α(r1+d), E*=(a2d(r-a1y2*)-α(r2+βy2*)(r1+d))/αq(r1+d), and y2*=c/(p-σ)q. It is easy to show that interior equilibrium exists, provided the following conditions are satisfied:
(3.1)0<σ<min⁡{p-ca1rq,p-c[a1a2d+αβ(r1+d)]q[a2dr-αr2(r1+d)]},
which provides the range of taxation for the existence of interior equilibrium. This range of taxation may be utilized when the regulatory agency establishes relevant agencies for harvesting.The model system (2.1) can be interpreted as the matrix form:
(3.2)X˙(t)=H(X(t)),
where X(t)=(x(t),y1(t),y2(t),E(t))T∈ℝ4, and H(X(t)) is given as follows,
(3.3)H(X(t))=(H1(X(t))H2(X(t))H3(X(t))H4(X(t)))=(x(t)(r-αx(t)-a1y2(t))a2x(t-τ)y2(t-τ)-r1y1(t)-dy1(t)dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t)α0E(t)[(p-σ)qy2(t)-c]).
Let ℝ+4=[0,∞)4 be the nonnegative octant in ℝ4, then G:ℝ+4+1→ℝ4 is locally Lipschitz and satisfies the condition Hi(X(t))|X∈ℝ4≥0.Due to lemma in [27] and Theorem A.4 in [28], any solution of the model system (2.1) with positive initial conditions exist uniquely, and each component of the solution remains within the interval [0,b) for some b>0. Furthermore, if b<∞, then limsup[x(t)+y1(t)+y2(t)+E(t)]=∞. Hence, this completes the positivity for the solutions of model system (2.1).Now we consider boundedness of positive solutions(x(t),y1(t),y2(t),E(t)), and firstly choose the function W1(t)=x(t-τ)+y1(t).Fort>T1+τ,a1-a2>0 (T1 is some fixed positive time), by calculating the time derivative of W1(t) along the solutions of model system (2.1), we get
(3.4)W˙1(t)=rx(t-τ)-αx2(t-τ)-(a1-a2)x(t-τ)y2(t-τ)-(r1+d)y1(t).By virtue of positiveness of solutionx(t-τ), it is easy to show that
(3.5)rx(t-τ)-αx2(t-τ)≤r24α.Based on the positiveness of solutionx(t-τ),  y2(t-τ) and assumption a1-a2>0, it is easy to show that
(3.6)(a1-a2)x(t-τ)y2(t-τ)>0.Hence, it is easy to show that(3.7)W˙1(t)<r24α-(r1+d)y1(t)<r24α-(r1+d)(x(t-τ)+y1(t)),
which follows that there exists a positive quantity M1 such that 0<W1(t)<M1 for all large t>T1+τ. It proves the boundedness of positive solution x(t), y1(t).LetW2(t)=y2(t), by calculating the time derivative of W2(t) along the solutions of model system (2.1), we have
(3.8)W˙2(t)≤dy1(t)-r2y2(t).By virtue of the positivity of the solutions of model system (2.1) and the boundedness of y1(t) mentioned above, it follows that there exists a positive quantity M2 such that W˙2(t)<M2-r2y2(t) for all large time t>T2 (T2 is some fixed positive time). From the above differential inequality it follows that, there exists a positive quantity M3 such that 0<W2(t)<M3 for all large t>T2, which proves the boundedness of positive solution y2(t).LetW3(t)=y2(t)+E(t), by calculating the time derivative of W3(t) along the solutions of model system (2.1), we have
(3.9)W˙3(t)=dy1(t)-r2y2(t)-βy22(t)-(1-α0(p-σ))qE(t)y2(t)-cα0E(t).By virtue of positiveness and boundedness of solutiony1(t) and y2(t), it follows that there exists a positive quantity M4 such that dy1(t)-r2y2(t)-βy22(t)≤M4. Furthermore, under the following assumption:
(3.10)1-α0(p-σ)>0
it is easy to show that W˙3(t)≤M4-cα0E(t) for for all large t>T3 (T3 is some fixed positive time), which derives that there exists a positive quantity M5 such that 0<W3(t)<M5 for all large t>T3, which proves boundedness of positive solution E(t).Remark 3.1.
Since the components(x(t),y1(t),y2(t)) of solution of model system (2.1) represent the population in the prey-predator system, the positivity implies that the population survives, and the boundedness reveals a natural restriction to growth as a consequence of limited resources. Furthermore, with the purpose of maintaining the sustainable development of prey-predator system, the harvesting cannot increase without any restriction. As analyzed above, the assumption (3.10) provides the range of taxation for the boundedness of harvesting effort. It is an inspiration for people to regulate the harvesting effort by means of economic instrument.
The Jacobian of model system (2.1) evaluated at the only interior equilibrium P* leads to the following characteristic equations:
(3.11)|λ+αx*0a1x*0-a2y2*e-λτλ+(r1+d)-a2x*e-λτ00-dλ+(dy1*y2*+βy2*)qy2*00-α0E*(p-σ)qλ|=0,(3.12)M(λ)+N(λ)e-λτ=0,
where
(3.13)M(λ)=λ4+m1λ3+m2λ2+m3λ+m4,N(λ)=n2λ2+n3λ,m1=r1+d+dy1*y2*+βy2*+αx*,m2=(r1+d)(dy1*y2*+βy2*)+αx*(r1+d)(dy1*y2*+βy2*)+cα0qE*,m3=αx*(r1+d)(dy1*y2*+βy2*)+cα0qE*(r1+d+αx*),m4=cα0qE*αx*(r1+d),n2=-a2dx*,n3=a2dx*(a1y2*-αx*).


### 3.1. Case I: Gestation Delayτ=0

In absence of gestation delay, stability of interior equilibriumP* is investigated, and an optimal harvesting policy with taxation control is also investigated.


#### 3.1.1.  Local Stability Analysis

In the absence of gestation delay (τ=0), model system (2.1) is written as follows:
(3.14)x˙(t)=x(t)(r-αx(t)-a1y2(t)),y˙1(t)=a2x(t)y2(t)-r1y1(t)-dy1(t),y˙2(t)=dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t),E˙(t)=α0E(t)[(p-σ)qy2(t)-c],
and (3.12) can be written as follows:
(3.15)λ4+m1λ3+(m2+n2)λ2+(m3+n3)λ+m4=0.It can be shown that(3.16)m1>0,m4>0,m1(m2+n2)-(m3+n3)=αx*[(r1+d)2+(dy1*y2*+βy2*)2+(r1+d)(dy1*y2*+βy2*)]+(αx*)2(r1+d+dy1*y2*+βy2*)+cα0qE*(dy1*y2*+βy2*)+αβ(r1+d)y2*2(r1+d+dy1*y2*+βy2*)>0,(m1+n1)[(m2+n2)(m3+n3)-(m3+n3)]-m12m4=αx*(r1+d)(dy1*y2*+βy2*)22×[(r1+d+dy1*y2*+βy2*)(r1+d)(dy1*y2*+βy2*+αx*(r1+d+dy1*y2*+βy2*+αx*))22+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+cα0qE*(r1+d+αx*)22×[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+αx*(dy1*y2*+βy2*)2+a2dx*(a1y2*-αx*)[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+αx*(r1+d+dy1*y2*+βy2*+αx*)(r1+d+dy1*y2*+βy2*)]22+(a2dx*)2(r1+d+dy1*y2*+βy2*+αx*)>0.Based on the above analysis, it can be concluded that the roots of (3.15) have negative real parts by using the Routh-Hurwitz criteria [10]. Consequently, the interior equilibrium P* is locally asymptotically stable in absence of gestation delay.Furthermore, letJ* represent the variational matrix of the model system (3.14) at P*, then
(3.17)∫0T‍tr⁡(J(x*,y1*,y2*,E*))dt=∫0T‍tr⁡(-αx*0-a1x*0a2y2*-(r1+d)a2x*00d-(dy1*y2*+βy2*)-qy2*00α0E*(p-σ)q0)dt=-∫0T‍(αx*+r1+d+dy1*y2*+βy2*)dt<0.It can be easily verified that-(αx*+r1+d+(dy1*/y2*)+βy2*)<0 based on the positivity of the solutions of model system.Hence,∫0T‍tr(J(x*,y1*,y2*,E*))dt<0, which eliminates the existence of Hopf bifurcating periodic solution in the vicinity of P*.Subsequently, we will show the nonexistence of periodic orbit encirclingP*. Let h(x(t),y1(t),y2(t),E(t))=1/x(t)y1(t)y2(t)E(t). According to the positivity of solutions of the model system (3.14), it is obvious that h(x(t),y1(t),y2(t),E(t))>0.DefineΔ(x(t),y1(t),y2(t),E(t))=(∂/∂x)(H1h)+(∂/∂y1)(H2h)+(∂/∂y2)(H3h)+(∂/∂E)(H4h), where Hi,  i=1,2,3,4 have been defined before, then we have
(3.18)Δ(x(t),y1(t),y2(t),E(t))=-αy1y2E-a2y12E-dxy22E-βxy1E<0
for x(t), y1(t), y2(t), E(t)>0, since all other parameters are strictly positive.Therefore, there will be no periodic orbit within the interior of the first quadrant of state space aroundP* based on Benedixon-Dulac criterion [29].


#### 3.1.2. Optimal Harvesting Policy

With the purpose of planning harvesting and keeping sustainable development of ecosystem, we design an optimal harvesting policy to maximize the total discounted net revenue from the harvesting using taxation as a control instrument. The path traced out by(x(t),y1(t),y2(t),E(t)) with optimal taxation σ(t) is also investigated.Net economic revenue to the societyπ(x(t),y1(t),y2(t),E(t),σ,t) = Net economic revenue of harvesting + Net economic revenue to the regulatory agency = (p-σ(t))qy2(t)E(t)-cE(t)+σqy2(t)E(t) = (pqy2(t)-c)E(t).Our objective is to maximize the following optimization problem:(3.19)max⁡∫0∞‍e-δt(pqy2(t)-c)dt,
where δ is the instantaneous annual rate of discount, and the optimization problem is subject to the model system (3.14).By using the Pontryagin's Maximum Principle [10], the associated Hamiltonian function is constructed by
(3.20)H(x(t),y1(t),y2(t),E(t),σ(t),t)=e-δt(pqy2(t)-c)E(t)+λ1(t)[x(t)(r-αx(t)-a1y2(t))]+λ2(t)[a2x(t)y2(t)-r1y1(t)-dy1(t)]+λ3(t)[dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t)]+λ4(t)α0E(t)[(p-σ(t))qy2(t)-c],
where λ1, λ2, λ3, λ4 are adjoint variables. σ is the control variable satisfying the constraints σmin≤σ≤σmax. σmax and σmin represent a feasible upper and lower limit of the taxation for harvesting effort, respectively. Specially, σmin<0 implies that subsidies have the effect of increasing the rate of expansion of the harvesting. According to [29], the condition for a singular control to be optimal can be obtained, that is, ∂H/∂σ=0, from which we get
(3.21)λ4(t)=0.For adjoint variablesλi(t), i=1,2,3,4, we have
(3.22)dλ1dt=-∂H∂x,dλ2dt=-∂H∂y1,dλ3dt=-∂H∂y2,dλ4dt=-∂H∂E,(3.23)dλ1dt=(2αx+a1y2-r)λ1-a2y2λ2,(3.24)dλ2dt=(r1+d)λ2-dλ3,(3.25)dλ3dt=a1xλ1-a2xλ2+(2βy2+r2+qE)λ3-pqEe-δt,(3.26)dλ4dt=qy2λ3-e-δt(pqy2-c).Based on (3.21), it follows from (3.26) that
(3.27)λ3(t)=e-δt(p-cqy2).In order to obtain an optimal equilibrium solution, by considering the interior equilibriumP* and solving (3.25),
(3.28)λ1(t)=a2λ2-A1e-δta1,
where A1=(p(2βy2*+r2)-(p-(c/qy2*))(2βy2*+r2+qE*+δ))/x*. By virtue of (3.28), (3.23) can be rewritten as follows:
(3.29)dλ2dt=(2αx*-r)λ2-A2e-δt,
where A2=(2αx*+a1y2*-r+δ)A1/a2.It is easy to obtain the solution of the above linear differential equation(3.30)λ2(t)=A3e-δt,
where A3=(2αx*+a1y2*-r+δ)A1/a2(δx*+2αx*+r), based on (3.30), by solving (3.24) it follows that
(3.31)λ3(t)=r1+d+δdA3e-δt.
Substituting (3.31) into (3.27), we have
(3.32)p-cqy2*=(1+r1+δd)A3,
which provides an equation to the singular path and gives the optimal equilibrium levels of population x*=xδ,y1*=y1δ,y2*=y2δ. Then the optimal equilibrium levels of harvesting effort and taxation can be obtained as follows:
(3.33)Eδ=a2d(r-a1y2δ)-α(r2+βy2δ)(r1+d)αq(r1+d),σδ=p-cqy2δ.Remark 3.2.
According to [30], λi(t)eδt(i=1,2,3,4) represent unusual shadow prices along the singular path. From (3.30), (3.32), and (3.36), it may be concluded that these shadow prices remain constant over time interval in an optimum equilibrium when they strictly satisfy the transversality condition at ∞ [31]. Furthermore, they remain bounded as t→∞.
Considering the interior equilibrium, (3.27) can be written as
(3.34)λ3(t)qy2*=e-δt(pqy2*-c)=e-δt∂π∂E,
which implies that the user's total  cost of harvesting per unit effort is equal to the discounted values of the future price at the steady state effort level.


### 3.2.  Case II: Gestation Delayτ>0

In this section, a stability switch in model system (2.1) due to gestation delay is investigated. Furthermore, a phenomenon of Hopf bifurcation occurs, and a stable limit cycle corresponding to the periodic solution of model system (2.1) is observed.


#### 3.2.1. Local Stability Analysis

Letλ=iω be a root of (3.12), where ω is positive. Substitute λ=iω into (3.12), and separate the real and imaginary parts, then two transcendental equations can be obtained as follows:
(3.35)ω4-m2ω2+m4=-n3ωsin⁡(ωτ)+n2ω2cos⁡(ωτ),(3.36)m1ω3-m3ω=n3ωcos⁡(ωτ)+n2ω2sin⁡(ωτ).By squaring and adding these two equations, it can be obtained that,(3.37)ω8+B1ω6+B2ω4+B3ω2+B4=0,
where B1=m12-2m2, B2=m22+2m4-2m1m3-n22, B3=m32-2m2m4-n32, B4=m42 and mi, nj (i=1,2,3,4;j=2,3) have been defined in (3.12).According to the values ofBi (i=1,2,3,4) and the Routh-Hurwitz criteria [10], a simple assumption of the existence of a positive root for (3.37) is B3<0. If B3<0 holds, then (3.37) has a positive root ω0, and (3.12) has a pair of purely imaginary roots of the form ±iω0. Consequently, it can be obtained by eliminating sin(ωτ) from (3.35) and (3.36) that
(3.38)cos⁡(ωτ)=ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω,
where the τk corresponding to ω0 is as follows,
(3.39)τk=1ω0arccos[ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω]+2kπω0,k=0,1,2,….
By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for τ<τ0, as k=0.


#### 3.2.2. Hopf Bifurcation

In this section, the condition for Hopf bifurcation in [29] is utilized to investigate whether there is a phenomenon of Hopf bifurcation as τ increases through τ0. As stated above, λ=iω0 represents a purely imaginary root of (3.12), and it follows from the above analysis that |M(λ)|=|N(λ)|, which determines a set of possible values of ω0. We will determine the direction of motion of λ=iω0 as τ is varied, namely,
(3.40)Θ=sign⁡[d(Reλ)dτ]λ=iω0=sign⁡[Re(dλdτ)-1]λ=iω0.Theorem 3.3.
Model system (2.1) undergoes Hopf bifurcation at the interior equilibrium P* when tau=τk, k=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P* when τ>τ0 and ∥τ-τ0∥≪1.Proof.
Differentiating (3.12) with respect to τ, we get
(3.41)(dλdτ)-1=4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ,Θ=sign⁡[Re(dλdτ)-1]λ=iω0=sign⁡[Re4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ]λ=iω0=1ω02sign⁡[m1ω06+(m1m2-2m3)ω04+(m2m3-3m1m4)ω02+m3m4(ω04-m2ω02+m4)2+(m3ω0-m1ω03)2].
According to the values of m1,m2,andn1 defined in (3.12), it is easy to show that m1m2-2m3>0 and m2m3-3m1m4>0. Consequently, it follows that sign[d(Reλ)/dτ]λ=iω0>0, which implies there exists at least one eigenvalue with positive real part for τ>τ0, and the condition for Hopf bifurcation in the reference [29] is also satisfied yielding the required periodic solution.
Furthermore, an attracting invariant closed curve bifurcates from interior equilibriumP* when τ>τ0 and ∥τ-τ0∥≪1.


## 3.1. Case I: Gestation Delayτ=0

In absence of gestation delay, stability of interior equilibriumP* is investigated, and an optimal harvesting policy with taxation control is also investigated.


### 3.1.1.  Local Stability Analysis

In the absence of gestation delay (τ=0), model system (2.1) is written as follows:
(3.14)x˙(t)=x(t)(r-αx(t)-a1y2(t)),y˙1(t)=a2x(t)y2(t)-r1y1(t)-dy1(t),y˙2(t)=dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t),E˙(t)=α0E(t)[(p-σ)qy2(t)-c],
and (3.12) can be written as follows:
(3.15)λ4+m1λ3+(m2+n2)λ2+(m3+n3)λ+m4=0.It can be shown that(3.16)m1>0,m4>0,m1(m2+n2)-(m3+n3)=αx*[(r1+d)2+(dy1*y2*+βy2*)2+(r1+d)(dy1*y2*+βy2*)]+(αx*)2(r1+d+dy1*y2*+βy2*)+cα0qE*(dy1*y2*+βy2*)+αβ(r1+d)y2*2(r1+d+dy1*y2*+βy2*)>0,(m1+n1)[(m2+n2)(m3+n3)-(m3+n3)]-m12m4=αx*(r1+d)(dy1*y2*+βy2*)22×[(r1+d+dy1*y2*+βy2*)(r1+d)(dy1*y2*+βy2*+αx*(r1+d+dy1*y2*+βy2*+αx*))22+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+cα0qE*(r1+d+αx*)22×[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+αx*(dy1*y2*+βy2*)2+a2dx*(a1y2*-αx*)[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+αx*(r1+d+dy1*y2*+βy2*+αx*)(r1+d+dy1*y2*+βy2*)]22+(a2dx*)2(r1+d+dy1*y2*+βy2*+αx*)>0.Based on the above analysis, it can be concluded that the roots of (3.15) have negative real parts by using the Routh-Hurwitz criteria [10]. Consequently, the interior equilibrium P* is locally asymptotically stable in absence of gestation delay.Furthermore, letJ* represent the variational matrix of the model system (3.14) at P*, then
(3.17)∫0T‍tr⁡(J(x*,y1*,y2*,E*))dt=∫0T‍tr⁡(-αx*0-a1x*0a2y2*-(r1+d)a2x*00d-(dy1*y2*+βy2*)-qy2*00α0E*(p-σ)q0)dt=-∫0T‍(αx*+r1+d+dy1*y2*+βy2*)dt<0.It can be easily verified that-(αx*+r1+d+(dy1*/y2*)+βy2*)<0 based on the positivity of the solutions of model system.Hence,∫0T‍tr(J(x*,y1*,y2*,E*))dt<0, which eliminates the existence of Hopf bifurcating periodic solution in the vicinity of P*.Subsequently, we will show the nonexistence of periodic orbit encirclingP*. Let h(x(t),y1(t),y2(t),E(t))=1/x(t)y1(t)y2(t)E(t). According to the positivity of solutions of the model system (3.14), it is obvious that h(x(t),y1(t),y2(t),E(t))>0.DefineΔ(x(t),y1(t),y2(t),E(t))=(∂/∂x)(H1h)+(∂/∂y1)(H2h)+(∂/∂y2)(H3h)+(∂/∂E)(H4h), where Hi,  i=1,2,3,4 have been defined before, then we have
(3.18)Δ(x(t),y1(t),y2(t),E(t))=-αy1y2E-a2y12E-dxy22E-βxy1E<0
for x(t), y1(t), y2(t), E(t)>0, since all other parameters are strictly positive.Therefore, there will be no periodic orbit within the interior of the first quadrant of state space aroundP* based on Benedixon-Dulac criterion [29].


### 3.1.2. Optimal Harvesting Policy

With the purpose of planning harvesting and keeping sustainable development of ecosystem, we design an optimal harvesting policy to maximize the total discounted net revenue from the harvesting using taxation as a control instrument. The path traced out by(x(t),y1(t),y2(t),E(t)) with optimal taxation σ(t) is also investigated.Net economic revenue to the societyπ(x(t),y1(t),y2(t),E(t),σ,t) = Net economic revenue of harvesting + Net economic revenue to the regulatory agency = (p-σ(t))qy2(t)E(t)-cE(t)+σqy2(t)E(t) = (pqy2(t)-c)E(t).Our objective is to maximize the following optimization problem:(3.19)max⁡∫0∞‍e-δt(pqy2(t)-c)dt,
where δ is the instantaneous annual rate of discount, and the optimization problem is subject to the model system (3.14).By using the Pontryagin's Maximum Principle [10], the associated Hamiltonian function is constructed by
(3.20)H(x(t),y1(t),y2(t),E(t),σ(t),t)=e-δt(pqy2(t)-c)E(t)+λ1(t)[x(t)(r-αx(t)-a1y2(t))]+λ2(t)[a2x(t)y2(t)-r1y1(t)-dy1(t)]+λ3(t)[dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t)]+λ4(t)α0E(t)[(p-σ(t))qy2(t)-c],
where λ1, λ2, λ3, λ4 are adjoint variables. σ is the control variable satisfying the constraints σmin≤σ≤σmax. σmax and σmin represent a feasible upper and lower limit of the taxation for harvesting effort, respectively. Specially, σmin<0 implies that subsidies have the effect of increasing the rate of expansion of the harvesting. According to [29], the condition for a singular control to be optimal can be obtained, that is, ∂H/∂σ=0, from which we get
(3.21)λ4(t)=0.For adjoint variablesλi(t), i=1,2,3,4, we have
(3.22)dλ1dt=-∂H∂x,dλ2dt=-∂H∂y1,dλ3dt=-∂H∂y2,dλ4dt=-∂H∂E,(3.23)dλ1dt=(2αx+a1y2-r)λ1-a2y2λ2,(3.24)dλ2dt=(r1+d)λ2-dλ3,(3.25)dλ3dt=a1xλ1-a2xλ2+(2βy2+r2+qE)λ3-pqEe-δt,(3.26)dλ4dt=qy2λ3-e-δt(pqy2-c).Based on (3.21), it follows from (3.26) that
(3.27)λ3(t)=e-δt(p-cqy2).In order to obtain an optimal equilibrium solution, by considering the interior equilibriumP* and solving (3.25),
(3.28)λ1(t)=a2λ2-A1e-δta1,
where A1=(p(2βy2*+r2)-(p-(c/qy2*))(2βy2*+r2+qE*+δ))/x*. By virtue of (3.28), (3.23) can be rewritten as follows:
(3.29)dλ2dt=(2αx*-r)λ2-A2e-δt,
where A2=(2αx*+a1y2*-r+δ)A1/a2.It is easy to obtain the solution of the above linear differential equation(3.30)λ2(t)=A3e-δt,
where A3=(2αx*+a1y2*-r+δ)A1/a2(δx*+2αx*+r), based on (3.30), by solving (3.24) it follows that
(3.31)λ3(t)=r1+d+δdA3e-δt.
Substituting (3.31) into (3.27), we have
(3.32)p-cqy2*=(1+r1+δd)A3,
which provides an equation to the singular path and gives the optimal equilibrium levels of population x*=xδ,y1*=y1δ,y2*=y2δ. Then the optimal equilibrium levels of harvesting effort and taxation can be obtained as follows:
(3.33)Eδ=a2d(r-a1y2δ)-α(r2+βy2δ)(r1+d)αq(r1+d),σδ=p-cqy2δ.Remark 3.2.
According to [30], λi(t)eδt(i=1,2,3,4) represent unusual shadow prices along the singular path. From (3.30), (3.32), and (3.36), it may be concluded that these shadow prices remain constant over time interval in an optimum equilibrium when they strictly satisfy the transversality condition at ∞ [31]. Furthermore, they remain bounded as t→∞.
Considering the interior equilibrium, (3.27) can be written as
(3.34)λ3(t)qy2*=e-δt(pqy2*-c)=e-δt∂π∂E,
which implies that the user's total  cost of harvesting per unit effort is equal to the discounted values of the future price at the steady state effort level.


## 3.1.1.  Local Stability Analysis

In the absence of gestation delay (τ=0), model system (2.1) is written as follows:
(3.14)x˙(t)=x(t)(r-αx(t)-a1y2(t)),y˙1(t)=a2x(t)y2(t)-r1y1(t)-dy1(t),y˙2(t)=dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t),E˙(t)=α0E(t)[(p-σ)qy2(t)-c],
and (3.12) can be written as follows:
(3.15)λ4+m1λ3+(m2+n2)λ2+(m3+n3)λ+m4=0.It can be shown that(3.16)m1>0,m4>0,m1(m2+n2)-(m3+n3)=αx*[(r1+d)2+(dy1*y2*+βy2*)2+(r1+d)(dy1*y2*+βy2*)]+(αx*)2(r1+d+dy1*y2*+βy2*)+cα0qE*(dy1*y2*+βy2*)+αβ(r1+d)y2*2(r1+d+dy1*y2*+βy2*)>0,(m1+n1)[(m2+n2)(m3+n3)-(m3+n3)]-m12m4=αx*(r1+d)(dy1*y2*+βy2*)22×[(r1+d+dy1*y2*+βy2*)(r1+d)(dy1*y2*+βy2*+αx*(r1+d+dy1*y2*+βy2*+αx*))22+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+cα0qE*(r1+d+αx*)22×[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+a2dx*(r1+d+dy1*y2*+βy2*+αx*+a1y2*)]+αx*(dy1*y2*+βy2*)2+a2dx*(a1y2*-αx*)[(dy1*y2*+βy2*)((r1+d)(r1+d+dy1*y2*+βy2*)+cα0qE*)2222+αx*(r1+d+dy1*y2*+βy2*+αx*)(r1+d+dy1*y2*+βy2*)]22+(a2dx*)2(r1+d+dy1*y2*+βy2*+αx*)>0.Based on the above analysis, it can be concluded that the roots of (3.15) have negative real parts by using the Routh-Hurwitz criteria [10]. Consequently, the interior equilibrium P* is locally asymptotically stable in absence of gestation delay.Furthermore, letJ* represent the variational matrix of the model system (3.14) at P*, then
(3.17)∫0T‍tr⁡(J(x*,y1*,y2*,E*))dt=∫0T‍tr⁡(-αx*0-a1x*0a2y2*-(r1+d)a2x*00d-(dy1*y2*+βy2*)-qy2*00α0E*(p-σ)q0)dt=-∫0T‍(αx*+r1+d+dy1*y2*+βy2*)dt<0.It can be easily verified that-(αx*+r1+d+(dy1*/y2*)+βy2*)<0 based on the positivity of the solutions of model system.Hence,∫0T‍tr(J(x*,y1*,y2*,E*))dt<0, which eliminates the existence of Hopf bifurcating periodic solution in the vicinity of P*.Subsequently, we will show the nonexistence of periodic orbit encirclingP*. Let h(x(t),y1(t),y2(t),E(t))=1/x(t)y1(t)y2(t)E(t). According to the positivity of solutions of the model system (3.14), it is obvious that h(x(t),y1(t),y2(t),E(t))>0.DefineΔ(x(t),y1(t),y2(t),E(t))=(∂/∂x)(H1h)+(∂/∂y1)(H2h)+(∂/∂y2)(H3h)+(∂/∂E)(H4h), where Hi,  i=1,2,3,4 have been defined before, then we have
(3.18)Δ(x(t),y1(t),y2(t),E(t))=-αy1y2E-a2y12E-dxy22E-βxy1E<0
for x(t), y1(t), y2(t), E(t)>0, since all other parameters are strictly positive.Therefore, there will be no periodic orbit within the interior of the first quadrant of state space aroundP* based on Benedixon-Dulac criterion [29].


## 3.1.2. Optimal Harvesting Policy

With the purpose of planning harvesting and keeping sustainable development of ecosystem, we design an optimal harvesting policy to maximize the total discounted net revenue from the harvesting using taxation as a control instrument. The path traced out by(x(t),y1(t),y2(t),E(t)) with optimal taxation σ(t) is also investigated.Net economic revenue to the societyπ(x(t),y1(t),y2(t),E(t),σ,t) = Net economic revenue of harvesting + Net economic revenue to the regulatory agency = (p-σ(t))qy2(t)E(t)-cE(t)+σqy2(t)E(t) = (pqy2(t)-c)E(t).Our objective is to maximize the following optimization problem:(3.19)max⁡∫0∞‍e-δt(pqy2(t)-c)dt,
where δ is the instantaneous annual rate of discount, and the optimization problem is subject to the model system (3.14).By using the Pontryagin's Maximum Principle [10], the associated Hamiltonian function is constructed by
(3.20)H(x(t),y1(t),y2(t),E(t),σ(t),t)=e-δt(pqy2(t)-c)E(t)+λ1(t)[x(t)(r-αx(t)-a1y2(t))]+λ2(t)[a2x(t)y2(t)-r1y1(t)-dy1(t)]+λ3(t)[dy1(t)-r2y2(t)-βy22(t)-qE(t)y2(t)]+λ4(t)α0E(t)[(p-σ(t))qy2(t)-c],
where λ1, λ2, λ3, λ4 are adjoint variables. σ is the control variable satisfying the constraints σmin≤σ≤σmax. σmax and σmin represent a feasible upper and lower limit of the taxation for harvesting effort, respectively. Specially, σmin<0 implies that subsidies have the effect of increasing the rate of expansion of the harvesting. According to [29], the condition for a singular control to be optimal can be obtained, that is, ∂H/∂σ=0, from which we get
(3.21)λ4(t)=0.For adjoint variablesλi(t), i=1,2,3,4, we have
(3.22)dλ1dt=-∂H∂x,dλ2dt=-∂H∂y1,dλ3dt=-∂H∂y2,dλ4dt=-∂H∂E,(3.23)dλ1dt=(2αx+a1y2-r)λ1-a2y2λ2,(3.24)dλ2dt=(r1+d)λ2-dλ3,(3.25)dλ3dt=a1xλ1-a2xλ2+(2βy2+r2+qE)λ3-pqEe-δt,(3.26)dλ4dt=qy2λ3-e-δt(pqy2-c).Based on (3.21), it follows from (3.26) that
(3.27)λ3(t)=e-δt(p-cqy2).In order to obtain an optimal equilibrium solution, by considering the interior equilibriumP* and solving (3.25),
(3.28)λ1(t)=a2λ2-A1e-δta1,
where A1=(p(2βy2*+r2)-(p-(c/qy2*))(2βy2*+r2+qE*+δ))/x*. By virtue of (3.28), (3.23) can be rewritten as follows:
(3.29)dλ2dt=(2αx*-r)λ2-A2e-δt,
where A2=(2αx*+a1y2*-r+δ)A1/a2.It is easy to obtain the solution of the above linear differential equation(3.30)λ2(t)=A3e-δt,
where A3=(2αx*+a1y2*-r+δ)A1/a2(δx*+2αx*+r), based on (3.30), by solving (3.24) it follows that
(3.31)λ3(t)=r1+d+δdA3e-δt.
Substituting (3.31) into (3.27), we have
(3.32)p-cqy2*=(1+r1+δd)A3,
which provides an equation to the singular path and gives the optimal equilibrium levels of population x*=xδ,y1*=y1δ,y2*=y2δ. Then the optimal equilibrium levels of harvesting effort and taxation can be obtained as follows:
(3.33)Eδ=a2d(r-a1y2δ)-α(r2+βy2δ)(r1+d)αq(r1+d),σδ=p-cqy2δ.Remark 3.2.
According to [30], λi(t)eδt(i=1,2,3,4) represent unusual shadow prices along the singular path. From (3.30), (3.32), and (3.36), it may be concluded that these shadow prices remain constant over time interval in an optimum equilibrium when they strictly satisfy the transversality condition at ∞ [31]. Furthermore, they remain bounded as t→∞.
Considering the interior equilibrium, (3.27) can be written as
(3.34)λ3(t)qy2*=e-δt(pqy2*-c)=e-δt∂π∂E,
which implies that the user's total  cost of harvesting per unit effort is equal to the discounted values of the future price at the steady state effort level.


## 3.2.  Case II: Gestation Delayτ>0

In this section, a stability switch in model system (2.1) due to gestation delay is investigated. Furthermore, a phenomenon of Hopf bifurcation occurs, and a stable limit cycle corresponding to the periodic solution of model system (2.1) is observed.


### 3.2.1. Local Stability Analysis

Letλ=iω be a root of (3.12), where ω is positive. Substitute λ=iω into (3.12), and separate the real and imaginary parts, then two transcendental equations can be obtained as follows:
(3.35)ω4-m2ω2+m4=-n3ωsin⁡(ωτ)+n2ω2cos⁡(ωτ),(3.36)m1ω3-m3ω=n3ωcos⁡(ωτ)+n2ω2sin⁡(ωτ).By squaring and adding these two equations, it can be obtained that,(3.37)ω8+B1ω6+B2ω4+B3ω2+B4=0,
where B1=m12-2m2, B2=m22+2m4-2m1m3-n22, B3=m32-2m2m4-n32, B4=m42 and mi, nj (i=1,2,3,4;j=2,3) have been defined in (3.12).According to the values ofBi (i=1,2,3,4) and the Routh-Hurwitz criteria [10], a simple assumption of the existence of a positive root for (3.37) is B3<0. If B3<0 holds, then (3.37) has a positive root ω0, and (3.12) has a pair of purely imaginary roots of the form ±iω0. Consequently, it can be obtained by eliminating sin(ωτ) from (3.35) and (3.36) that
(3.38)cos⁡(ωτ)=ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω,
where the τk corresponding to ω0 is as follows,
(3.39)τk=1ω0arccos[ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω]+2kπω0,k=0,1,2,….
By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for τ<τ0, as k=0.


### 3.2.2. Hopf Bifurcation

In this section, the condition for Hopf bifurcation in [29] is utilized to investigate whether there is a phenomenon of Hopf bifurcation as τ increases through τ0. As stated above, λ=iω0 represents a purely imaginary root of (3.12), and it follows from the above analysis that |M(λ)|=|N(λ)|, which determines a set of possible values of ω0. We will determine the direction of motion of λ=iω0 as τ is varied, namely,
(3.40)Θ=sign⁡[d(Reλ)dτ]λ=iω0=sign⁡[Re(dλdτ)-1]λ=iω0.Theorem 3.3.
Model system (2.1) undergoes Hopf bifurcation at the interior equilibrium P* when tau=τk, k=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P* when τ>τ0 and ∥τ-τ0∥≪1.Proof.
Differentiating (3.12) with respect to τ, we get
(3.41)(dλdτ)-1=4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ,Θ=sign⁡[Re(dλdτ)-1]λ=iω0=sign⁡[Re4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ]λ=iω0=1ω02sign⁡[m1ω06+(m1m2-2m3)ω04+(m2m3-3m1m4)ω02+m3m4(ω04-m2ω02+m4)2+(m3ω0-m1ω03)2].
According to the values of m1,m2,andn1 defined in (3.12), it is easy to show that m1m2-2m3>0 and m2m3-3m1m4>0. Consequently, it follows that sign[d(Reλ)/dτ]λ=iω0>0, which implies there exists at least one eigenvalue with positive real part for τ>τ0, and the condition for Hopf bifurcation in the reference [29] is also satisfied yielding the required periodic solution.
Furthermore, an attracting invariant closed curve bifurcates from interior equilibriumP* when τ>τ0 and ∥τ-τ0∥≪1.


## 3.2.1. Local Stability Analysis

Letλ=iω be a root of (3.12), where ω is positive. Substitute λ=iω into (3.12), and separate the real and imaginary parts, then two transcendental equations can be obtained as follows:
(3.35)ω4-m2ω2+m4=-n3ωsin⁡(ωτ)+n2ω2cos⁡(ωτ),(3.36)m1ω3-m3ω=n3ωcos⁡(ωτ)+n2ω2sin⁡(ωτ).By squaring and adding these two equations, it can be obtained that,(3.37)ω8+B1ω6+B2ω4+B3ω2+B4=0,
where B1=m12-2m2, B2=m22+2m4-2m1m3-n22, B3=m32-2m2m4-n32, B4=m42 and mi, nj (i=1,2,3,4;j=2,3) have been defined in (3.12).According to the values ofBi (i=1,2,3,4) and the Routh-Hurwitz criteria [10], a simple assumption of the existence of a positive root for (3.37) is B3<0. If B3<0 holds, then (3.37) has a positive root ω0, and (3.12) has a pair of purely imaginary roots of the form ±iω0. Consequently, it can be obtained by eliminating sin(ωτ) from (3.35) and (3.36) that
(3.38)cos⁡(ωτ)=ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω,
where the τk corresponding to ω0 is as follows,
(3.39)τk=1ω0arccos[ω4+m1ω3-m2ω2-m3ω+m4n2ω2+n3ω]+2kπω0,k=0,1,2,….
By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for τ<τ0, as k=0.


## 3.2.2. Hopf Bifurcation

In this section, the condition for Hopf bifurcation in [29] is utilized to investigate whether there is a phenomenon of Hopf bifurcation as τ increases through τ0. As stated above, λ=iω0 represents a purely imaginary root of (3.12), and it follows from the above analysis that |M(λ)|=|N(λ)|, which determines a set of possible values of ω0. We will determine the direction of motion of λ=iω0 as τ is varied, namely,
(3.40)Θ=sign⁡[d(Reλ)dτ]λ=iω0=sign⁡[Re(dλdτ)-1]λ=iω0.Theorem 3.3.
Model system (2.1) undergoes Hopf bifurcation at the interior equilibrium P* when tau=τk, k=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P* when τ>τ0 and ∥τ-τ0∥≪1.Proof.
Differentiating (3.12) with respect to τ, we get
(3.41)(dλdτ)-1=4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ,Θ=sign⁡[Re(dλdτ)-1]λ=iω0=sign⁡[Re4λ3+3m1λ2+2m2λ+m3λ4+m1λ3+m2λ2+m3λ+m4+2n2λ+n3n2λ2+n3λ-τλ]λ=iω0=1ω02sign⁡[m1ω06+(m1m2-2m3)ω04+(m2m3-3m1m4)ω02+m3m4(ω04-m2ω02+m4)2+(m3ω0-m1ω03)2].
According to the values of m1,m2,andn1 defined in (3.12), it is easy to show that m1m2-2m3>0 and m2m3-3m1m4>0. Consequently, it follows that sign[d(Reλ)/dτ]λ=iω0>0, which implies there exists at least one eigenvalue with positive real part for τ>τ0, and the condition for Hopf bifurcation in the reference [29] is also satisfied yielding the required periodic solution.
Furthermore, an attracting invariant closed curve bifurcates from interior equilibriumP* when τ>τ0 and ∥τ-τ0∥≪1.


## 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.


### 4.1. Numerical Simulation for Optimal Harvesting Policy

In this subsection, values of parameters are taken from [9] which are used in Example 1 of [9] and set in appropriate units, r=5, α=4, β=1, a1=3, a2=2, r1=0.1, r2=0.1, d=2, p=13, q=0.18, α0=0.08, c=1. For the model system (3.14), the range of the taxation σ∈(0.5,7.0848) can be obtained based on (3.1) and (3.10). According to [10], we take the instantaneous annual rate of discount δ=0.05 in appropriate units. According to the given values of parameters, (3.32) can be numerically computed and three roots can be obtained as follows: y2*=0.5212,1.6323,2.2935. Based on (p-σ)q=y2*, the corresponding σi(i=1,2,3) can be calculated, σ1=2.3405,σ2=9.5965, and σ3=10.5777, respectively. It is obvious that only σ1=2.3405 satisfies the range (0.5,7.0848). Consequently, the optimal taxation is σδ=σ1=2.3405, then the optimal equilibrium levels of population and harvest effort can be also obtained (xδ,y1δ,y2δ,Eδ)=(0.8591,0.4245,0.5212,5.6399), which are indicated in Figure 1.Figure 1
Dynamical responses of model system (3.14) with the optimal taxation.


### 4.2. Numerical Simulation for the Hopf Bifurcation

In this subsection, values of parameters are taken from [9], which are used in Example 3 of [9] and set in appropriate units, r=2, α=0.5, β=0.5, a1=3, a2=2, r1=1, r2=0.1, d=1, p=13, q=0.18, α0=0.08, andc=1. It follows from (3.1) and (3.10) that the taxation range is (0.5, 3.7407), and it can be obtained that population densities in model system (2.1) is (x*,y1*,y2*)=(0.4,0.24,0.6) with σ=2.3. It should be noted that σ=2.3 is arbitrarily selected from the interval (0.5,3.7407), which can guarantee the existence of interior equilibrium of model system (2.1). Furthermore, it can be also calculated that B3<0, which satisfies the assumption of the existence of a positive root for (3.37), and then τ0=0.8734 is calculated based on (3.39). By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for τ<τ0, which can be seen in Figure 2. It should be noted that τ=0.3 is randomly selected in the interval (0,0.8734), which is enough to merit the above mathematical study.Figure 2
Dynamical responses of model system (2.1) with discrete time delay τ=0.3.According to Theorem3.3 in this paper, a periodic solution caused by the phenomenon of Hopf bifurcation and a limit cycle corresponding to this periodic solution occurs as τ increases through τ0, which are shown in Figures 3 and 4, respectively.Figure 3
Dynamical responses of model system (2.1) with discrete time delay τ=0.8734.Figure 4
A limit cycle for the model system (2.1) corresponding to the periodic solution shown in Figure 3.


## 4.1. Numerical Simulation for Optimal Harvesting Policy

In this subsection, values of parameters are taken from [9] which are used in Example 1 of [9] and set in appropriate units, r=5, α=4, β=1, a1=3, a2=2, r1=0.1, r2=0.1, d=2, p=13, q=0.18, α0=0.08, c=1. For the model system (3.14), the range of the taxation σ∈(0.5,7.0848) can be obtained based on (3.1) and (3.10). According to [10], we take the instantaneous annual rate of discount δ=0.05 in appropriate units. According to the given values of parameters, (3.32) can be numerically computed and three roots can be obtained as follows: y2*=0.5212,1.6323,2.2935. Based on (p-σ)q=y2*, the corresponding σi(i=1,2,3) can be calculated, σ1=2.3405,σ2=9.5965, and σ3=10.5777, respectively. It is obvious that only σ1=2.3405 satisfies the range (0.5,7.0848). Consequently, the optimal taxation is σδ=σ1=2.3405, then the optimal equilibrium levels of population and harvest effort can be also obtained (xδ,y1δ,y2δ,Eδ)=(0.8591,0.4245,0.5212,5.6399), which are indicated in Figure 1.Figure 1
Dynamical responses of model system (3.14) with the optimal taxation.


## 4.2. Numerical Simulation for the Hopf Bifurcation

In this subsection, values of parameters are taken from [9], which are used in Example 3 of [9] and set in appropriate units, r=2, α=0.5, β=0.5, a1=3, a2=2, r1=1, r2=0.1, d=1, p=13, q=0.18, α0=0.08, andc=1. It follows from (3.1) and (3.10) that the taxation range is (0.5, 3.7407), and it can be obtained that population densities in model system (2.1) is (x*,y1*,y2*)=(0.4,0.24,0.6) with σ=2.3. It should be noted that σ=2.3 is arbitrarily selected from the interval (0.5,3.7407), which can guarantee the existence of interior equilibrium of model system (2.1). Furthermore, it can be also calculated that B3<0, which satisfies the assumption of the existence of a positive root for (3.37), and then τ0=0.8734 is calculated based on (3.39). By virtue of Butler's lemma [32], it can be concluded that the interior equilibrium remains stable for τ<τ0, which can be seen in Figure 2. It should be noted that τ=0.3 is randomly selected in the interval (0,0.8734), which is enough to merit the above mathematical study.Figure 2
Dynamical responses of model system (2.1) with discrete time delay τ=0.3.According to Theorem3.3 in this paper, a periodic solution caused by the phenomenon of Hopf bifurcation and a limit cycle corresponding to this periodic solution occurs as τ increases through τ0, which are shown in Figures 3 and 4, respectively.Figure 3
Dynamical responses of model system (2.1) with discrete time delay τ=0.8734.Figure 4
A limit cycle for the model system (2.1) corresponding to the periodic solution shown in Figure 3.


## 5. Conclusion

In this paper, a bioeconomic model is proposed to investigate dynamics of the effects of a stage-structured prey-predator system with harvesting effort and gestation delay. Theoretical analysis shows that the interior equilibrium is locally asymptotically stable around interior equilibrium when the model system is in absence of discrete time delay. By using Pantryagin's Maximum Principle, an optimal harvesting policy with taxation is derived to ensure the sustainable development of biological resource and prosperous commercial harvesting. It reveals that the user's total  cost of harvest per unit effort must be equal to the discounted value of the future price at the steady state level. In the case of gestation delay, the stability analysis reveals that gestation delay is responsible for the stability switch of model system. A phenomenon of Hopf bifurcation occurs as the discrete time delay increases through a certain threshold.It should be noted that taxation instrument is discussed to control overharvesting from prey-predator system with gestation delay in [26]. Compared with work done in [26], stage structure of predator population is considered, and the periodic orbit within the interior of the first quadrant of state space around interior equilibrium is also investigated in this paper. The work done in [9] is extended by incorporating the harvesting effort into the prey-predator system, and taxation is adopted as a controlling instrument to regulate harvesting of predator. From the qualitative analysis of the model, the effect of harvesting effort is extensively investigated with and without discrete time delay. Compared with model system investigated in this paper, the harvesting effort is not considered in [9], the interior equilibrium becomes unstable as τ=0.7. However, the interior equilibrium of the model system (2.1) remains stable as τ=0.7 based on the analysis in this paper. It implies that the harvesting effort has an effect of stabilizing the interior equilibrium, and the cyclic behavior can be prevented by applying the harvesting effort into the model system.



---

*Source: 803270-2012-07-10.xml*