# Global Stability Analysis and Optimal Control of a Harvested Ecoepidemiological Prey Predator Model with Vaccination and Taxation

Authors: Chao Liu; Qingling Zhang; Jinna Li
Journal: Abstract and Applied Analysis
 (2013)

---

## Abstract

We propose an ecoepidemiological prey predator model, where selective harvest effort on predator population is considered. Vaccination and taxation are introduced as control instruments, which are utilized to control number of susceptible prey population and protect predator population from overexploitation, respectively. Conditions which influence nonnegativity and boundedness of solutions are studied. Global stability analysis around disease-free equilibrium is discussed based on robust Bendixson criterion, which is theoretically beneficial to studying coexistence and interaction mechanism of population within harvested ecoepidemiological system. By using Pontryagin’s maximum principle, an optimal control strategy is derived to maximize the total discounted net economic revenue to society as well as protect prey population from infectious disease. Numerical simulations are carried out to show the consistency with theoretical analysis.


---

## Body

## 1. Introduction

In the natural world, prey predator ecosystem with infectious disease is a well-studied ecological phenomenon. Generally, species within prey predator ecosystem does not exist alone. While species spreads the disease, it also competes with the other species for space or food, or is predated by other species. Kermack-Mc kendrick made the pioneering work on (susceptible-infective-removal-susceptible) SIRS [1]. Along with the line of this research, mathematical epidemiology has become an important subject of research in recent decades, where the evolution of a disease which gets transmitted upon contact is described. In the work done by Hadeler and Freedman [2], a prey predator model where the prey is infected by a parasite and in turn infects the predator with the parasite. Thereafter, plenty of papers were made to investigate prey predator system with infection in the prey. Arino et al. [3, 4] observed that the introduction of an infected population in the classical ratio-dependent predator-prey model may act as a biological control to save the population from extinction. Hethcote et al. [5, 6] considered the effects of disease-induced mortality or disease-reduced reproduction in regulating natural populations, decreasing their population sizes, reducing their natural fluctuations, or causing destabilization of equilibria into oscillations of the population states. In the work done by Haque and Venturino [7], the mixed terms of the underlying demographic model describing the mass action law of populations affect each other, and Holling-type functions are taken into account [8]. In the above mentioned work, by performing dynamic analysis with the help of linear approximation to model systems that consist of nonlinear equations, it concludes that it is of importance from ecological and mathematical points of view to study ecological systems subject to epidemiological factors.Recently, an eco-epidemiological model consisting of susceptible prey, infected prey, and predator is proposed and analyzed in [9]. It is assumed that the predation functional response (trophic function) allows Holling type II (Michaelis-Menten kinetics). From the above assumptions model system can be expressed as follows:
(1)ds(t)dt=rs(t)(1-s(t)+i(t)k)-a1s(t)p(t)1+b1s(t)+b2i(t)-cs(t)i(t),di(t)dt=cs(t)i(t)-a2i(t)p(t)1+b1s(t)+b2i(t)-ei(t),dp(t)dt=(d1a1s(t)+d2a2i(t))p(t)1+b1s(t)+b2i(t)-fp(t),
where s(t) is the susceptible prey population densities per hectare and i(t) is the infected prey population densities per hectare, so that the total prey population densities at any time t per hectare is s(t)+i(t), and it is assumed that the infected prey population contributes with the susceptible prey population to the prey population growth towards the environment carrying capacity. Consequently, the susceptible prey population grows according to a logistic law with carrying capacity k>0 and intrinsic growth rate constant r>0 involving the whole prey population (susceptible prey and infected prey population), which is best regarded as a purely descriptive equation. The transmission rate among the susceptible prey population and the infected prey population follows the simple law of mass action cs(t)i(t) for all t>0, where c>0 denotes the transmission coefficient. The infectious disease spreads among the prey population only and that disease is not genetically inherited. The infected prey population does not recover or become immune. p(t) is the predator population densities per hectare which is assumed to feed on both the susceptible and the infected prey population. e denotes the natural death rate of the infected prey populations; f is the death rate of the predator population; d1 and d2 which are within the interval (0,1) are the conversion rate of the susceptible and the infected prey population to the predator population, respectively. a1 and a2 are the searching efficiency constants or equivalently the predation rates on the susceptible and the infected prey populations, respectively. b1 and b2 are the positive constants.It is well known that harvesting has a strong impact on the dynamic evolution of a population. Evidence shows that many species have already become extinct and many others are at the verge of extinction due to several natural or man-made reasons like over exploitation, indiscriminate harvesting and mismanagement of natural resources, and so forth. The severity of this impact that may range from rapid depletion to complete preservation of a population depends on the implemented harvesting agency [10, 11]. Due to practical and economic utilization, biological resources in the prey predator system are extensively harvested nowadays. Furthermore, exploitation of biological resources has been increased by people’s multifarious material needs, which attracts a global concern to protect the limited biological resources. Therefore, regulation of exploitation of biological resources has become a problem of major concern in view of dwindling resource stocks and the deteriorating environment. It is necessary to establish a constructive management of commercial exploitation of the biological resources. The techniques and issues associated with bioeconomic exploitation have been discussed in details by Clark [12]. Taxation, license fees, lease of property rights, seasonal harvesting, and so forth are usually considered as possible governing instruments in regulation for harvesting as well as regulatory mechanisms to keep the damage to the ecosystem minimal. Out of these regulating options, taxation is considered to be superior because of its economic flexibility [13], and economists are particularly attracted to taxation because a competitive system can be better maintained under taxation rather than other regulatory methods. There has been considerable interest in the modeling of harvesting of biological resources in recent decades, where harvest effort is considered to be a dynamic variable, and optimal harvesting policies with taxation are discussed. Ganguly and Chaudhuri [14], Krishna et al. [15], and Dubey et al. [16, 17] investigated the optimal harvesting of a class of models of a single fishery species with taxation as a control. Kar [18, 19] studied the optimal tax policies for harvesting of the prey predator system.Recently, stimulated by new analytical and numerical findings, scholars have found interest to discuss the applications of optimal control theory to eco-epidemiological model, which can be found in [20–22]. Mathematical modelling of infectious disease with application of optimal control theory is investigated in [23–25], where vaccination cost is considered to discuss associated optimal problems. It should be noted that harvesting issues are also considered in eco-epidemiological model [26, 27], where dynamical behaviour and impulsive harvesting control are investigated. However, to authors’ best knowledge, optimal control strategy associated with maximum net economic revenue of harvesting and vaccination cost in eco-epidemiological model has not been investigated in related research.In this paper, a harvested eco-epidemiological prey predator system is established in the second section. We extend the work done in [9] by incorporating the harvest effort on predator population. Based on superiority of taxation control in regulating harvesting and effectiveness of vaccination in protecting population from infectious disease introduced above, vaccination is introduced to control number of susceptible prey population, and taxation is chosen to regulate exploitation of biological resource within the harvested eco-epidemiological prey predator system. In the third section, nonnegativity and boundedness of solutions of the proposed model are discussed. In the fourth section, global stability analysis of model system around disease-free equilibrium is studied based on robust Bendixson criterion. Furthermore, by taking taxation and vaccination control as time dependent, an optimal strategy associated with net harvesting revenue and vaccination cost is investigated by using Pontryagin’s maximum principle in the fifth section. It is assumed that vaccination cost is considered and exploitation of the predator species is regulated by an authority by imposing taxation per unit biomass of the harvested predator, and the imposition of taxation acts as a deterrent to control harvesting. We aim to find an optimal strategy which guarantees a lasting exploitation of the biological resource and maximize the total discounted net economic revenue to society as well as protect prey population from infectious disease. Numerical simulations are provided to support the analytical findings obtained in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

In 1954, Gordon [28] proposed the economic theory of a common property resource, which studies the effect of harvest effort on ecosystem from an economic perspective. It is well known that the harvested predator is usually sold as commodity in the market in order to achieve economic interest, which implies that the unit price of harvested population is influenced by the fluctuation of supply and demand for harvested population in the market, and the unit cost of harvest effort depends on fecundity of harvested population [12, 28]. Consequently, the unit price and unit cost are not always constant, but time varying. Assuming that there is always constant demand for harvested predator, then v1(t) the unit price of harvested predator and v2(t) the unit cost of harvest effort can be expressed as
(2)v1(t)=h-σβ+E(t)p(t),v2(t)=gp(t),
respectively, where E(t) represents the harvest effort on predator population p(t) at any time t>0, h, σ, β and g are all positive constants. To conserve the population in the harvested eco-epidemiological prey predator system, the regulatory agency imposes a tax σ>0 per unit biomass of the predator, and σ<0 denotes the subsidies given to the harvest effort. E(t)p(t) denotes supply amount of harvested predator into the market.By simple computation, it can be found that the unit price will decrease if the supply of harvested predator is larger than the demand for harvested predator; on the other hand, the unit price will increase if the supply of harvested predator cannot meet the demand for harvested predator, and the maximum unit price is(h-σ)/β. It is also obvious that
(3)limp(t)→0v2(t)=limp(t)→0gp(t)=∞,limp(t)→∞v2(t)=limp(t)→∞gp(t)=0,
which imply that the unit cost of harvest effort is inversely proportional to the population density of harvested predator. More richer the population density of harvested predator becomes, more lower the unit cost of harvest effort is; on the other hand, the unit cost of harvest effort will increase if the population density of harvested predator becomes relatively rare.According to the above analysis, a harvested eco-epidemiological prey predator model with vaccination and taxation is proposed based on (1) and (2), which can be governed by the following differential equations:
(4)ds(t)dt=rs(t)(1-s(t)+i(t)k)-a1s(t)p(t)1+b1s(t)+b2i(t)-cs(t)i(t)-us(t),di(t)dt=cs(t)i(t)-a2i(t)p(t)1+b1s(t)+b2i(t)-ei(t),dp(t)dt=(d1a1s(t)+d2a2i(t))p(t)1+b1s(t)+b2i(t)-fp(t)-E(t)p(t),dE(t)dt=αE(t)p(t)(h-σβ+E(t)p(t)-gp(t)).
The initial conditions of model system (4) are given as follows:
(5)s(0)≥0,i(0)≥0,p(0)≥0,E(0)≥0,
where α describes stiffness parameter measuring the strength of reaction of harvest effort and susceptible prey population is vaccinated at vaccination rate u. Other parameters that share the same interpretations are mentioned in (1) and (2).Remark 1.
It should be noted that it is assumed that harvesting issue in model system (4) does not affect the growth of prey population (susceptible prey and infected prey) directly. Effect of harvest effort and harvesting regulation on prey population is biologically functioned through the predator population.Remark 2.
It is further assumed that vaccination rateu such that constraint 0≤u<1 for any time, which vividly reflects the fact that vaccination on the entire susceptible prey population is impossible.


## 3. Nonnegativity and Boundedness of Solution

Theorem 3.
All solutions of model system (4) with initial conditions are nonnegative.Proof.
The model system (4) can be interpreted as the matrix form:
(6)dX(t)dt=H(X(t)),
where X(t)=(s(t),i(t),p(t),E(t))T∈ℝ4 and H(X(t)) is given as follows:(7)H(X(t))=(H1(X(t))H2(X(t))H3(X(t))H4(X(t)))=(rs(t)(1-s(t)+i(t)k)-a1s(t)p(t)1+b1s(t)+b2i(t)-cs(t)i(t)-us(t)cs(t)i(t)-a2i(t)p(t)1+b1s(t)+b2i(t)-ei(t)(d1a1s(t)+d2a2i(t))p(t)1+b1s(t)+b2i(t)-fp(t)-E(t)p(t)αE(t)p(t)(h-σβ+E(t)p(t)-gp(t))).Let ℝ+4=[0,∞)4 be the nonnegative octant in ℝ4, then L:ℝ+4+1→ℝ4 is locally Lipschitz and satisfy the condition Hi(X(t))|X∈ℝ4≥0,i=1,2,3,4.
Due to lemma in [29] and [30, Theorem A.4], any solution of the model system (4) with nonnegative initial conditions exist uniquely, and each component of the solution remains within the interval [0,A0) for some A0>0. Standard and simple arguments show that solutions of model system (4) always exist and stay nonnegative. Hence, this completes the nonnegativity of the solutions of model system (4).Theorem 4.
Ifh-σ>0 and d1≥d2, then all solutions of model system (4) are bounded within a region Ω for all t>0,
(8)Ω={kd1r4m1(s(t),i(t),p(t),E(t))∣s(t)≤k,E(t)≤h-σg,w1(t)≤k(r+e)24er,w2(t)≤kd1r4m1}.
On the other hand, ifh-σ>0 and d1<d2, then all solutions of model system (4) are bounded within a region Ω for all t>0,
(9)Ω={ked1r2+ck2m1(d2-d1)(r+e)24erm1(s(t),i(t),p(t),E(t))∣s(t)≤k,E(t)≤h-σg,w1(t)≤k(r+e)24er,w2(t)≤ked1r2+ck2m1(d2-d1)(r+e)24erm1},
where w1(t)=s(t)+i(t), w2(t)=d1s(t)+d2i(t)+p(t), and m1=min{e,u,f}.Proof.
Firstly, we consider the first equation of model system (4). We have
(10)ds(t)dt≤rs(t)(1-s(t)k),
which implies that
(11)limt→∞sups(t)≤k.
Now, defining the functionw1(t)=s(t)+i(t) by calculating the time derivative of w1(t) along the solutions of model system (4), we get
(12)dw1(t)dt≤(r+e)s(t)-rs2(t)k-ew1(t).
By utilizing the standard comparison principle, it follows that
(13)limt→∞supw1(t)≤k(r+e)24er.
Furthermore, by using nonnegativity of solution of model system (4), if h-σ>0, then it derives that
(14)dE(t)dt≤α(h-σ)-αgE(t),
which implies that
(15)limt→∞supE(t)≤h-σg.
Letw2(t)=d1s(t)+d2i(t)+p(t), and by calculating the time derivative of w2(t) along the solutions of model system (4), it can be obtained that
(16)dw2(t)dt≤rd1s(t)(1-s(t)k)+c(d2-d1)s(t)i(t)-m1(d1s(t)+d2i(t)+p(t))=rd1s(t)(1-s(t)k)+c(d2-d1)s(t)i(t)-m1w2(t).
Hence, it derives that(17)limt→∞supw2(t)={kd1r4m1,d1≥d2,ked1r2+cm1k2(d2-d1)(r+e)24erm1,d2>d1,
where m1=min{e,u,f}.
Consequently, Ifh-σ>0 and d1≥d2, then all solutions of model system (4) are bounded within a region Ω for all t>0,
(18)Ω={kd1r4m1(s(t),i(t),p(t),E(t))∣s(t)≤k,E(t)≤h-σg,w1(t)≤k(r+e)24er,w2(t)≤kd1r4m1}.
On the other hand, ifh-σ>0 and d1<d2, then all solutions of model system (4) are bounded within a region Ω for all t>0,
(19)Ω={ked1r2+ck2m1(d2-d1)(r+e)24erm1(s(t),i(t),p(t),E(t))∣s(t)≤k,E(t)≤h-σg,w1(t)≤k(r+e)24er,w2(t)≤ked1r2+ck2m1(d2-d1)(r+e)24erm1}.Remark 5.
Since the components(s(t),  i(t),  and  p(t)) of solution of model system (4) represent the population in the harvested eco-epidemiological prey predator system, boundedness reveals a natural restriction to growth as a consequence of limited resources. Especially, boundedness of i(t) epidemiologically interprets infectious disease in prey population, it can not spread without limitation, and prey population may not die out due to rampant spread of infectious disease. Furthermore, with the purpose of maintaining the sustainable development of prey predator system, the harvesting can not increase without any restriction, which is reflected by the boundedness of E(t). It is of inspiration for people to regulate the harvesting effort by means of economic instrument.


## 4. Global Stability Analysis around Disease-Free Equilibrium

In this section, global stability analysis around disease-free equilibrium is discussed. It is utilized to discuss stability mechanism of model system due to variation of taxation and vaccination control. It should be noted that we only concentrate on the dynamical behavior and stability analysis around disease-free equilibrium, since biological interpretation of disease-free equilibrium implies that susceptible prey, predator population all exist under infectious disease-free environment, only susceptible prey population is predated and harvest effort on predator population exists, which are relevant to our main study in this paper.For simplicity, model system (4) is nondimensionalized with the following scaling: S=s/k,  I=i/k,  P=a1p/r,  Y=E/r,  T=rt. By using these quantities, model system (4) can be transformed into a dimensionless form as follows:
(20)dS(T)dT=S(T)(1-S(T)-I(T))-S(T)P(T)1+b3S(T)+b4I(T)-c1S(T)I(T)-u1S(T),dI(T)dT=c1S(T)I(T)-a3I(T)P(T)1+b3S(T)+b4I(T)-e1I(T),dP(T)dT=(d3S(T)+d4I(T))P(T)1+b3S(T)+b4I(T)-f1P(T)-P(T)Y(T),dY(T)dT=α1P(T)Y(T)(h1-σ1β1+P(T)Y(T)-g1P(T)),
where c1=ck/r, e1=e/r, a3=a2/a1, b3=b1k, b4=b2k, f1=f/r, d3=(d1a1k)/r, d4=d2a2k/r, α1=α/a1, β1=(a1β)/r2, h1-σ1=(a1(h-σ))/r2, u1=u/r, and g1=(a1g)/r.By simple computation, disease-free equilibriumM*(S*,0,P*,Y*) can be obtained as follows.P
*
=
(
1
-
u
1
-
S
*
)
(
1
+
b
3
S
*
), Y*=((d3-b3f1)S*-f1)/(1+b3S*), and S* satisfies the following equation:
(21)[g1(d3-b3f1)-b3(h1-σ1)]S*2+[(1-u1)b3-1]×(h1-σ1+f1g1)S*-(1-u1)d3g1S*+(h1-σ1+g1f1)(1-u1)-g1β1=0.Conditions for existence of unique positive real root of (21) are as follows:
(22)g1(β1-f1)≤σ1≤h1+g1f1-g1d3b3,1-b3(β1+f1)+d3≤u1<r.Remark 6.
The above inequalities associated with taxation and vaccination control in (22) also guarantee existence and uniqueness of disease-free equilibrium M*.Theorem 7.
If the following inequalities hold(23)(a3+e1)>(a3+c1)S*,(S*-1)[b3(h1-σ1)+2g1(b3f1-d3)]>0,u1>b3-1-2d3b3,
then model system (20) is locally stable around disease-free equilibrium M*.Proof.
For disease-free equilibriumM*(S*,0,P*,Y*), the variational matrix of model system (20) around M* takes the following form:
(24)V(S*,0,P*,Y*)=(V11V12V13V14V21V22V23V24V31V32V33V34V41V42V43V44),
where
(25)V11=S*[(1-u1)b3-1-2b3S*]1+b3S*,V12=S*[b4-c1-1-(b4+b3+b3c1)S*]1+b3S*,V13=-S*1+b3S*,V22=(c1+a3)S*-a3-e1,V31=d3(1-S*)1+b3S*,V32=(1-S*)[d4+(b3d4-b4d3)S*]1+b3S*,V33=-d3S*1+b3S*,V34=(S*-1)(1+b3S*),V43=Δ*[(d3-b3f1)S*-f1]1+b3S*,V44=Δ*(1-S*)(1+b3S*)-α1g1,V14=V21=V23=V24=V41=V42=0.
It follows fromV(S*,0,P*,Y*) that the characteristic equation around M* is as follows:
(26)(λ+(a3+e1)-(a3+c1)S*)×(λ3+A1*λ2+A2*λ+A3*)=0,where
(27)A1*=α1[(S*-1)(1+b3S*)Δ*+g1],Δ*=β1(h1-σ1)[β1+(1-S*)(S*(d3-b3f1)-f1)]2,A2*=(2b3(b3+d3)S*3+(3b3-b32-2d3)S*2+(2d3-b3+u1b3+1)S*S*2)×((1+b3S*)2)-1+α1Δ*(S*-1)(1+b3S*)2[(d3-b3f1)S*-f1](1+b3S*)2,A3*=(S*A1*[2b3(b3+d3)S*2+(3b3-b32-2d3)S*b3+fd3+2d3-b3+u1b3+12b3(b3+d3)S*2])×((1+b3S*)2)-1+((b3+d3)S*2α1S*A1*(A1*-g1)(2b3S*-b3+1)×[(d3-b3f1)S*-f12b3(b3+d3)S*2])×((1+b3S*)2)-1.
It follows from Routh-Hurwitz criteria [13] that there are four eigenvalues with negative real part provided that A1*>0,A3*>0, and A1*A2*>A3*, which derives that
(28)(a3+e1)>(a3+c1)S*,(S*-1)[b3(h1-σ1)+2g1(b3f1-d3)]>0,u1>b3-1-2d3b3.In the following part, global stability analysis of model system (20) around disease-free equilibrium M* is investigated. Firstly, some preliminaries are discussed as follows.Lemma 8 (see [31]).
The unique interior equilibrium is globally stable in the domainG if it is locally stable and all trajectories in G converge to the unique interior equilibrium, where domain G satisfies the following conditions:(i)
the domainG is simply connected;(ii)
there is a compact absorbing setJ⊂G.LetV|2| be the second additive compound matrix [31, 32]. In the case of n=4, V|2| can be written based on variational matrix of model system (20) around M*, which is as follows:(29)V|2|=(V11+V22V23V24-V13-V140V32V11+V33V34V120-V14V42V43V11+V440V12V13-V31V210V22+V33V34-V24-V410V21V43V22+V44V230-V41V31-V42V32V33+V44),where Vij,  i,j=1,2,3,4, can be found in Theorem 7. Since V14=0, V21=0, V23=0, V24=0, V41=0, and V42=0, V|2| can be deduced as follows:(30)V|2|=(V11+V2200-V1300V32V11+V33V34V12000V43V11+V440V12V13-V3100V22+V33V340000V43V22+V44000V310V32V33+V44).Furthermore, we have(31)N-1=diag{YS,YS,YS,YS},NF(Z)=dNdZ=diag{S˙Y-SY˙Y2,S˙Y-SY˙Y2,S˙Y-SY˙Y2},
where F(Z(T)) and Z(T) are defined based on model system (20), and the model system (20) can be interpreted as the matrix form(32)dZ(T)dT=F(Z(T)),
where Z(T)=(S(T),I(T),P(T),Y(T))T and F(Z(T)) is given as follows:(33)F(Z(T))=(S(T)(1-S(T)-I(T))-S(T)P(T)1+b3S(T)+b4I(T)-c1S(T)I(T)-u1S(T)c1S(T)I(T)-a3I(T)P(T)1+b3S(T)+b4I(T)-e1I(T)(d3S(T)+d4I(T))P(T)1+b3S(T)+b4I(T)-f1P(T)-P(T)Y(T)α1P(T)Y(T)(h1-σ1β1+P(T)Y(T)-g1P(T))).By simple computation, it can be obtained that(34)NFN-1=diag{S˙S-Y˙Y,S˙S-Y˙Y,S˙S-Y˙Y},NV|2|N-1=V|2|.Let(35)D=NFN-1+NV|2|N-1=(D11D12D13D21D22D23D31D32D33),where(36)D11=S˙S-Y˙Y+V11+V22,D12=(0,0),D13=(-V13,0,0),D21=(V32,0)T,D22=(S˙S-Y˙Y+V11+V33V34V43S˙S-Y˙Y+V11+V44),D23=(V12000V12V13),D31=(-V31,0,0)T,D32=(00000V31),D33=(S˙S-Y˙Y+V22+V33V340V43S˙S-Y˙Y+V22+V4400V32S˙S-Y˙Y+V33+V44).Consider the Lozinskii measureΓ of D with respect to a vector norm |·| in ℝ4 (see, [32]) and I (identity matrix), then it can be obtained that
(37)Γ(D)=limϵ→0+|I+ϵD|-1ϵ.Hence, it follows from simple computation that(38)Γ(D)≤sup{f1,f2,f3},
where fi=Γ1(Dii)+|Dij| for i,j=1,2,3, and i≠j. |D12|, |D13|, |D21|, |D23|, |D31|, and |D32| are matrix norms with respect to the L1 vector norm, and Γ1 is the Lozinskii measure with respect to this vector norm.Furthermore, letΔ=β1(h1-σ1)/[β1+(1-S)(S(d3-b3f1)-f1)]2, then it can be computed that
(39)Γ1(D11)=S˙S-Y˙Y+|(1-u1)b3S-S-2b3S21+b3S+(c1+a3)S-(a3+e1)(1-u1)b3S-S-2b3S21+b3S|,|D12|=0,|D13|=S1+b3S,Γ1(D22)=S˙S-Y˙Y+max{d3S1+b3S,α1Δ|(1-S)(1+b3S)|d3S1+2b3S+α1g1d3S1+b3S}+|S[(1-u1)b3-1-2b3S]1+b3S|+max{α1Δ|(d3-b3f1)S-f1|1+b3S(1-S)(1+b3S),ffffffffα1Δ|(d3-b3f1)S-f1|1+b3S},|D21|=d3S1+b3S,|D23|=max{S|b4-1-c1-(b4+b3+b3c1)S|1+b3S,S1+b3S},|D31|=|D32|=d3(1-S)1+b3S,Γ1(D33)=S˙S-Y˙Y+|(α1Δ+1)(1-S)(1+b3S)|+d3S1+b3S+|(c1+a3)S-(a3+e1)|+α1max{Δ|S(d3-b3f1)-f1|,g1}.According to Theorem3 and definition of fi,i=1,2,3, it follows that
(40)f1≤S˙S-Y˙Y+2+c1+e1,f2≤S˙S-Y˙Y+2+max{1,(b4-1-c1)24[b4+b3(1+c1)]}+max{(1+b3)24b3,(16α1β1(h1-σ1)(d3-b3f1)2mmmmmmmmmmm×[d3+f1(b3+1)](d3-b3f1)2)ffggfff×([((1+b3)f1-d3)2ffffgggfff+4β1(d3-b3f1)]2)-1(1+b3)24b3}+max{+4β1(d3-b3f1)((1+b3)f1-d3)2]2)-11,(16α1β1(h1-σ1)(1+b3)(d3-b3f1)2)ggggggg×(+4β1(d3-b3f1)((1+b3)f1-d3)2]2[((1+b3)f1-d3)2gggggggggg+4β1(d3-b3f1)((1+b3)f1-d3)2]2)-1+α1g1},f3≤S˙S-Y˙Y+1+c1+e1+(1+16α1β1(h1-σ1)(d3-b3f1)2[((1+b3)f1-d3)2+4β1(d3-b3f1)]2)×[1+(b3+1)24b3]+max{+4β1(d3-b3f1)((1+b3)f1-d3)2]2)-1α1g1,((d3-b3f1)216α1β1(h1-σ1)mmmmmmml×f1(d3+b3+1)(d3-b3f1)2)gggggg×(+4β1(d3-b3f1)((1+b3)f1-d3)2]2[((1+b3)f1-d3)2mmmmmmm+4β1(d3-b3f1)((1+b3)f1-d3)2]2)-1}.Lemma 9.
For elementP(T) of the disease-free equilibrium M*, some inequalities about liminfT→∞P(T) hold as follows:  (a)
ifa3d3+c1(a3d3-d4)≤0, then
(41)liminfT→∞P(T)≥d3(m2+1-u1)24m2-d32+d4m3d3;(b)
ifa3d3+c1(a3d3-d4)>0, then
(i)
liminfT→∞P(T)≥[a3d3(1-u1+m2-m3)-c1m3(a3d3-d4)]2/4d3m2a32  −  (d32+d4m3)/d3, provided that a3d3(1-u1+m2-m3)≥c1m3(a3d3-d4),(ii)
liminfT→∞P(T)≥(a3d3(m2-m3)-c1m3(a3d3-d4))/a3m2-(d32+d4m3)/d3, provided that a3d3(1-u1+m2-m3)<c1m3(a3d3-d4), where m2=max{e1,f1+r(h1-σ1)/g1} and m3=k(e1+1)2/4e1.Proof.
Based on Theorem4 and nondimensionalized scaling of model system (20) mentioned in Section 4, it derives that
(42)Y(T)≤r(h1-σ1)g1,I(T)≤k(1+e1)24e1,S(T)≤1.
According to model system (20), let W(T)=d3S(T)+(d4/a3)I(T)+P(T). By calculating the time derivative of W(T) along the solutions of model system (20), it can be obtained that
(43)dWdT=d3S(1-S-I)-c1SI(d3-d4a3)-d4e1Ia3-f1P-PY-d3u1S=d3S(1-u1-S)-[d3+c1(d3-d4a3)]SI-d4e1Ia3-f1P-PY≥d3(1-u1+m2)S-d3S2-[d3+c1(d3-d4a3)]SI-m2W,
where m2=max{e1,f1+r(h1-σ1)/g1}.
Ifa3d3+c1(a3d3-d4)≤0, then
(44)dWdT≥d3(1-u1+m2)S-d3S2-m2W.
According toW=d3S+(d4/a3)I+P, it follows from the standard comparison principle that
(45)liminfT→∞P(T)≥d3(m2+1-u1)24m2-d32+d4m3d3:=n1.
On the other hand, ifa3d3+c1(a3d3-d4)>0, then
(46)dWdT≥[d3(1-u1+m2-m3)-c1m3(d3-d4a3)]S-d3S2-m2W,
where m3=k(e1+1)2/4e1.
According toW=d3S+(d4/a3)I+P, it follows from the standard comparison principle that
(47)liminfT→∞P(T)≥[a3d3(1-u1+m2-m3)-c1m3(a3d3-d4)]24d3m2a32-d32+d4m3d3:=n2,
provided that a3d3(1-u1+m2-m3)≥c1m3(a3d3-d4).
Furthermore,(48)liminfT→∞P(T)≥a3d3(m2-m3)-c1m3(a3d3-d4)a3m2-d32+d4m3d3:=n3,
provided that a3d3(1-u1+m2-m3)<c1m3(a3d3-d4).Let(49)l1=α1g1+1+c1+e1,l2=1+max{1,(b4-1-c1)24[b4+b3(1+c1)]}+max{(1+b3)24b3,(16α1β1(h1-σ1)(d3-b3f1)2gggggbbbbbgbbbggg×[d3+f1(b3+1)](d3-b3f1)2)fffggfff×([((1+b3)f1-d3)2mmmmmmmm+4β1(d3-b3f1)]2)-1(1+b3)24b3}+max{+4β1(d3-b3f1)]2((1+b3)f1-d3)2)-11,(16α1β1(h1-σ1)(1+b3)(d3-b3f1)2)dddfddd×([((1+b3)f1-d3)2cccccccccccc+4β1(d3-b3f1)((1+b3)f1-d3)2]2((1+b3)f1-d3)2)-1+α1g1},l3=c1+e1+(1+16α1β1(h1-σ1)(d3-b3f1)2[((1+b3)f1-d3)2+4β1(d3-b3f1)]2)×[1+(b3+1)24b3]+max{([((1+b3)f1-d3)2+4β1(d3-b3f1)]2)-1+4β1(d3-b3f1)((1+b3)f1-d3)2]2((1+b3)f1-d3)2)-1α1g1,((d3-b3f1)216α1β1(h1-σ1)ggggggffgggg+f1(d3+b3+1)(d3-b3f1)2)dddffddd×(+4β1(d3-b3f1)((1+b3)f1-d3)2]2((1+b3)f1-d3)2[((1+b3)f1-d3)2gggffffvvvvv+4β1(d3-b3f1)((1+b3)f1-d3)2]2((1+b3)f1-d3)2)-1}.Theorem 10.
If the model system is locally stable aroundM* and inequality (22) holds, then(a)
ifa3d3+c1(a3d3-d4)≤0 and α1g1n1(h1-σ1)/(β1g1+rn1(h1-σ1))>l1+l2+l3, then model system (20) is globally stable around M*;(b)
ifa3d3+c1(a3d3-d4)>0, then
(i)
model system (20) is globally stable around M* provided that  a3d3(1-u1+m2-m3)≥c1m3(a3d3-d4),  α1g1n2(h1-σ1)/(β1g1+rn2(h1-σ1))>l1  +  l2+l3,(ii)
model system (20) is globally stable around M* provided that  a3d3(1-u1+m2-m3)<c1m3(a3d3-d4),  α1g1n3(h1-σ1)/(β1g1+rn3(h1-σ1))>l1  +  l2+l3, where ni,i=1,2,3, and mi,i=2,3, have been defined in Lemma 9. li,i=1,2,3, have been defined in (49)–(9).Proof.
According to model system (20), if inequality (22) holds, then there exists a unique disease-free equilibrium M*.
Ifa3d3+c1(a3d3-d4)≤0 and α1g1n1(h1-σ1)/(β1g1+rn1(h1-σ1))>l1+l2+l3, then
(50)Γ(D)≤sup{f1,f2,f3}≤S˙S-α1g1n1(h1-σ1)β1g1+rn1(h1-σ1)+l1+l2+l3,
which derives that
(51)1T∫0T‍Γ(D)dS≤1TlogS(T)S(0)-α1g1n1(h1-σ1)β1g1+rn1(h1-σ1)+l1+l2+l3,limsupT→+∞1T∫0T‍Γ(D)dS<-[α1g1n1(h1-σ1)β1g1+rn1(h1-σ1)-l1-l2-l3]<0.
Based on robust Bendixson criterion [31, 32], inequality (51) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ4, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around M*, then it follows from Lemma 8 that model system (20) is globally stable around M*.
Ifa3d3+c1(a3d3-d4)>0, a3d3(1-u1+m2-m3)≥c1m3(a3d3-d4), and α1g1n2(h1-σ1)/(β1g1+rn2(h1-σ1))>l1+l2+l3, then
(52)Γ(D)≤sup{f1,f2,f3}≤S˙S-α1g1n2(h1-σ1)β1g1+rn2(h1-σ1)+l1+l2+l3,
which derives that
(53)1T∫0T‍Γ(D)dS≤1TlogS(T)S(0)-α1g1n2(h1-σ1)β1g1+rn2(h1-σ1)+l1+l2+l3,limsupT→+∞1T∫0T‍Γ(D)dS<-[α1g1n2(h1-σ1)β1g1+rn2(h1-σ1)-l1-l2-l3]<0.
Based on robust Bendixson criterion [31, 32], inequality (53) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ4, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around M*, then it follows from Lemma 8 that model system (20) is globally stable around M*.
Ifa3d3+c1(a3d3-d4)>0, a3d3(1-u1+m2-m3)<c1m3(a3d3-d4), and α1g1n3(h1-σ1)/(β1g1+rn3(h1-σ1))>l1+l2+l3, then
(54)Γ(D)≤sup{f1,f2,f3}≤S˙S-α1g1n3(h1-σ1)β1g1+rn3(h1-σ1)+l1+l2+l3.
It derives that
(55)1T∫0T‍Γ(D)dS≤1TlogS(T)S(0)-α1g1n3(h1-σ1)β1g1+rn3(h1-σ1)+l1+l2+l3,limsupT→+∞1T∫0T‍Γ(D)dS<-[α1g1n3(h1-σ1)β1g1+rn3(h1-σ1)-l1-l2-l3]<0.
Based on robust Bendixson criterion [31, 32], inequality (55) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ4, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around M*, then it follows from Lemma 8 that model system (20) is globally stable around M*.


## 5. Optimal Control of Model System

With the purpose of regulating harvesting and protecting population from infectious disease within harvested eco-epidemiological system, we design an optimal control to maximize total discounted net revenue to society by using taxation and vaccination as control instruments. The path traced out by(S,I,P,Y) with optimal taxation σ1 and optimal vaccination u1 is also investigated.Total discounted net economic revenue to societyπ(S,I,P,Y,σ1,u1,T)  = discounted economic revenue of harvesting + discounted economic revenue to regulatory agency − vaccination cost.It follows from model system (20), that the following can be obtained:
(56)π(S,I,P,Y,σ1,u1,T)=e-δT[(h1-σ1)PYβ1+PY-g1Y]+e-δTσ1PYβ1+PY-12Bu12=e-δT(h1PYβ1+PY-g1Y)-12Bu12,
where δ is the instantaneous annual rate of discount and the optimization problem is subject to the model system (20). B is a positive weight parameter associated with the square of control variable to balance size of the terms; see [20, 21, 24, 25] and references therein. (1/2)Bu12 represents vaccination cost of susceptible prey population.Remark 11.
Vaccination cost discussed in this section is a quadratic cost, which is the simplest and widest used nonlinear representation of vaccination cost; see [20, 21, 24, 25] and references therein. Actually, when nonlinear description is adopted, it may be not completely clear which nonlinear form must be appropriately chosen. In this case, choosing the simplest form compatible with mathematical requirements (existence, well-posedness of maximization problems, etc.) is a possible guideline; see [20, 21, 24] and references therein. Based on the above introduction, we consider a quadratic cost on the optimal control problem in this section. It should be noted that other nonlinear functions might provide a better description of the actual vaccination cost due to the increase of vaccination cost when most of population is already removed (vaccinated or immune); see [21, 24] and references therein.


### 5.1. Design of Optimal Control

Our objective is to maximize the following optimization problem:(57)max∫0∞‍[e-δT(h1PYβ1+PY-g1Y)-12Bu12]dT.By using the Pontryagin’s maximum principle [12], the associated Hamiltonian function is constructed by
(58)R(S,I,P,Y,σ1,u1,T)=e-δT(h1PYβ1+PY-g1Y)-12Bu12+λ1(T)×[S(1-S-I)-SP1+b3S+b4I-c1SI-u1S]+λ2(T)(c1SI-a3IP1+b3S+b4I-e1I)+λ3(T)[(d3S+d4I)P1+b3S+b4I-f1P-PY]+λ4(T)α1PY(h1-σ1β1+PY-g1P),
where λi(T),i=1,2,3,4, are adjoint variables. σ1 and u1 are the control variables. According to [33], the conditions for a singular control to be optimal can be obtained by  ∂R/∂σ1=0 and ∂R/∂u1=0, from which we get
(59)λ4(T)α1PYβ1+PY=0,λ1(T)S+Bu1=0.According to Theorem3, it follows that
(60)λ4(T)=0,u1=-λ1(T)SB.For adjoint variablesλi(T),i=1,2,3,4, we have
(61)dλ1(T)dT=-∂R∂S,dλ2(T)dT=-∂R∂I,dλ3(T)dT=-∂R∂P,dλ4(T)dT=-∂R∂Y.By substitutingλ4(T)=0 into the Equation dλ4(T)/dT=-∂R/∂Y, it follows that
(62)λ3(T)=e-δT(h1β1(β1+PY)2-g1P).By virtue of (62), it gives that
(63)dλ1(T)dT=[2S+(c1+1)I+P(1+b4I)(1+b3S+b4I)2-1+u1]λ1-I[c1+a3b3P(1+b3S+b4I)2]λ2+e-δT[g1-h1β1P(β1+PY)2](b4d3-b3d4)I+d3(1+b3S+b4I)2.By associating withM*(S*,0,P*,Y*), the above equation can be deduced as follows:
(64)dλ1(T)dT=Q2λ1-d3(1-S*)Q11+b3S*e-δT,
where Q1=h1β1/(β1+P*Y*)2-g1/P* and Q2=[2b3S*-b3+1+b3u1]S*/(1+b3S*).By solving (64), it can be obtained that
(65)λ1(T)=d3(1-S*)Q1(δ+Q2)(1+b3S*)e-δT.By virtue of (62), it gives that
(66)dλ2(T)dT=S[(1+c1)-b4P(1+b3S+b4I)2]λ1+[a3P(1+b3S)(1+b3S+b4I)2+e1-c1S]λ2+e-δT(g1-h1β1P(β1+PY)2)d4+(b3d4-b4d3)S(1+b3S+b4I)2.By associating withM*(S*,0,P*,Y*), the above equation can be deduced as follows:
(67)dλ2(T)dT=Q3λ1+Q5λ2+(S*+u1-1)[d4+(b3d4-b4d3)S*]Q11+b3S*e-δT,
where Q3=[(b4+b3+b3c1)S*-b4+c1+1+b4u1]S*/(1+b3S*) and Q5=[a3+e1-(c1+a3)S*].According to (65), (67) can be rewritten as follows:
(68)dλ2(T)dT=Q5λ2-(S*-1)Q1Q41+b3S*e-δT,
where Q4=d3Q3/(δ+Q2)(1+b3S*)-[d4+(b3d4-b4d3)S*].By solving (68), it can be obtained that
(69)λ2(T)=(S*+u1-1)Q1Q4(1+b3S*)(δ+Q5)e-δT.Based on (65) and (69), it can be obtained that
(70)dλ3(T)dT=h1β1Ye-δT(β1+PY)2+Sλ11+b3S+b4I+a3Iλ21+b3S+b4I+(f1+Y)λ3.By associating withM*(S*,0,P*,Y*), (70) can be deduced as follows:
(71)dλ3(T)dT=h1β1Y*e-δT(β1+P*Y*)2+S*λ11+b3S*+d3S*1+b3S*λ3=d3S*1+b3S*λ3+[h1β1Y*(β1+P*Y*)2+d3S*(1-S*)Q1(δ+Q2)(1+b3S*)2]e-δT.By solving the above equation, it can be obtained that(72)λ3(T)=h1β1d3S*(S*-1)-h1β1Y*(δ+Q2)(1+b3S*)2(δ+Q2)(1+b3S*)(β1+P*Y*)2[δ+(b3δ+d3)S*]×e-δT.In order to obtain an optimal equilibrium solution, it can be obtained by considering the disease-free equilibriumM* and substituting (72) into (62), and it derives that
(73)1-g1(β+P*Y*)2P*h1β1=d3S*(S*-1)-Y*(δ+Q2)(1+b3S*)2(δ+Q2)(1+b3S*)[δ+(b3δ+d3)S*].Equation (73) provides an equation to the singular path and gives the optimal equilibrium levels of population S*=Sδ,  P*=(1-u1-Sδ)(1+b3Sδ), Y*=((d3-b3f1)Sδ-f1)/(1+b3Sδ). Then, the optimal equilibrium levels of taxation and vaccination control can be obtained as follows:
(74)u1δ=(4b3d3Sδ(Sδ-1)e-δTBb3-1-(b3δ+2b3)Sδ+((4b3d3Sδ(Sδ-1)e-δTB)1/2[b3-1-(b3δ+2b3)Sδ]2+(4b3d3Sδ(Sδ-1)e-δTB)1/2)×(2b3)-1,σδ=([g1(d3-b3f1)-h1b3]Sδ2+h1(1-u1)+g1[f1(1-u1)-β1][g1(d3-b3f1)-h1b3]Sδ2)×((1+b3Sδ)(1-u1-Sδ))-1+([(b3-b3u1-1)h1-g1[f1+(1-u1)mmmmmmmmmmmmmn×(d3-b3f1)]]Sδ)×((1+b3Sδ)(1-u1-Sδ))-1.Remark 12 (according to [13]).
λ
i
(
T
)
e
δ
T
(
i
=
1,2
,
3,4
) represent unusual shadow prices along the singular path. From (65), (69), and (72), 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 ∞ [34]. Furthermore, they remain bounded as T→∞.Generally, there is a strong ecological relationship between predator population and susceptible prey population in the harvested eco-epidemiological system, and these two species are usually influenced by each other in several aspects. The imposition of taxation and vaccination defined in (74) act as a deterrent to regulate harvest effort on predator population and control the spread of infectious disease within eco-epidemiological prey predator system, which can not only guarantee the sustainable development of predator population, but also benefit the ecological balance in the ecosystem and sustainable development of other relevant species in the harvested eco-epidemiological system.


### 5.2. Numerical Simulation

With the help of MATLAB, simulation work with hypothetical set of parameters is performed to understand the theoretical results which have been established in Section5.1 of this paper. The values of parameters are cited from [9, 12], which are given as follows: r=1,  k=1,  a1=1,  b1=2.34,  b2=3,  c=0.01,  a2=0.344,  e=0.051,  d1=0.21,  d2=0.6774,  f=0.001,  α=0.4,  h=1.67,  β=1.1,  g=0.41,  δ=0.05,  and   B=1.According to the parameters given above, it follows from (22) that the range of the taxation σ1∈[0.4496,1.6336] and vaccination u1∈[-1.3663,1), which guarantees that there exists a unique disease-free equilibrium. Based on (74), it can be computed that the optimal taxation and optimal vaccination control is σδ=1.0653 and u1δ=0.483, respectively. Since σδ∈[0.4496,1.6336] and u1δ∈[-1.3663,1), it can be obtained that there exists a unique disease-free equilibrium M* corresponding to the optimal taxation and vaccination control, that is, M*(0.5648,0,0.9406,0.0501). According to Theorem 10, model system (20) is globally stable around M*, whose dynamical responses can be shown in Figure 1, and the corresponding phase portraits are shown in Figure 2 (S-I-P space), Figure 3 (S-Y-P space), and Figure 4 (Y-P-I space) with different initial values.Figure 1
Dynamical responses of model system (20) with the optimal taxation and vaccination σδ=1.0653 and u1δ=0.483.Figure 2
Phase portrait (S-I-P space) of model system (20) with σδ=1.0653, and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).Figure 3
Phase portrait (S-Y-P space) of model system (20) with σδ=1.0653 and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).Figure 4
Phase portrait (Y-P-I space) of model system (20) with σδ=1.0653 and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).


## 5.1. Design of Optimal Control

Our objective is to maximize the following optimization problem:(57)max∫0∞‍[e-δT(h1PYβ1+PY-g1Y)-12Bu12]dT.By using the Pontryagin’s maximum principle [12], the associated Hamiltonian function is constructed by
(58)R(S,I,P,Y,σ1,u1,T)=e-δT(h1PYβ1+PY-g1Y)-12Bu12+λ1(T)×[S(1-S-I)-SP1+b3S+b4I-c1SI-u1S]+λ2(T)(c1SI-a3IP1+b3S+b4I-e1I)+λ3(T)[(d3S+d4I)P1+b3S+b4I-f1P-PY]+λ4(T)α1PY(h1-σ1β1+PY-g1P),
where λi(T),i=1,2,3,4, are adjoint variables. σ1 and u1 are the control variables. According to [33], the conditions for a singular control to be optimal can be obtained by  ∂R/∂σ1=0 and ∂R/∂u1=0, from which we get
(59)λ4(T)α1PYβ1+PY=0,λ1(T)S+Bu1=0.According to Theorem3, it follows that
(60)λ4(T)=0,u1=-λ1(T)SB.For adjoint variablesλi(T),i=1,2,3,4, we have
(61)dλ1(T)dT=-∂R∂S,dλ2(T)dT=-∂R∂I,dλ3(T)dT=-∂R∂P,dλ4(T)dT=-∂R∂Y.By substitutingλ4(T)=0 into the Equation dλ4(T)/dT=-∂R/∂Y, it follows that
(62)λ3(T)=e-δT(h1β1(β1+PY)2-g1P).By virtue of (62), it gives that
(63)dλ1(T)dT=[2S+(c1+1)I+P(1+b4I)(1+b3S+b4I)2-1+u1]λ1-I[c1+a3b3P(1+b3S+b4I)2]λ2+e-δT[g1-h1β1P(β1+PY)2](b4d3-b3d4)I+d3(1+b3S+b4I)2.By associating withM*(S*,0,P*,Y*), the above equation can be deduced as follows:
(64)dλ1(T)dT=Q2λ1-d3(1-S*)Q11+b3S*e-δT,
where Q1=h1β1/(β1+P*Y*)2-g1/P* and Q2=[2b3S*-b3+1+b3u1]S*/(1+b3S*).By solving (64), it can be obtained that
(65)λ1(T)=d3(1-S*)Q1(δ+Q2)(1+b3S*)e-δT.By virtue of (62), it gives that
(66)dλ2(T)dT=S[(1+c1)-b4P(1+b3S+b4I)2]λ1+[a3P(1+b3S)(1+b3S+b4I)2+e1-c1S]λ2+e-δT(g1-h1β1P(β1+PY)2)d4+(b3d4-b4d3)S(1+b3S+b4I)2.By associating withM*(S*,0,P*,Y*), the above equation can be deduced as follows:
(67)dλ2(T)dT=Q3λ1+Q5λ2+(S*+u1-1)[d4+(b3d4-b4d3)S*]Q11+b3S*e-δT,
where Q3=[(b4+b3+b3c1)S*-b4+c1+1+b4u1]S*/(1+b3S*) and Q5=[a3+e1-(c1+a3)S*].According to (65), (67) can be rewritten as follows:
(68)dλ2(T)dT=Q5λ2-(S*-1)Q1Q41+b3S*e-δT,
where Q4=d3Q3/(δ+Q2)(1+b3S*)-[d4+(b3d4-b4d3)S*].By solving (68), it can be obtained that
(69)λ2(T)=(S*+u1-1)Q1Q4(1+b3S*)(δ+Q5)e-δT.Based on (65) and (69), it can be obtained that
(70)dλ3(T)dT=h1β1Ye-δT(β1+PY)2+Sλ11+b3S+b4I+a3Iλ21+b3S+b4I+(f1+Y)λ3.By associating withM*(S*,0,P*,Y*), (70) can be deduced as follows:
(71)dλ3(T)dT=h1β1Y*e-δT(β1+P*Y*)2+S*λ11+b3S*+d3S*1+b3S*λ3=d3S*1+b3S*λ3+[h1β1Y*(β1+P*Y*)2+d3S*(1-S*)Q1(δ+Q2)(1+b3S*)2]e-δT.By solving the above equation, it can be obtained that(72)λ3(T)=h1β1d3S*(S*-1)-h1β1Y*(δ+Q2)(1+b3S*)2(δ+Q2)(1+b3S*)(β1+P*Y*)2[δ+(b3δ+d3)S*]×e-δT.In order to obtain an optimal equilibrium solution, it can be obtained by considering the disease-free equilibriumM* and substituting (72) into (62), and it derives that
(73)1-g1(β+P*Y*)2P*h1β1=d3S*(S*-1)-Y*(δ+Q2)(1+b3S*)2(δ+Q2)(1+b3S*)[δ+(b3δ+d3)S*].Equation (73) provides an equation to the singular path and gives the optimal equilibrium levels of population S*=Sδ,  P*=(1-u1-Sδ)(1+b3Sδ), Y*=((d3-b3f1)Sδ-f1)/(1+b3Sδ). Then, the optimal equilibrium levels of taxation and vaccination control can be obtained as follows:
(74)u1δ=(4b3d3Sδ(Sδ-1)e-δTBb3-1-(b3δ+2b3)Sδ+((4b3d3Sδ(Sδ-1)e-δTB)1/2[b3-1-(b3δ+2b3)Sδ]2+(4b3d3Sδ(Sδ-1)e-δTB)1/2)×(2b3)-1,σδ=([g1(d3-b3f1)-h1b3]Sδ2+h1(1-u1)+g1[f1(1-u1)-β1][g1(d3-b3f1)-h1b3]Sδ2)×((1+b3Sδ)(1-u1-Sδ))-1+([(b3-b3u1-1)h1-g1[f1+(1-u1)mmmmmmmmmmmmmn×(d3-b3f1)]]Sδ)×((1+b3Sδ)(1-u1-Sδ))-1.Remark 12 (according to [13]).
λ
i
(
T
)
e
δ
T
(
i
=
1,2
,
3,4
) represent unusual shadow prices along the singular path. From (65), (69), and (72), 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 ∞ [34]. Furthermore, they remain bounded as T→∞.Generally, there is a strong ecological relationship between predator population and susceptible prey population in the harvested eco-epidemiological system, and these two species are usually influenced by each other in several aspects. The imposition of taxation and vaccination defined in (74) act as a deterrent to regulate harvest effort on predator population and control the spread of infectious disease within eco-epidemiological prey predator system, which can not only guarantee the sustainable development of predator population, but also benefit the ecological balance in the ecosystem and sustainable development of other relevant species in the harvested eco-epidemiological system.


## 5.2. Numerical Simulation

With the help of MATLAB, simulation work with hypothetical set of parameters is performed to understand the theoretical results which have been established in Section5.1 of this paper. The values of parameters are cited from [9, 12], which are given as follows: r=1,  k=1,  a1=1,  b1=2.34,  b2=3,  c=0.01,  a2=0.344,  e=0.051,  d1=0.21,  d2=0.6774,  f=0.001,  α=0.4,  h=1.67,  β=1.1,  g=0.41,  δ=0.05,  and   B=1.According to the parameters given above, it follows from (22) that the range of the taxation σ1∈[0.4496,1.6336] and vaccination u1∈[-1.3663,1), which guarantees that there exists a unique disease-free equilibrium. Based on (74), it can be computed that the optimal taxation and optimal vaccination control is σδ=1.0653 and u1δ=0.483, respectively. Since σδ∈[0.4496,1.6336] and u1δ∈[-1.3663,1), it can be obtained that there exists a unique disease-free equilibrium M* corresponding to the optimal taxation and vaccination control, that is, M*(0.5648,0,0.9406,0.0501). According to Theorem 10, model system (20) is globally stable around M*, whose dynamical responses can be shown in Figure 1, and the corresponding phase portraits are shown in Figure 2 (S-I-P space), Figure 3 (S-Y-P space), and Figure 4 (Y-P-I space) with different initial values.Figure 1
Dynamical responses of model system (20) with the optimal taxation and vaccination σδ=1.0653 and u1δ=0.483.Figure 2
Phase portrait (S-I-P space) of model system (20) with σδ=1.0653, and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).Figure 3
Phase portrait (S-Y-P space) of model system (20) with σδ=1.0653 and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).Figure 4
Phase portrait (Y-P-I space) of model system (20) with σδ=1.0653 and u1δ=0.483 and differential initial values, which indicates model system (20) is globally stable around M*(0.5648,0,0.9406,0.0501).


## 6. Conclusion

It is well known that studies on ecology and epidemiology share some common features. Nowadays, taxation is extensively used as economic control instrument to protect the population of the prey predator system from overexploitation, and vaccination is widely adopted to control spread of infectious disease within harvested eco-epidemiological system. Based on superiority of taxation control in regulating harvesting and effectiveness of vaccination in protecting population from infectious disease, it motives the introduction of taxation and vaccination control into model proposed in this paper. In this paper, a harvested eco-epidemiological model is established to investigate the dynamics effects of vaccination and taxation control on the population dynamics. We extend the work in [9] by incorporating the harvest effort into eco-epidemiological prey predator system. Vaccination and taxation are chosen as controlling instruments to protect population from infectious disease and regulate harvesting on predator, respectively.Since it is of more biological significance to consider the effect of interacting species when it comes to discussing dynamics of epidemiological model, conditions that influence nonnegativeness and boundedness of solutions are discussed, which can be found in Theorems3 and 4. Local stability analysis around disease-free equilibrium is studied in Theorem 7. Based on robust Bendixson criterion, global stability of model system around the unique disease-free equilibrium is discussed in Theorem 10, which is theoretically beneficial to studying the coexistence and interaction mechanism of harvested eco-epidemiological prey predator system. By using Pontryagin’s maximum principle, optimal vaccination and taxation control are derived to ensure the sustainable development of biological resource and protect population from infectious disease in the fifth section of this paper. It reveals that the total user’s cost of harvest per unit effort must be equal to the discounted value of the future price at the steady-state level, and infectious disease among prey population can be removed with optimal vaccination strategy on susceptible prey population.It should be noted that the optimal control is obtained based on an interesting reflection, which discusses the impact and efficiency of optimal control and deducts a tendency for the management of biological resource. Since infectious disease spreads within the eco-epidemiological prey predator system and the commercial harvesting exists extensively in the natural world, it is necessary to discuss optimal strategy associated with net harvesting revenue and vaccination cost. Hence, theoretical analysis obtained in the fifth section is constructive for regulatory agency to formulate constructive policies and maintain the sustainable harvesting and effective vaccination by means of optimal taxation and vaccination control; it makes the work studied in this paper have some new and positive features.



---

*Source: 950396-2013-07-30.xml*