# Dynamical Behavior and Stability Analysis in a Hybrid Epidemiological-Economic Model with Incubation

Authors: Chao Liu; Wenquan Yue; Peiyong Liu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A hybrid SIR vector disease model with incubation is established, where susceptible host population satisfies the logistic equation and the recovered host individuals are commercially harvested. It is utilized to discuss the transmission mechanism of infectious disease and dynamical effect of commercial harvest on population dynamics. Positivity and permanence of solutions are analytically investigated. By choosing economic interest of commercial harvesting as a parameter, dynamical behavior and local stability of model system without time delay are studied. It reveals that there is a phenomenon of singularity induced bifurcation as well as local stability switch around interior equilibrium when economic interest increases through zero. State feedback controllers are designed to stabilize model system around the desired interior equilibria in the case of zero economic interest and positive economic interest, respectively. By analyzing corresponding characteristic equation of model system with time delay, local stability analysis around interior equilibrium is discussed due to variation of time delay. Hopf bifurcation occurs at the critical value of time delay and corresponding limit cycle is also observed. Furthermore, directions of Hopf bifurcation and stability of the bifurcating periodic solutions are studied. Numerical simulations are carried out to show consistency with theoretical analysis.


---

## Body

## 1. Introduction

In recent decades, plenty of mathematical models describing the population dynamics of infectious disease have been extensively utilized to understand the transmission mechanism of infectious disease within population ecosystem (see [1–4] and references therein). Much research efforts have been paid to susceptible-infective-recovered (SIR) vector disease model and corresponding model dynamics (see [5–12] and references therein). Generally, in modelling of communicable disease, the incidence rate (the rate of new infections) is considered to play a vital role in ensuring that the model can provide a reasonable qualitative description of the infectious disease dynamics [3, 4].In order to discuss the spread of an infectious disease transmitted by a vector (e.g., mosquitoes and rats), Takeuchi et al. [7] formulated a delayed SIR epidemic model with a bilinear incidence rate. Beretta et al. [8] considered the global stability of disease free equilibrium and endemic equilibrium of model system; it was shown that the disease free equilibrium is globally stable for any time delay while the endemic equilibrium is not feasible. By constructing a suitable Lyapunov functional, sufficient conditions were derived to guarantee that if the endemic equilibrium is feasible, it is also globally stable for the delay being sufficiently small. Ruan and Wang [13] studied the global dynamics of an SIR model with vital dynamics and nonlinear incidence rate of saturated mass action and global qualitative and bifurcation analyses are carried out. Ma et al. [14] derived an explicit expression of lower bound of the infective individual of solution of model system, which was proposed as an open problem. They therefore gave an estimation of the length of time delay ensuring global asymptotic stability of the endemic equilibrium. Xu and Ma [15] proposed an SIR epidemic model with nonlinear incidence rate and time delay. By analyzing the corresponding characteristic equations, local stability of an endemic equilibrium and a disease free equilibrium are discussed. An SIR model with distributed delay and a general incidence function is studied in McCluskey [9], and the global dynamics for the SIR epidemiological system is analyzed in Zhou and Cui [10]. Wang et al. [11] considered the asymptotic behavior of the following SIR vector model:
(1)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t),
where S(t), I(t), and R(t) represent the population density of susceptible, infective, and recovered host individuals at time t, respectively. It is assumed that the population growth of susceptible host individuals is governed by the logistic growth with a carrying capacity k>0 as well as intrinsic birth rate constant r>0. β>0 is the average number of constants per infective per unit time and τ≥0 denotes the incubation time, and μ1>0 and μ2>0 stand for the death rate of infective and recovered host individuals, respectively. m>0 represents the recovery rate of infective host individuals. The local stability of endemic equilibrium is investigated, and conditions for Hopf bifurcation to occur are derived in [11]. Along with the line of this research, Enatsu et al. [12] analyze stability of equilibria for a delayed SIR epidemic model, in which population growth is subject to logistic growth in absence of disease and the proposed model with a nonlinear incidence rate satisfying suitable monotonicity conditions.Nowadays, biological resource within ecosystem is commercially harvested and sold with aim of achieving economic interest [16, 17]. It is well known that harvesting has a strong impact on the dynamic evolution of population and several mathematical models have been established to discuss dynamic effects of harvest effort on population in ecological-epidemiological system, which can be found in [18–21] and the references therein. The role of harvesting in a predator-prey-parasite system is discussed in [18]; theoretical results show that, using impulsive harvesting effort as control parameter, it is not only possible to control the cyclic behavior of the system populations leading to the persistence of all species but other desired stable equilibrium including disease free can be obtained. A ratio-dependent eco-epidemiological system is proposed in [19] where prey population is subject to harvesting. Positive invariance, boundedness, stability of equilibria, and permanence of system have been established. In [20], an eco-epidemiological model is studied where prey disease is modeled by a susceptible-infective scheme, and the role of harvesting and switching on the dynamics of disease propagation and/or eradication is discussed. An eco-epidemiological model with distributed time delay and impulsive control strategy is investigated in [21]; local stability and complex dynamical behavior are discussed. Under the system of market economy, harvest effort is usually influenced by variation of economic interest of commercial harvesting [16, 22]. It should be noted that the above mentioned related work [18–21] only concentrate on the role of harvest effort on population dynamics, while the dynamic effect of economic interest on commercial harvesting and indirect dynamic effect on ecosystem are not considered. The work done in [12] is an extension of [11] with nonlinear incidence rate, while dynamic effect of harvest effort on population dynamics is not considered.Recently, some hybrid dynamical models are proposed in [23–28], which are utilized to discuss the interaction mechanism of harvested ecosystem from an economic perspective. Compared with the traditional mathematical models (differential equations or difference equations) discussing the population dynamics in ecosystem, the hybrid mathematical models proposed in [23–28] are made up of differential equations and algebraic equations, where differential equations concentrate on coexistence and interaction mechanism of population and algebraic equations offer a simpler way to study the effect of harvest effort on ecosystem from an economic perspective. Complex dynamical behavior and stability analysis in prey-predator ecosystems with stage-structured population and gestation delay are considered in [23–28]. In general, differential-algebraic models exhibit more complicated dynamics than ordinary differential models. The differential-algebraic models have been applied widely in power systems, aerospace engineering, chemical processes, social management systems, biological systems, network analysis and oil catalysis, and cracking process (see [29–31] and references therein). With the help of differential-algebraic model for the power systems and bifurcation theory, complex dynamical behaviors of the power systems, especially the bifurcation phenomena that reveal the instability mechanism of power systems have been extensively studied, which can be found in [32–34] and the references therein. Furthermore, some applications of differential-algebraic models in the field of economy, which can be found in [35, 36].It is well known that the recovered host individuals are naturally immune to vector disease [1], and its potential economic interest can be commercially exploited. Furthermore, harvest effort is usually influenced by variation of economic interest of commercial harvesting [16, 22] under the system of market economy. Consequently, it is necessary to discuss the coexistence and interaction mechanism of population within harvested epidemiological ecosystem as well as dynamical effect of harvest effort due to variation of economic interest. However, as far as knowledge goes, nobody has explicitly proposed a mathematical model to discuss the dynamic effect of commercial harvest on epidemiological system under the system of market economy. The main objective of this paper is to investigate the transmission mechanism of infectious disease and dynamical effect of commercial harvest on population dynamics, especially the complex dynamical behavior and stability switch due to variation of incubation and commercial harvest economic interest. The organisation of the rest section of this paper is as follows. By introducing commercial harvest effort into model system (1), a hybrid epidemiological-economic model is established in Section 2. Positivity and permanence of solutions of model system are discussed in Section 3. In Section 4, qualitative analyses of model system are performed. Conditions for existence of interior equilibrium of model system are studied. Dynamical behavior of model system without incubation around the interior equilibrium is investigated due to variation of economic interest, and state feedback controllers are designed to stabilize model system around the desired interior equilibria. Furthermore, local stability analysis of model system with incubation is analyzed due to variation of time delay; directions of Hopf bifurcation and stability of the bifurcating periodic solutions are also studied. Numerical simulations are made in Section 5, which are utilized to support the theoretical findings obtained in this paper. Finally, this paper ends with a conclusion.


## 2. Model Formulation

In 1954, Gordon [22] proposed the economic theory of a common-property resource, which studies the effect of harvest effort on ecosystem from an economic perspective. In [22], an algebraic equation is proposed to investigate the economic interest of yield of the harvest effort, which takes form as follows:
(2)NetEconomicRevenue(NER)=TotalRevenue(TR)-TotalCost(TC).Associated with model (1), an algebraic equation, which considers the economic interest v of the harvest effort on recovered host individuals in epidemiological system, that is, R(t), is established as follows:
(3)E(t)(wR(t)-c)=v,
where E(t) represents the harvest effort on recovered host individuals at time t. v represents the economic interest of harvest effort on the recovered host individuals. w and c represent unit price of harvested population and cost of harvest effort, respectively.Based on (1) and (3), a delayed hybrid model which consists of three differential equations and an algebraic equation can be established as follows:
(4)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-E(t)R(t),0=E(t)(wR(t)-c)-v,
where S(t), I(t), R(t), E(t), and other parameters share the same interpretations mentioned in (1) and (3), and initial conditions ψ=(ψ1,ψ2,ψ3,ψ4) for model system (4) are defined in the Banach space:
(5){ψ∈C+([-τ,0],ℝ+4)∣ψ1(θ)=S(θ),ψ2(θ)=I(θ),mmmmmmmmmmmmψ3(θ)=R(θ),ψ4(θ)=E(θ)C+([-τ,0],ℝ+4)},
where ℝ+4={(S,I,R,E)∈ℝ4:S≥0,I≥0,R≥0,E≥0}. It is also assumed that ψi(0)>0 (i=1,2,3,4) for a biological reason.Model system (4) can be expressed in the following form:
(6)Ξ(t)X˙(t)=F(X(t)),
where
(7)X(t)=(S(t),I(t),R(t),E(t))T,Ξ(t)=[1000010000100000],F(X(t))=[F1(X(t))F2(X(t))F3(X(t))F4(X(t))]=[r(1-S(t)k)-βS(t)I(t-τ)βS(t)I(t-τ)-μ1I(t)-mI(t)mI(t)-μ2R(t)-E(t)R(t)E(t)(wR(t)-c)-v].Remark 1.
The algebraic equation in model system (6) contains no differentiated variables; hence, the leading matrix Ξ(t) in model system (6) has a corresponding zero row.Remark 2.
The model proposed in [11], which composed of differential equations, only discusses the interaction and coexistence mechanism of susceptible, infective, and recovered host individuals. Compared with the model proposed in [11], algebraic equations are incorporated into the model system (4), which focus on the economic interest of harvesting on recovered host individuals. Hence, the established model not only investigates interaction and coexistence mechanism of population in harvested ecosystem but also studies the dynamical behavior due to the variation of economic interest of commercial harvesting and incubation.


## 3. Positivity and Permanence

Theorem 3.
Any solutions of model system (4) with initial conditions are positive.Proof.
For any solutions of model system (4), it is easy to show that Fi:ℝ+4+1→ℝ4 is locally Lipschitz and satisfies the condition, Fi(X(t))|X∈ℝ4>0, where Fi(X(t)) (i=1,2,3,4) have been defined in model system (4).
Due to the lemma in [37] and Theorem A.4 in [38], any solution of the model system (4) with positive initial conditions exists 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 positive. Hence, this completes the positivity of the solutions of model system (4).From a viewpoint of biological and economic interest perspective, persistence of solutions of model system (4) in the case of economic interest v≥0 will be investigated in this section. Some preliminaries are introduced as follows.Definition 4 (see [39]).
Model system (4) is said to be permanent if there exists a compact region Ω0∈int⁡Ω such that every solution of model system (4) with initial conditions will eventually enter and remain in region Ω0.Definition 5 (see [39]).
Consider a metric spaceQ with metric d. The distance d(x,y) of a point x∈Q from a subset Y of Q is defined by
(8)d(x,y)=inf⁡y∈Yd(x,y).It is further assumed thatQ is the closure of an open set Q0, and Q0=∂Q0 is nonempty and is the boundary of Q0. Consequently, Q0∪Q0=Q,Q0∩Q0=Ø. We will also suppose that T(t) is a Q0 semigroup on Q satisfying
(9)T(t):Q0⟶Q0,T(t):Q0⟶Q0.LetT∂(t)=T(t)|Q0 and A∂ be the global attractor for T∂(t).Lemma 6 (see [39]).
Suppose thatT(t) satisfies (9) and the following conditions hold.(i)
There is at0≥0 such that T(t) is compact for t>t0;(ii)
T(t) is point dissipative in Q;(iii)
A~∂=⋃x∈A∂ω(x) is isolated and has an acyclic covering Z.Then T(t) is uniformly persistent if and only if for each Zi∈Z, Ws(Zi)∩Q0=Ø for i=1,2,…,n.Lemma 7 (see [40]).
Consider the following equation:(10)u˙(t)=au(t-τ)-bu(t),
where a,b,τ>0 and u(t)>0 for all -τ≤t≤0; it derives the following:(i)
Ifa<b, then lim⁡t→+∞u(t)=0,(ii)
Ifa>b, then lim⁡t→+∞u(t)=+∞.Lemma 8.
For any solutions of model system (4), we have
(11)limsup⁡t→+∞S(t)≤k,limsup⁡t→+∞(S(t)+I(t))≤k(r+μ1+m)24r(μ1+m).Proof.
By using Theorem3, it follows from the first equation of model system (4) that
(12)S˙(t)≤r(1-S(t)k)S(t),
which derives that limsup⁡t→+∞S(t)≤k.
According to Theorem3 and the first and second equation of model system (4), it gives that
(13)S˙(t)+I˙(t)≤r(1-S(t)k)S(t)-(μ1+m)I(t),
which derives that limsup⁡t→+∞(S(t)+I(t))≤k(r+μ1+m)2/4r(μ1+m).Lemma 9.
Ifμ1+m<1, then (S(t),I(t)) of solution of model system (4) with initial conditions satisfies
(14)liminf⁡t→+∞S(t)≥Sη,liminf⁡t→+∞I(t)≥Iη,
where Sη>0 and Iη>0 are independent of corresponding initial values of model system (4).Proof.
Firstly, letC+([-τ,0],ℝ+2) denote space of continuous functions mapping [-τ,0] into ℝ+2, where ℝ+2={(x,y)∣x≥0,y≥0}:
(15)Q1={(ϕ1,ϕ2)∈C+([-τ,0],ℝ+2)∣ϕ1(θ)=0,θ∈[-τ,0]},Q2={(ϕ1,ϕ2)∈C+([-τ,0],ℝ+2)∣ϕ1(θ)>0,mϕ2(θ)=0,θ∈[-τ,0]C+([-τ,0],ℝ+2)}.
DenoteQ0=Q1∪Q2,  Q=C+([-τ,0],ℝ+2) and Q0=int⁡C+([-τ,0],ℝ+2).
Next, all conditions in Lemma6 will be checked. In order to facilitate the proof, we consider the following subsystem of model system (4):
(16)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),
where S(θ)≥0,  I(θ)≥0 are continuous on θ∈[-τ,0] and S(0)>0,  I(0)>0.
By Definition5 and model system (16), it is easy to see that Q0 and Q0 are positively invariant, and conditions (i) and (ii) of Lemma 6 clearly hold.
Since model system (16) possesses two constant solutions in Q0: P~0∈Q1, P~1∈Q2 with the following form:
(17)P~0={(ϕ1,ϕ2)∈C+([-τ,0],ℝ+2)∣ϕ1(θ)=ϕ2(θ)=0,mmmmmmmmmmmmmmmθ∈[-τ,0]C+([-τ,0],ℝ+2)},P~1={(ϕ1,ϕ2)∈C+([-τ,0],ℝ+2)∣ϕ1(θ)=1,ϕ2(θ)=0,mmmmmmmmmmmmmmmθ∈[-τ,0]C+([-τ,0],ℝ+2)}.
It follows from simple computation that(18)S˙(t)|(ϕ1,ϕ2)∈Q1=0,S(t)S˙(t)|(ϕ1,ϕ2)∈Q1=0fort≥0.
Moreover, it follows from the second equation of model system (16) that
(19)I˙(t)|(ϕ1,ϕ2)∈Q1=-(μ1+m)I(t)≤0,
which reveals that all points in Q1 approach to P~0; that is, Q1=Ws(P~0). By using the similar analysis mentioned above, it can be also concluded that all points in Q2 approach to P~1; that is, Q2=Ws(P~1). Based on the above analysis, it shows that invariant sets P~0 and P~1 are isolated invariant, and {P~0,P~1} is isolated and an acyclic covering. It can be concluded that condition (iii) of Lemma 6 holds.
Finally, we will show thatWs(P~i)∩Q0=Ø for i=0,1. Based on the definition of P~0, it is easy to show that Ws(P~0)∩Q0=Ø. We will show Ws(P~1)∩Q0=Ø in the following part.
IfWs(P~1)∩Q0≠Ø, then there exists a positive solution (S(t),I(t)) to model system (16) with lim⁡t→+∞(S(t),I(t))=(1,0). If μ1+m<1, then μ1+m<1-ϵ holds for sufficiently small ϵ>0 and there exists a positive constant T=T(ϵ) such that S(t)>1-ϵ>0, and 0<I(t)<ϵ for all t≥T.
By the second equation of model system (16), it derives that
(20)I˙(t)≥(1-ϵ)I(t-τ)-(μ1+m)I(t)
holds for all t≥T+τ.
Consider the following equation:(21)x˙(t)=(1-ϵ)x(t-τ)-(μ1+m)x(t),t≥T+τ,x(t)=I(t),T≤t≤T+τ.
Based on (21) and the comparison principle, it derives that I(t)≥x(t) for all t>T.
On the other hand, ifμ1+m<1, then it follows from Lemma 7 that lim⁡t→+∞x(t)=+∞ for all solutions of (21). It can be concluded that lim⁡t→+∞I(t)=∞, which is a contradiction to I(t)<ϵ. Consequently, it can be derived that Ws(P~1)∩Q0=Ø.
According to the above analysis, all conditions of Lemma6 hold. By using Lemma 6, it can be obtained that
(22)liminf⁡t→+∞S(t)≥Sη,liminf⁡t→+∞I(t)≥Iη,
where Sη>0 and Iη>0 are independent of the corresponding initial values of model system (4).Theorem 10.
Ifμ1+m<1, cμ2<mwIη, and 0≤v<cμ2+mwIη, then all solutions of model system (4) with initial conditions are persistent.Proof.
According to Lemmas8 and 9, it can be obtained that
(23)Sη≤S(t)≤k,S(t)+I(t)≤k(r+μ1+m)24r(μ1+m),
hold for all t>0, which derive that
(24)Iη≤I(t)≤k(r+μ1+m)24r(μ1+m)-Sη.
When the economic interestv=0, it follows from Theorem 3 and the fourth equation of model system (4) that
(25)R(t)=cw
and R˙(t)=0. Based on the third equation of model system (4), it can be computed that E(t)=(mw/c)I(t)-μ2. According to (24), it derives that
(26)0<mwIη-cμ2c≤E(t)≤mw[k(r+μ1+m)2-4Sηr(μ1+m)]4cr(μ1+m)-μ2
provided that cμ2<mwIη.
In the case ofv>0, it derives that E(t)=v/(wR(t)-c) based on implicit function theory [41]. According to the third equation of model system (4), it can be obtained that
(27)R˙(t)≥mIη+cμ2-vw-μ2R(t),
which derives that
(28)liminf⁡t→+∞R(t)≥mwIη+cμ2-vwμ2:=R_>0,
provided that 0<v<cμ2+mwIη.
It follows from Theorem3 and the third equation of model system (4) that
(29)R˙(t)≤mI(t)-μ2R(t),
which derives that
(30)limsup⁡t→+∞R(t)≤m[k(r+μ1+m)2-4Sηr(μ1+m)]4rμ2(μ1+m):=R¯.
Hence, it gives thatv/(wR¯-c)≤E(t)≤v/(wR_-c), and it can be rewritten as follows:
(31)4rμ2v(μ1+m)wm[k(r+μ1+m)2-4Sηr(μ1+m)]-4rcμ2(μ1+m)≤E(t)≤μ2vmwIη-v.
Based on (23), (24), (25), and (26), it can be concluded that all solutions of model system (4) with initial conditions are persistent in the case of v=0, and it follows from (23), (24), (28), (30) and (31) that all solutions of model system (4) with initial conditions are persistent in the case of 0<v<cμ2+mwIη.


## 4. Qualitative Analysis of Model System

Dynamical effects of harvest effort and time delay on population dynamics are discussed in this section. It should be noted that the interior equilibrium biologically interprets that susceptible, infective, and recovered host individuals survive as well as harvest on recovered host individuals exists. Bifurcation phenomenon around the interior equilibria can reveal instability mechanism of model system, which are theoretically relevant to infectious disease control and sustainable yield on recovered host individuals in the real world. Consequently, we will mainly concentrate on dynamical behavior and local stability analysis around interior equilibrium of model system (4) in this paper.


### 4.1. Model System without Time Delay

In this section, dynamical behavior of model system (4) without time delay is investigated, and local stability analysis around the interior equilibrium is discussed due to variation of economic interest of commercial harvesting. Furthermore, state feedback controllers are designed to stabilize model system around the desired interior equilibria in the case of zero economic interest and positive economic interest, respectively.


#### 4.1.1. Singularity Induced Bifurcation

Theorem 11.
Model system (4) without time delay has a singularity induced bifurcation around the interior equilibrium, and v=0 is a bifurcation value. Furthermore, local stability switch occurs as v increases through 0.Proof.
Based on the economic theory of a common-property resource [22], there is a phenomenon of bioeconomic equilibrium in the case of zero harvest economic interest; that is, v=0. An interior equilibrium can be obtained as follows: P*(S*,I*,R*,E*), where S*=(μ1+m)/β, I*=r(kβ-μ1-m)/kβ2, R*=c/w, and E*=(wmr(kβ-μ1-m)-ckμ2β2)/ckβ2.
According to biological interpretation of the interior equilibrium, it follows thatS*>0,  I*>0,R*>0 and E*>0. In order to guarantee the existence of interior equilibrium, some inequalities are satisfied:
(32)kβ-μ1-m>0,wmr(kβ-μ1-m)-ckμ2β2>0.
Letv be a bifurcation parameter, H(t)=(S(t),I(t),R(t))T,
(33)h1(H(t),E(t),v)=[r(1-S(t)k)-βS(t)I(t-τ),βS(t)I(t-τ)-μ1I(t)-mI(t),mI(t)-μ2R(t)-E(t)R(t)],h2(H(t),E(t),v)=E(t)(wR(t)-c)-v.
It can be calculated that(34)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=-wmr(kβ-μ1-m)-ckμ2β2kβ2.
By virtue of (32), it can be obtained that
(35)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0.
Furthermore, it can be also calculated that(36)|DHh1DEh1DHh2DEh2|P*=r(μ1+m)(kβ-μ1-m)[wmr(kβ-μ1-m)-ckμ2β2]k2β3.
It follows from (32) that
(37)|DHh1DEh1DHh2DEh2|P*≠0.
Based on Section IV(A) in [42], h3(H(t),E(t),v) can be defined as follows:
(38)h3(H(t),E(t),v)=det⁡(DEg)=wR(t)-c.
By simple computing,(39)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*=cr(μ1+m)(kβ-μ1+m)kβ.
According to (32), it derives that
(40)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*≠0.
Based on the above analysis, four items (i–iv) can be obtained as follows.(i)
It is easy to show thatDEh2 has a simple zero eigenvalue:
(41)h1(H(t),E(t),v)|P*=0,h2(H(t),E(t),v)|P*=0,
and trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0 based on (35).(ii)
It follows from (37) that [DHh1DEh1DHh2DEh2] is nonsingular around P*.(iii)
By virtue of (40), it can be shown that [DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3] is nonsingular around P*; hence rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=5.(iv)
It is easy to showrank⁡(h1(H(t),E(t),v))=3 and rank⁡(h2(H(t),E(t),v))=1, which follows
(42)rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=rank⁡(h1(H(t),E(t),v))+rank⁡(h2(H(t),E(t),v))+1.
It should be noted that the conditions for singularity induced bifurcation, which is introduced in Section III (A) in [42], consist of three conditions, that is, SI1, SI2, and SI3. According to the above items (i)–(iv), SI1, SI2, and SI3 are all satisfied; hence model (4) without time delay has a singularity induced bifurcation around the interior equilibrium P* and the bifurcation value is v=0.
Along with the line of the above proof, for model (4) without time delay, it follows from simple computing that
(43)M=-trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=wmr(kβ-μ1-m)-ckμ2β2kβ2,N=[Dvh3-[DHh3,DEh3][DHh1DEh1DHh2DEh2]-1[Dvh1Dvh2]]|P*=cr(μ1+m)(kβ-μ1-m)kβ.
It follows from (32) that
(44)MN=mrw(kβ-μ1-m)-ckμ2β2crβ(μ1+m)(kβ-μ1-m)>0.
Inequality (44) satisfies Theorem 3 of [42]. According to Theorem 3 of [42], when v increases through 0, one eigenvalue (denoted by λ1) of model system (4) without time delay moves from ℂ- to ℂ+ along the real axis by diverging through infinity; the movement behavior of this eigenvalue influences the stability of model system (4) without time delay.
Since the Jacobian of model system (4) without time delay evaluated around P* takes the following form:
(45)JP*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*0],
according to the leading matrix Ξ(t) in model system (4) and JP*, the characteristic equation of the model system (4) without time delay around P* is
(46)det⁡(λΞ-JP*)=0.
By virtue of simple computation, the characteristic equation is as follows:(47)λ2+rS*kλ+β2S*I*=0.
It can be concluded that the rest eigenvalues of model system (4) without time delay (denoted by λ2 and λ3) have negative real parts by using the Routh-Hurwitz criteria [43]. It follows from Theorem  3 in [42] that there is only one eigenvalue diverging to infinity as v increases through 0, and the rest eigenvalues are continuous and nonzero and cannot jump from one half open complex plane to another one as v increases through 0. It has been shown that λ1 moves from ℂ- to ℂ+ along the real axis by diverging through infinity. Therefore, λ2 and λ3 are continuous and bounded in the ℂ- half plane as v increases through 0 and their movement behaviors have no influence on the stability of model system (4) without time delay around the interior equilibrium P*.According to Table1 and the stability theory, it can be concluded that model system (4) without time delay is stable around P* as v<0 and model system (4) without time delay is unstable around P* as v>0. Consequently, a stability switch occurs as v increases through 0.Table 1
Signs of real parts of eigen values of model (4) without time delay around interior equilibrium P*.
Reλ1
Reλ2
Reλ3
υ
<
0
−
−
−
υ
>
0
+
−
−Remark 12.
Some preliminaries of singularity induced bifurcation are introduced below. Parameter dependent differential-algebraic hybrid system of the form(48)x˙(t)=h(x(t),y(t),λ),h:Rn×Rm×Rp⟶Rn,0=g(x(t),y(t),λ),g:Rn×Rm×Rp⟶Rm,
where x(t), y(t), and λ have appropriate dimensions. It has been shown recently that there are generically three types of codimension one local bifurcation associated with the differential-algebraic model (48), namely, saddle-node bifurcation, Hopf bifurcation, and singularity induced bifurcation (see [42]).The singularity induced bifurcation is firstly introduced and analyzed in [42, 44]. It is a new type of bifurcation and does not occur in usual ordinary differential equation system, which has been characterized for differential-algebraic system, and later improved in [45, 46]. Roughly speaking, the singularity induced bifurcation refers to a stability change of the differential-algebraic hybrid model (48) owing to some eigenvalues of related linearization hx-hygy-1gx diverging to infinity when Jacobian gy is singular.One of the important consequences of the singularity induced bifurcation is that it leads to an impulse phenomenon, which may result in the collapse of the differential-algebraic system (see [45]). More detailed introductions of the singularity induced bifurcation can be found in [42, 44–46].Remark 13.
It follows from Theorem11 that there is a phenomenon of singularity induced bifurcation around the interior equilibrium when economic interest increases through zero, which can cause local stability switch of model system (4). As stated in Remark 12, the singularity induced bifurcation can result in impulse phenomenon, which may lead to the collapse of the proposed model. In the harvested epidemiological-economic system, the impulse phenomenon is vividly reflected with the outbreak of infectious disease during a short period in the real world. Under this climate, the infected population will be beyond the carrying capacity of environment, which is disastrous for sustainable development of the harvested ecosystem as well as prosperous yield on recovered host individuals.


#### 4.1.2. State Feedback Controller

In order to maintain the sustainable yield on recovered host individuals biological resource as well as economic interest of commercial harvesting at an ideal level, some corresponding control strategies should be taken to eliminate the impulse phenomenon caused by singularity induced bifurcation and stabilize model (4) without time delay. In this subsection, state feedback controllers are designed to stabilize model system (4) without time delay around corresponding interior equilibria in the case of v=0 and v>0, respectively.According to the leading matrixΞ(t) in model system (4) and JP* in (45) (the Jacobian of model system (4) without time delay around the interior equilibrium P*), it can be calculated that rank⁡(JP*,ΞJP*,Ξ2JP*,Ξ3JP*)=4. By using Theorem 2-2.1 in [47], it is easy to show that the model system (4) without time delay is locally controllable around the interior equilibrium P* in the case of v=0. Consequently, a state feedback controller can be applied to stabilize the model system (4) without time delay around P*. By using Theorem 3-1.2 in [47], a state feedback controller u(t)=l(E(t)-E*) (l is a feedback gain and E* is the component of the interior equilibrium P*) can be applied to stabilize model system (4) without time delay around P*.Furthermore, the controlled model system (4) without time delay takes the following form:
(49)S˙(t)=r(1-S(t)k)-βS(t)I(t),I˙(t)=βS(t)I(t)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-E(t)R(t),0=E(t)(wR(t)-c)-v+u(t).Theorem 14.
When economic interest of harvesting is zero,v=0, if the feedback gain l satisfies the following inequality:
(50)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mn+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmn+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)},
then singularity induced bifurcation is eliminated and model system (49) is stable around P*.Proof.
The Jacobian of the model system (49) evaluated at the interior equilibrium P* takes the form
(51)J~P*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*l].
According to the leading matrixΞ(t) in the model system (4) and J~P*, the characteristic equation of model system (49) around P* is det⁡(λΞ-J~P*)=0, which can be expressed as follows:
(52)λ3+B1λ2+B2λ+B3=0,
where B1=μ2+E*+rS*/k-wE*R*/l, B2=β2I*S*+(rS*/k)(μ2+E*-wE*R*/l), and B3=β2I*S*(μ2+E*-wE*R*/l).
By using the Routh-Hurwitz criteria [43], model (49) is locally stable around P* if and only if l satisfies
(53)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*24k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mm+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmm+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)}.
Consequently, if the feedback gain satisfies the above inequality, then model system (49) is stable around P* in the case of zero interest of harvesting.LetP~*(S~*,I~*,R~*,E~*) denote interior equilibrium of model (4) in the case of positive economic interest of harvesting (v>0), where S~*=(μ1+m)/β, I~*=r(kβ-μ1-m)/kβ2, E~*=v/(wR*-c), and R~* satisfies the following equation:
(54)R~*2+B~2R~*+B~3=0,
where B~2=(kβ2(v-cμ2)+wmr(m+μ1-kβ))/kwμ2β2, B~3=cmr(kβ-mu1-m)/kwμ2β2.Based on Routh-Hurwitz criteria [43], (54) has two positive roots if economic interest v satisfies the following inequalities:
(55)0<v<min⁡{cμ2+wmr(kβ-m-μ1)kβ2,mmm(1-c)cμ2+wmr(c+1)(kβ-μ1-m)kβ2}∶=v~.As analyzed above, there are two interior equilibria (denoted byP~1* and P~2*) when 0<v*<v~. In this subsection, we only design the controller for the model (4) around the interior equilibrium P~1*. Some symmetric results about P~2* can be also obtained, and P~1* is denoted as P~* for simplicity in the following part.Theorem 15.
When economic interest of harvesting is positive,0<v*<v~, if feedback gain l~ of controller u(t)=l~(E(t)-E~*) satisfies following inequality:
(56)l~>max⁡{kwE~*R~*rS~*+k(μ2+E~*),wE~*R~*μ2+E~*,rwE~*R~*kβ2I~*+r(μ2+E~*),m(+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2wE~*R~*[2kr(μ2+E~*)+2k2β2I~*+r2S~*]mm+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2)m×(+k(μ2+E~*)2]2r[(μ2+E~*)2r(μ2+E~*)S~*+kβ2I~*S~*mmmmm+k(μ2+E~*)2])kwE~*R~*rS~*+k(μ2+E~*)},
then model system (49) is stable around the interior equilibrium P~*(S~*,I~*,R~*,E~*).Proof.
The proof is similar to the proof of Theorem14 of this paper.Remark 16.
It follows from (55) and Theorem 15 that economic interest of commercial harvesting should be regulated within certain interval v∈(0,v~), which guarantees the existence of interior equilibrium in the case of positive economic interest. After applying the state feedback controller into model system (4) without time delay, model system can be stabilized around the corresponding interior equilibrium, respectively. The elimination of the singularity induced bifurcation means the harvested epidemiological-economic system restores to ecological balance and avoidance of infectious disease outbreak.


### 4.2. Model System with Time Delay

By analyzing corresponding characteristic equation of model system with time delay, local stability analysis around the interior equilibrium due to variation of time delay is discussed. Conditions for existence of Hopf bifurcation are studied. Furthermore, directions of Hopf bifurcation and stability of periodic solutions are investigated.


#### 4.2.1. Local Stability and Hopf Bifurcation

As analyzed in the above subsection, in the case of time delay and positive economic interest of harvesting0<v*<v~ where v~ is defined in (55), there are two interior equilibria P~1* and P~2* for model system (4) with respect to the positive economic interest v*.In this subsection, we only investigate dynamical behavior of model system (4) around the interior equilibrium P~1*. Some symmetric results about the interior equilibrium P~2* can be also obtained, and P~1* is denoted as P~* for simplicity. According to Jacobian evaluated at the interior equilibrium P~* and the leading matrix Ξ(t) in model system (4), we can obtain the characteristic equation of model system (4) around P~*, which can be expressed as follows:
(57)|λ+rS~*kβS~*e-λτ00-βI~*λ-βS~*e-λτ+μ1+m000-mλ+μ1+E~*R~*00-wE~*-v*E~*|=0.
⟹(58)M(λ)+N(λ)e-λτ=0,
where
(59)M(λ)=λ3+m1λ2+m2λ+m3,N(λ)=n1λ2+n2λ+n3,m1=μ2+(μ1+m)(1+rkβ)-cv*(wR~*-c)2,m2=μ1+mkβ[r(μ1+m+μ2-cv*(wR~*-c)2)mmmmm+kβ(μ2-cv*(wR~*-c)2)],m3=r(μ1+m)2kβ[μ2-cv*(wR~*-c)2],n1=-μ1-m,n2=-(μ1+m)×[μ2+2r(μ1+m)-krβkβ-cv*(wR~*-c)2],n3=-r(μ1+m)[2(μ1+m)kβ-1]×[μ2-cv*(wR~*-c)2].Now substitutingλ=iσ (σ is a positive real number) into (58) and separating the real and imaginary parts, two transcendental equations can be obtained as follows:
(60)σ3-m2σ=n2σcos⁡(στ)-(n3-n1σ2)sin(στ),(61)m1σ2-m3=(n3-n1σ2)cos⁡(στ)+n2σsin(στ).By squaring and adding (60) and (61), it can be calculated that
(62)(n3-n1σ2)2+n22σ2=(m1σ2-m3)2+(σ3-m2σ)2,
⟹(63)σ6+C1σ4+C2σ2+C3=0,
where C1=m12-2m2-n12, C2=m22-2m1m3+2n1n3-n22, C3=m32-n32.According to the values ofCj,(j=1,2,3) and the Routh-Hurwitz criteria [43], a simple assumption that (58) has at least one positive real root σ0 is C3<0, which derives that kβ>3(μ1+m). Hence, under this assumption, (58) will have a pair of purely imaginary roots of the form ±iσ0.By eliminatingsin(στ) from (60) and (61), it can be calculated that the τj* corresponding to σ0 is as follows:
(64)τj*=1σ0×arccos[n2σ02(σ02-m2)+(n3-n1σ02)(m1σ02-m3)(n3-n1σ02)2+(n2σ0)2]+2jπσ0,
where j=0,1,2,….By using Butler’s lemma [48], model system (4) is locally stable around P~* for τ<τ0*. Subsequently, conditions for existence of Hopf bifurcation in [39] are utilized to investigate whether Hopf bifurcation occurs as τ increases through τj*.Theorem 17.
Ifkβ>3(μ1+m), then model system (4) undergoes Hopf bifurcation around the interior equilibrium P~* when τ=τj*,  j=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Proof.
As mentioned above, letλ=iσ0 represent the purely imaginary root of (58). It follows from (58) that |M(iσ0)|=|N(iσ0)|, which determines a set of possible values of σ0.
In the following part, we determine the direction of motion ofλ=iσ0 as τ is varied; namely, we determine
(65)Θ=sign⁡[d(Reλ)dτ]λ=iσ0=sign⁡[Re(dλdτ)-1]λ=iσ0.
By differentiating (58) with respect to τ, it can be obtained that
(66)(dλdτ)-1=3λ2+2m1λ+m2λe-λτ(n1λ2+n2λ+n3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=3λ2+2m1λ+m2-λ(λ3+m1λ2+m2λ+m3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=2λ3+m1λ2-m3-λ2(λ3+m1λ2+m2λ+m3)+n1λ2-n3λ2(n1λ2+n2λ+n3)-τλ.
From (62) and the above equation, it can be obtained that
(67)Θ=sign⁡[Re(dλdτ)-1]λ=iσ0=1σ02sign⁡[(m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)(m1σ02-m3)2+(σ03-m2σ0)2mmmml+(n1σ02+n3)(n3-n1σ02)(n3-n1σ02)2+(n2σ0)2]=sign⁡[((m1σ02-m3)2+(σ03-m2σ0)2)-1((m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)mmmm+(n1σ02+n3)(n3-n1σ02))mmm×((m1σ02-m3)2+(σ03-m2σ0)2)-1]=sign⁡[2σ06+(m12-2m2-n12)σ04+n32-m32(m1σ02-m3)2+(σ03-m2σ0)2],=sign⁡[2σ06+C1σ04-C3(m1σ02-m3)2+(σ03-m2σ0)2].
According to the values ofCj,(j=1,2,3) given in (58) of this paper, it is easy to show that C1=m12-2m2-n12=r2(μ1+m)2/k2β2+[μ2-cv*/(wR~-c)2]2>0.
Furthermore, ifkβ>3(μ1+m), then it can be shown that C3>0. Hence, it can be concluded that 2σ06+C1σ04-C3>0, which derives sign⁡[d(Reλ)/dτ]τ=τj*,σ=σ0>0. Consequently, the transversality condition holds and Hopf bifurcation occurs at σ=σ0, τ=τj*. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Remark 18.
sign
⁡
[
d
(
Re
λ
)
/
d
τ
]
τ
=
τ
j
*
>
0 signifies that there exists at least one eigenvalue with positive real part for τ=τj*, and the conditions for Hopf bifurcation in [39] are also satisfied yielding the required periodic solution.


#### 4.2.2. Properties of Hopf Bifurcation

By using normal theory and center manifold theorem [49], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 4.1.2, when economic interest of harvesting 0<v*<v~ (v¯ is defined in (55)), it follows from implicit function theorem [41] and the fourth equation of model system (4) that E(t)=v*/(wR(t)-c). Furthermore, model system (4) can be transformed into the following form:
(68)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-v*R(t)wR(t)-c.Firstly, some transformations associated with component(S~*,I~*,R~*) of interior equilibrium P~* are given as follows:
(69)y1=S-S~*,y2=I-I~*,y3=R-R~*,y¯i(t)=yi(τt),τ=ρ+τj,ρ∈ℝ=(-∞,+∞).
Then ρ=0 is the Hopf bifurcation value of model system (4). Bars of variables are dropped for simplicity of notations; model system (4) is transformed to a functional differential equation in C=C([-1,0],ℝ3) as
(70)y˙(t)=Lρ(yt)+f(ρ,yt),
where C=C([-1,0],ℝ3) is the Banach space of continuous functions mapping the interval [-τ,0] into ℝ3, y(t)=(y1(t),y2(t),y3(t))T∈ℝ3, yt(θ)=y(t+θ) for θ∈[-τ,0] and Lρ:C→ℝ3, f:ℝ×C→ℝ3 are defined as follows, respectively:
(71)Lρ(ϕ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×(ϕ1(0)ϕ2(0)ϕ3(0))+(τj+ρ)(0-βS~*00βS~*0000)×(ϕ1(-1)ϕ2(-1)ϕ3(-1)),(72)f(ρ,ϕ)=(τj+ρ)(-rkϕ12(0)-βϕ1(0)ϕ2(-1)βϕ1(0)ϕ2(-1)∑j=1∞fj(3)ϕ3j(0)),
where fj(3)=(di(-v*R/(wR-c))/dRi)|(S~*,I~*,R~*), j=1,2,3.It is easy to show thatLρ is a continuous linear function mapping C into ℝ3. According to Riesz representation theorem [40], there exists a 3×3 matrix function η(θ,ρ) of bounded variation for θ∈[-1,0] such that
(73)Lρ(ϕ)=∫-10‍dη(θ,ρ)ϕ(θ),
where ϕ∈C([-1,0],ℝ3).In fact, we can choose(74)η(θ,ρ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×δ(θ)-(τj+ρ)(0-βS~*00βS~*0000)δ(θ+1),
where δ denotes the Dirac delta function.Ifϕ is any given function in C([-1,0],ℝ3) and y(ϕ) is the unique solution of the linearized equation y˙(t)=Lρ(yt) of (70) with initial function ϕ at zero, then the solution operator T~(t):C→C is defined by
(75)T~(t)ϕ=yt(ϕ),t≥0.It follows from Lemma  7.1.1 in [39] that T~(t),  t≥0 is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator Aρ of T~(t), t≥0 is as follows:
(76)Aρ(ϕ)={dϕ(θ)dθ,θ∈[-1,0)∫-10‍dη(ρ,s)ϕ(s),θ=0,
for ϕ∈C1([-1,0],ℝ3), the space of functions mapping the interval [-1,0] into ℝ3 which have continuous first derivative and also define
(77)R(ρ)(ϕ)={0,θ∈[-1,0)f(ρ,ϕ),θ=0
then model system (70) is equivalent to
(78)y˙t=A(ρ)yt+R(ρ)yt.Forψ∈C1([0,1],(ℝ3)*), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative, define
(79)A*ϕ(s)={-dψ(s)ds,s∈(0,1]∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(80)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). It follows from the above analysis A(0) and A* are adjoint operators.By virtue of discussion in Section4.2.1, ±iω0τj are eigenvalues of A(0). Hence, they are also eigenvalues of A*. In the following, eigenvectors of A(0) and A* are corresponding to iω0τj and -iω0τj, respectively.Supposeq(θ)=(1,a,b)Teiω0τjθ is the eigenvectors of A(0) corresponding to iω0τj, which derives that A(0)q(θ)=iω0τjq(θ). By using the definition of A(0), (71) and (72), it gives that
(81)(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)q(0)+(0-βS~*00βS~*0000)q(-1)=iω0q(0).Forq(-1)=q(0)e-iω0τj, then it can be obtained that
(82)a=-r(μ1+m)+ikβω0kβ(μ1+m)e-iω0τj,b=m(wR~*-c)2[r(m+μ1)+ikβω0]kβ(μ1+m)[cv*-(μ2+iω0)(wR~*-c)2]e-iω0τj.Similarly, it follows from simple computation that eigenvectorq*(s)=J(1,a*,b*)eiω0τjs of A* is corresponding to -iω0τj, where
(83)a*=-r(μ1+m)-ikβω0kβ(μ1+m)eiω0τj,b*=m(wR~*-c)2[r(μ1+m)-ikβω0]kβ(μ1+m)[cv*-(μ2-iω0)(wR~*-c)2]eiω0τj.In order to assume〈q*(s),q(θ)〉=1, we need to determine the value of J in the following part.By virtue of (80), it derives that
(84)〈q*(s),q(θ)〉=J¯(1,a¯*,b¯*)(1,a,b)T-∫-10‍∫ξ=0θ‍J¯(1,a*,b*)e-iω0τj(ξ-θ)dη(θ)mmmmmmm×(1,a,b)Teiω0τjξdξ=J¯[∫-10‍1+aa¯*+bb¯*mn-∫-10‍(1,a¯*,b¯*)θeiω0τjθdη(θ)(1,a,b)T]=J¯[1+aa¯*+bb¯*+kb¯*τj(R~*+bS~*)eiω0τj].Hence, we can chooseJ as follows:
(85)J=11+a¯a*+b¯b*+kb*τj(R~*+bS~*)eiω0τj.Next, we will compute the coordinate to describe the centre manifoldC0 at ρ=0. Let yt be the solution of (78) when ρ=0.Define(86)z(t)=〈q*,yt〉,W(t,θ)=yt(θ)-2Re{z(t)q(θ)}.On the center manifoldC0, it derives that
(87)W(t,θ)=W(z(t),z¯(t),θ),
where
(88)W(z(t),z¯(t),θ)=W20(θ)z22+W11(θ)zz¯+W02z¯22+⋯,z and z¯ are local coordinates for center manifold C0 in the direction of q* and q¯*.It is noted thatW is real if yt is real, and we only consider real solutions. For solution yt∈C0 of (78), since ρ=0, it derives that
(89)z˙(t)=iω0τjz+q¯*(0)f(0,W(z,z¯,0)+2Re{zq(θ)})≜iω0τjz+q¯*(0)f0(z,z¯).The above equation can be rewritten as follows:(90)z˙(t)=iω0τjz(t)+g(z,z¯),
where
(91)g(z,z¯)=q¯*(0)f0(z,z¯)=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.It follows from (86) and (88) that
(92)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W20(θ)z22+W11(θ)zz¯+W02(θ)z¯22+(1,a,b)Teiω0τjθz+(1,a¯,b¯)Te-iω0τjθz¯+⋯.By virtue of (72), (91), and (92), it derives that(93)g(z,z¯)=q¯*(0)f0(z,z¯)=q¯*(0)f(0,yt)=τjJ¯(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0))=τjJ¯[∑j=1∞‍-rky1t2(0)-βy1t(0)y2t(-1)mml+a¯*βy1t(0)y2t(-1)+∑j=1∞‍b¯*fj(3)y3tj(0)]=τjJ¯[∑j=1∞‍-rky1t2(0)+(a¯*-1)βy1t(0)y2t(-1)mml+∑j=1∞‍b¯*fj(3)y3tj(0)]=-rkτjJ¯[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mmmmll+W02(1)(0)z¯22+o(|(z,z¯)|3)]2+(a¯*-1)βτjJ¯×[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mm+W02(1)(0)z¯22+o(|z,z¯|3)]×[e-iω0τjz+eiω0τjz¯+W20(2)(-1)z22mm+W11(2)(-1)zz¯+W02(2)(-1)z¯22+o(|z,z¯|3)]+τjJ¯b¯*[f1(3)(z+z¯+W20(3)(0)z22+W11(3)(0)zz¯mmmmmmmml+W02(1)(0)z¯22+o(|z,z¯|3))mmmmn+f2(3)(z+z¯+W20(3)(0)z22mmmmmmmmmm+W11(3)(0)zz¯+W02(1)(0)z¯22mmmmmmmmmm+W02(1)(0)z¯22+o(|z,z¯|3))+⋯].By comparing the coefficients with (91), it gives that
(94)g20=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W20(3)(0)2+f2(3))],g11=τjJ¯[(f1(3)W20(3)(0)2+f2(3))-2rk+2(a¯*-1)βcos⁡ω0τjmml+f1(3)W11(3)(0)+2f2(3)(f1(3)W20(3)(0)2+f2(3))],g02=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W02(3)(0)2+f2(3))],g21=2τjJ¯[W20(3)(-1)+W20(1)(0)eiω0τj2β(a¯*-1)mmll×(W11(2)(-1)+W20(3)(-1)+W20(1)(0)eiω0τj2mmmmml+W11(1)(0)e-iω0τjW20(3)(-1)+W20(1)(0)eiω0τj2)mmll-rk(W20(1)(0)+2W11(1)(0))mmll+f2(3)(W20(3)(0)+2W11(3)(0))W20(3)(-1)+W20(1)(0)eiω0τj2].Sinceg21 is associated with W20(θ) and W11(θ), further attempts should be carried out to compute W20(θ) and W11(θ).By virtue of (78) and (86), we have
(95)W˙=y˙t-z˙q-z¯˙q¯={AW-2Re{q¯*(0)f0q(θ)},θ∈[-1,0)AW-2Re{q¯*(0)f0q(0)+f0},θ=0.≜AW+H(z,z¯,θ),
where
(96)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯.By substituting the corresponding series into (95) and comparing the coefficients, we have
(97)(A-2iω0τj)W20(θ)=-H20(θ),AW11(θ)=-H11(θ),….It follows from (95) that for θ∈[-1,0)(98)H(z,z¯,θ)=-q¯*(0)f0q(θ)-q*(0)f¯0q¯(θ)=-g(z,z¯)q(θ)-g¯(z,z¯)q¯(θ).By comparing coefficients in (96) with those in (94), it derives that
(99)H20(θ)=-g20q(θ)-g¯02q¯(θ),(100)H11(θ)=-g11q(θ)-g¯11q¯(θ).Based on the definition ofA and (97) and (99), it can be obtained that
(101)W˙20(θ)=2iω0τjW20(θ)+g20q(θ)+g¯02q¯(θ).Forq(θ)=(1,a,b)Teiω0τjθ,
(102)W20(θ)=ig20ω0τjq(0)eiω0τjθ+ig¯023ω0τjq¯(0)e-iω0τjθ+G1e2iω0τjθ,
where G1=(G1(1),G1(2),G1(3)) is a constant vector.Similarly, it follows from (97) and (100) that
(103)W11(θ)=-ig11ω0τjq(0)eiω0τjθ+ig¯11ω0τjq¯(0)e-iω0τjθ+G2,
where G2=(G2(1),G2(2),G2(3)) is a constant vector.Subsequently, values ofG1 and G2 should be computed. By using the definition of A and (95), we have
(104)∫-10‍dηw20(θ)=2iω0τjW20(0)-H20(0),(105)∫-10‍dη(θ)W11(θ)=-H11(0),
where η(θ)=η(0,θ). Based on (95), it derives that in the case of θ=0,
(106)H(z,z¯,0)=-2Re{q¯*(0)f0q(0)}+f(0)=-q¯*(0)f0q(0)-q*(0)f¯0q¯(0)+f0=-g(z,z¯)q(0)-g¯(z,z¯)q¯(0)+f0,
which follows that
(107)H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯=-q(0)(g20z22+g11zz¯+g02z¯22+⋯)-q¯(0)(g¯20z¯22+g¯11zz¯+g¯02z22+⋯)+f0.By virtue of (72), it gives that
(108)f0=τk(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0)).By virtue of (86), it can be obtained that
(109)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W(t,θ)+z(t)q(θ)+z¯(t)q¯(θ)=W20(θ)z22+W21(θ)zz¯+z(t)q(θ)+z¯(t)q¯(θ)+⋯.
Then we have
(110)f0=τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3))z2+τk(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0))zz¯+⋯.According to (107) and (110), we have
(111)H20(0)=-g20q(0)-g¯02q¯(0)+2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),(112)H11(0)=-g11q(0)-g¯11q¯(0)+τj(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Sinceiω0τj is the eigenvalue of A(0) and q(0) is the corresponding eigenvector, we obtain that
(113)(iω0τjI-∫-10‍eiω0τjθdη(θ))q(0)=0,(-iω0τjI-∫-10‍e-iω0τjθdη(θ))q¯(0)=0,
where I is identity matrix.By substituting (102) and (111) into (104), it can be obtained that
(114)(2iω0τjI-∫-10‍e2iω0τjθdη(θ))G1=2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),
which can be rewritten as follows:
(115)(2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2)G1=2(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)).Based on Grammar Law [43], G1(1), G1(2), and G1(3) can be obtained as follows:(116)G1(1)=2U1|-rk-βe-iω0τjβS~*e-2iω0τj0βe-iω0τj2iω0+(μ1+m)-βS~*e-2iω0τj0f1(3)W20(3)(0)2+f2(3)-m2iω0+μ2-cv*(wR~*-c)2|,G1(2)=2U1|2iω0+rS~*k-rk-βe-iω0τj0-βI~*βe-iω0τj00f1(3)W20(3)(0)2+f2(3)2iω0+μ2-cv*(wR~*-c)2|,G1(3)=2U1|2iω0+rS~*kβS~*e-2iω0τj-rk-βe-iω0τj-βI~*2iω0+(μ1+m)-βS~*e-2iω0τjβe-iω0τj0-mf1(3)W20(3)(0)2+f2(3)|,
where
(117)U1=|2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2|.Similarly, substituting (103) and (112) into (105), it can be obtained that
(118)(rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2)G2=2(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Based on Grammar Law [43], G2(1), G2(2), and G2(3) can be obtained as follows:
(119)G2(1)=1U2|-2rk-2βcos⁡ω0τjβS~*02βcos⁡ω0τjμ1+m-βS~*02f2(3)+f1(3)W11(3)(0)-mμ2-cv*(wR~*-c)2|,G2(2)=1U2|rS~*k-2rk-2βcos⁡ω0τj0-βI~*2βcos⁡ω0τj002f2(3)+f1(3)W11(3)(0)μ2-cv*(wR~*-c)2|,G2(3)=1U2|rS~*kβS~*-2rk-2βcos⁡ω0τj-βI~*μ1+m-βS~*2βcos⁡ω0τj0-m2f2(3)+f1(3)W11(3)(0)|,
where
(120)U2=|rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2|.It follows from the above computation and analysis thatW20(θ) and W11(θ) can be determined based on (102) and (103).Furthermore, we can computeg21 based on (94). Hence, the following values can be computed as follows:
(121)d1(0)=i2ω0τj(g20g11-2|g11|2-|g02|23)+g212,δ2=-Re{d1(0)}Re{λ′(τj)},γ2=2Re(d1(0)),T2=Im⁡{d1(0)}+δ2Im⁡{λ′(τj)}ω0τj.By using similar arguments in [49], some properties of bifurcating periodic solutions of model (4) in the center manifold at the critical values are discussed in this paper. Based on the analysis in Section 4.2.2 of this paper, the following theorem can be concluded.Theorem 19.
The properties of Hopf bifurcation are determined by values in (121).(i)
δ2 determines directions of Hopf bifurcation: if δ2>0 (δ2<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τj (τ<τj).(ii)
γ2 determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if γ2<0 (γ2>0).(iii)
T2 determines the period of bifurcating periodic solutions: period increases (decreases) if T2>0 (T2<0).


## 4.1. Model System without Time Delay

In this section, dynamical behavior of model system (4) without time delay is investigated, and local stability analysis around the interior equilibrium is discussed due to variation of economic interest of commercial harvesting. Furthermore, state feedback controllers are designed to stabilize model system around the desired interior equilibria in the case of zero economic interest and positive economic interest, respectively.


### 4.1.1. Singularity Induced Bifurcation

Theorem 11.
Model system (4) without time delay has a singularity induced bifurcation around the interior equilibrium, and v=0 is a bifurcation value. Furthermore, local stability switch occurs as v increases through 0.Proof.
Based on the economic theory of a common-property resource [22], there is a phenomenon of bioeconomic equilibrium in the case of zero harvest economic interest; that is, v=0. An interior equilibrium can be obtained as follows: P*(S*,I*,R*,E*), where S*=(μ1+m)/β, I*=r(kβ-μ1-m)/kβ2, R*=c/w, and E*=(wmr(kβ-μ1-m)-ckμ2β2)/ckβ2.
According to biological interpretation of the interior equilibrium, it follows thatS*>0,  I*>0,R*>0 and E*>0. In order to guarantee the existence of interior equilibrium, some inequalities are satisfied:
(32)kβ-μ1-m>0,wmr(kβ-μ1-m)-ckμ2β2>0.
Letv be a bifurcation parameter, H(t)=(S(t),I(t),R(t))T,
(33)h1(H(t),E(t),v)=[r(1-S(t)k)-βS(t)I(t-τ),βS(t)I(t-τ)-μ1I(t)-mI(t),mI(t)-μ2R(t)-E(t)R(t)],h2(H(t),E(t),v)=E(t)(wR(t)-c)-v.
It can be calculated that(34)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=-wmr(kβ-μ1-m)-ckμ2β2kβ2.
By virtue of (32), it can be obtained that
(35)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0.
Furthermore, it can be also calculated that(36)|DHh1DEh1DHh2DEh2|P*=r(μ1+m)(kβ-μ1-m)[wmr(kβ-μ1-m)-ckμ2β2]k2β3.
It follows from (32) that
(37)|DHh1DEh1DHh2DEh2|P*≠0.
Based on Section IV(A) in [42], h3(H(t),E(t),v) can be defined as follows:
(38)h3(H(t),E(t),v)=det⁡(DEg)=wR(t)-c.
By simple computing,(39)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*=cr(μ1+m)(kβ-μ1+m)kβ.
According to (32), it derives that
(40)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*≠0.
Based on the above analysis, four items (i–iv) can be obtained as follows.(i)
It is easy to show thatDEh2 has a simple zero eigenvalue:
(41)h1(H(t),E(t),v)|P*=0,h2(H(t),E(t),v)|P*=0,
and trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0 based on (35).(ii)
It follows from (37) that [DHh1DEh1DHh2DEh2] is nonsingular around P*.(iii)
By virtue of (40), it can be shown that [DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3] is nonsingular around P*; hence rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=5.(iv)
It is easy to showrank⁡(h1(H(t),E(t),v))=3 and rank⁡(h2(H(t),E(t),v))=1, which follows
(42)rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=rank⁡(h1(H(t),E(t),v))+rank⁡(h2(H(t),E(t),v))+1.
It should be noted that the conditions for singularity induced bifurcation, which is introduced in Section III (A) in [42], consist of three conditions, that is, SI1, SI2, and SI3. According to the above items (i)–(iv), SI1, SI2, and SI3 are all satisfied; hence model (4) without time delay has a singularity induced bifurcation around the interior equilibrium P* and the bifurcation value is v=0.
Along with the line of the above proof, for model (4) without time delay, it follows from simple computing that
(43)M=-trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=wmr(kβ-μ1-m)-ckμ2β2kβ2,N=[Dvh3-[DHh3,DEh3][DHh1DEh1DHh2DEh2]-1[Dvh1Dvh2]]|P*=cr(μ1+m)(kβ-μ1-m)kβ.
It follows from (32) that
(44)MN=mrw(kβ-μ1-m)-ckμ2β2crβ(μ1+m)(kβ-μ1-m)>0.
Inequality (44) satisfies Theorem 3 of [42]. According to Theorem 3 of [42], when v increases through 0, one eigenvalue (denoted by λ1) of model system (4) without time delay moves from ℂ- to ℂ+ along the real axis by diverging through infinity; the movement behavior of this eigenvalue influences the stability of model system (4) without time delay.
Since the Jacobian of model system (4) without time delay evaluated around P* takes the following form:
(45)JP*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*0],
according to the leading matrix Ξ(t) in model system (4) and JP*, the characteristic equation of the model system (4) without time delay around P* is
(46)det⁡(λΞ-JP*)=0.
By virtue of simple computation, the characteristic equation is as follows:(47)λ2+rS*kλ+β2S*I*=0.
It can be concluded that the rest eigenvalues of model system (4) without time delay (denoted by λ2 and λ3) have negative real parts by using the Routh-Hurwitz criteria [43]. It follows from Theorem  3 in [42] that there is only one eigenvalue diverging to infinity as v increases through 0, and the rest eigenvalues are continuous and nonzero and cannot jump from one half open complex plane to another one as v increases through 0. It has been shown that λ1 moves from ℂ- to ℂ+ along the real axis by diverging through infinity. Therefore, λ2 and λ3 are continuous and bounded in the ℂ- half plane as v increases through 0 and their movement behaviors have no influence on the stability of model system (4) without time delay around the interior equilibrium P*.According to Table1 and the stability theory, it can be concluded that model system (4) without time delay is stable around P* as v<0 and model system (4) without time delay is unstable around P* as v>0. Consequently, a stability switch occurs as v increases through 0.Table 1
Signs of real parts of eigen values of model (4) without time delay around interior equilibrium P*.
Reλ1
Reλ2
Reλ3
υ
<
0
−
−
−
υ
>
0
+
−
−Remark 12.
Some preliminaries of singularity induced bifurcation are introduced below. Parameter dependent differential-algebraic hybrid system of the form(48)x˙(t)=h(x(t),y(t),λ),h:Rn×Rm×Rp⟶Rn,0=g(x(t),y(t),λ),g:Rn×Rm×Rp⟶Rm,
where x(t), y(t), and λ have appropriate dimensions. It has been shown recently that there are generically three types of codimension one local bifurcation associated with the differential-algebraic model (48), namely, saddle-node bifurcation, Hopf bifurcation, and singularity induced bifurcation (see [42]).The singularity induced bifurcation is firstly introduced and analyzed in [42, 44]. It is a new type of bifurcation and does not occur in usual ordinary differential equation system, which has been characterized for differential-algebraic system, and later improved in [45, 46]. Roughly speaking, the singularity induced bifurcation refers to a stability change of the differential-algebraic hybrid model (48) owing to some eigenvalues of related linearization hx-hygy-1gx diverging to infinity when Jacobian gy is singular.One of the important consequences of the singularity induced bifurcation is that it leads to an impulse phenomenon, which may result in the collapse of the differential-algebraic system (see [45]). More detailed introductions of the singularity induced bifurcation can be found in [42, 44–46].Remark 13.
It follows from Theorem11 that there is a phenomenon of singularity induced bifurcation around the interior equilibrium when economic interest increases through zero, which can cause local stability switch of model system (4). As stated in Remark 12, the singularity induced bifurcation can result in impulse phenomenon, which may lead to the collapse of the proposed model. In the harvested epidemiological-economic system, the impulse phenomenon is vividly reflected with the outbreak of infectious disease during a short period in the real world. Under this climate, the infected population will be beyond the carrying capacity of environment, which is disastrous for sustainable development of the harvested ecosystem as well as prosperous yield on recovered host individuals.


### 4.1.2. State Feedback Controller

In order to maintain the sustainable yield on recovered host individuals biological resource as well as economic interest of commercial harvesting at an ideal level, some corresponding control strategies should be taken to eliminate the impulse phenomenon caused by singularity induced bifurcation and stabilize model (4) without time delay. In this subsection, state feedback controllers are designed to stabilize model system (4) without time delay around corresponding interior equilibria in the case of v=0 and v>0, respectively.According to the leading matrixΞ(t) in model system (4) and JP* in (45) (the Jacobian of model system (4) without time delay around the interior equilibrium P*), it can be calculated that rank⁡(JP*,ΞJP*,Ξ2JP*,Ξ3JP*)=4. By using Theorem 2-2.1 in [47], it is easy to show that the model system (4) without time delay is locally controllable around the interior equilibrium P* in the case of v=0. Consequently, a state feedback controller can be applied to stabilize the model system (4) without time delay around P*. By using Theorem 3-1.2 in [47], a state feedback controller u(t)=l(E(t)-E*) (l is a feedback gain and E* is the component of the interior equilibrium P*) can be applied to stabilize model system (4) without time delay around P*.Furthermore, the controlled model system (4) without time delay takes the following form:
(49)S˙(t)=r(1-S(t)k)-βS(t)I(t),I˙(t)=βS(t)I(t)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-E(t)R(t),0=E(t)(wR(t)-c)-v+u(t).Theorem 14.
When economic interest of harvesting is zero,v=0, if the feedback gain l satisfies the following inequality:
(50)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mn+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmn+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)},
then singularity induced bifurcation is eliminated and model system (49) is stable around P*.Proof.
The Jacobian of the model system (49) evaluated at the interior equilibrium P* takes the form
(51)J~P*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*l].
According to the leading matrixΞ(t) in the model system (4) and J~P*, the characteristic equation of model system (49) around P* is det⁡(λΞ-J~P*)=0, which can be expressed as follows:
(52)λ3+B1λ2+B2λ+B3=0,
where B1=μ2+E*+rS*/k-wE*R*/l, B2=β2I*S*+(rS*/k)(μ2+E*-wE*R*/l), and B3=β2I*S*(μ2+E*-wE*R*/l).
By using the Routh-Hurwitz criteria [43], model (49) is locally stable around P* if and only if l satisfies
(53)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*24k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mm+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmm+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)}.
Consequently, if the feedback gain satisfies the above inequality, then model system (49) is stable around P* in the case of zero interest of harvesting.LetP~*(S~*,I~*,R~*,E~*) denote interior equilibrium of model (4) in the case of positive economic interest of harvesting (v>0), where S~*=(μ1+m)/β, I~*=r(kβ-μ1-m)/kβ2, E~*=v/(wR*-c), and R~* satisfies the following equation:
(54)R~*2+B~2R~*+B~3=0,
where B~2=(kβ2(v-cμ2)+wmr(m+μ1-kβ))/kwμ2β2, B~3=cmr(kβ-mu1-m)/kwμ2β2.Based on Routh-Hurwitz criteria [43], (54) has two positive roots if economic interest v satisfies the following inequalities:
(55)0<v<min⁡{cμ2+wmr(kβ-m-μ1)kβ2,mmm(1-c)cμ2+wmr(c+1)(kβ-μ1-m)kβ2}∶=v~.As analyzed above, there are two interior equilibria (denoted byP~1* and P~2*) when 0<v*<v~. In this subsection, we only design the controller for the model (4) around the interior equilibrium P~1*. Some symmetric results about P~2* can be also obtained, and P~1* is denoted as P~* for simplicity in the following part.Theorem 15.
When economic interest of harvesting is positive,0<v*<v~, if feedback gain l~ of controller u(t)=l~(E(t)-E~*) satisfies following inequality:
(56)l~>max⁡{kwE~*R~*rS~*+k(μ2+E~*),wE~*R~*μ2+E~*,rwE~*R~*kβ2I~*+r(μ2+E~*),m(+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2wE~*R~*[2kr(μ2+E~*)+2k2β2I~*+r2S~*]mm+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2)m×(+k(μ2+E~*)2]2r[(μ2+E~*)2r(μ2+E~*)S~*+kβ2I~*S~*mmmmm+k(μ2+E~*)2])kwE~*R~*rS~*+k(μ2+E~*)},
then model system (49) is stable around the interior equilibrium P~*(S~*,I~*,R~*,E~*).Proof.
The proof is similar to the proof of Theorem14 of this paper.Remark 16.
It follows from (55) and Theorem 15 that economic interest of commercial harvesting should be regulated within certain interval v∈(0,v~), which guarantees the existence of interior equilibrium in the case of positive economic interest. After applying the state feedback controller into model system (4) without time delay, model system can be stabilized around the corresponding interior equilibrium, respectively. The elimination of the singularity induced bifurcation means the harvested epidemiological-economic system restores to ecological balance and avoidance of infectious disease outbreak.


## 4.1.1. Singularity Induced Bifurcation

Theorem 11.
Model system (4) without time delay has a singularity induced bifurcation around the interior equilibrium, and v=0 is a bifurcation value. Furthermore, local stability switch occurs as v increases through 0.Proof.
Based on the economic theory of a common-property resource [22], there is a phenomenon of bioeconomic equilibrium in the case of zero harvest economic interest; that is, v=0. An interior equilibrium can be obtained as follows: P*(S*,I*,R*,E*), where S*=(μ1+m)/β, I*=r(kβ-μ1-m)/kβ2, R*=c/w, and E*=(wmr(kβ-μ1-m)-ckμ2β2)/ckβ2.
According to biological interpretation of the interior equilibrium, it follows thatS*>0,  I*>0,R*>0 and E*>0. In order to guarantee the existence of interior equilibrium, some inequalities are satisfied:
(32)kβ-μ1-m>0,wmr(kβ-μ1-m)-ckμ2β2>0.
Letv be a bifurcation parameter, H(t)=(S(t),I(t),R(t))T,
(33)h1(H(t),E(t),v)=[r(1-S(t)k)-βS(t)I(t-τ),βS(t)I(t-τ)-μ1I(t)-mI(t),mI(t)-μ2R(t)-E(t)R(t)],h2(H(t),E(t),v)=E(t)(wR(t)-c)-v.
It can be calculated that(34)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=-wmr(kβ-μ1-m)-ckμ2β2kβ2.
By virtue of (32), it can be obtained that
(35)trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0.
Furthermore, it can be also calculated that(36)|DHh1DEh1DHh2DEh2|P*=r(μ1+m)(kβ-μ1-m)[wmr(kβ-μ1-m)-ckμ2β2]k2β3.
It follows from (32) that
(37)|DHh1DEh1DHh2DEh2|P*≠0.
Based on Section IV(A) in [42], h3(H(t),E(t),v) can be defined as follows:
(38)h3(H(t),E(t),v)=det⁡(DEg)=wR(t)-c.
By simple computing,(39)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*=cr(μ1+m)(kβ-μ1+m)kβ.
According to (32), it derives that
(40)|DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3|P*≠0.
Based on the above analysis, four items (i–iv) can be obtained as follows.(i)
It is easy to show thatDEh2 has a simple zero eigenvalue:
(41)h1(H(t),E(t),v)|P*=0,h2(H(t),E(t),v)|P*=0,
and trace(DEh1adj(DEh2)(DHh2,DEh2))|P*≠0 based on (35).(ii)
It follows from (37) that [DHh1DEh1DHh2DEh2] is nonsingular around P*.(iii)
By virtue of (40), it can be shown that [DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3] is nonsingular around P*; hence rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=5.(iv)
It is easy to showrank⁡(h1(H(t),E(t),v))=3 and rank⁡(h2(H(t),E(t),v))=1, which follows
(42)rank⁡[DHh1DEh1Dvh1DHh2DEh2Dvh2DHh3DEh3Dvh3]=rank⁡(h1(H(t),E(t),v))+rank⁡(h2(H(t),E(t),v))+1.
It should be noted that the conditions for singularity induced bifurcation, which is introduced in Section III (A) in [42], consist of three conditions, that is, SI1, SI2, and SI3. According to the above items (i)–(iv), SI1, SI2, and SI3 are all satisfied; hence model (4) without time delay has a singularity induced bifurcation around the interior equilibrium P* and the bifurcation value is v=0.
Along with the line of the above proof, for model (4) without time delay, it follows from simple computing that
(43)M=-trace(DEh1adj(DEh2)(DHh2,DEh2))|P*=wmr(kβ-μ1-m)-ckμ2β2kβ2,N=[Dvh3-[DHh3,DEh3][DHh1DEh1DHh2DEh2]-1[Dvh1Dvh2]]|P*=cr(μ1+m)(kβ-μ1-m)kβ.
It follows from (32) that
(44)MN=mrw(kβ-μ1-m)-ckμ2β2crβ(μ1+m)(kβ-μ1-m)>0.
Inequality (44) satisfies Theorem 3 of [42]. According to Theorem 3 of [42], when v increases through 0, one eigenvalue (denoted by λ1) of model system (4) without time delay moves from ℂ- to ℂ+ along the real axis by diverging through infinity; the movement behavior of this eigenvalue influences the stability of model system (4) without time delay.
Since the Jacobian of model system (4) without time delay evaluated around P* takes the following form:
(45)JP*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*0],
according to the leading matrix Ξ(t) in model system (4) and JP*, the characteristic equation of the model system (4) without time delay around P* is
(46)det⁡(λΞ-JP*)=0.
By virtue of simple computation, the characteristic equation is as follows:(47)λ2+rS*kλ+β2S*I*=0.
It can be concluded that the rest eigenvalues of model system (4) without time delay (denoted by λ2 and λ3) have negative real parts by using the Routh-Hurwitz criteria [43]. It follows from Theorem  3 in [42] that there is only one eigenvalue diverging to infinity as v increases through 0, and the rest eigenvalues are continuous and nonzero and cannot jump from one half open complex plane to another one as v increases through 0. It has been shown that λ1 moves from ℂ- to ℂ+ along the real axis by diverging through infinity. Therefore, λ2 and λ3 are continuous and bounded in the ℂ- half plane as v increases through 0 and their movement behaviors have no influence on the stability of model system (4) without time delay around the interior equilibrium P*.According to Table1 and the stability theory, it can be concluded that model system (4) without time delay is stable around P* as v<0 and model system (4) without time delay is unstable around P* as v>0. Consequently, a stability switch occurs as v increases through 0.Table 1
Signs of real parts of eigen values of model (4) without time delay around interior equilibrium P*.
Reλ1
Reλ2
Reλ3
υ
<
0
−
−
−
υ
>
0
+
−
−Remark 12.
Some preliminaries of singularity induced bifurcation are introduced below. Parameter dependent differential-algebraic hybrid system of the form(48)x˙(t)=h(x(t),y(t),λ),h:Rn×Rm×Rp⟶Rn,0=g(x(t),y(t),λ),g:Rn×Rm×Rp⟶Rm,
where x(t), y(t), and λ have appropriate dimensions. It has been shown recently that there are generically three types of codimension one local bifurcation associated with the differential-algebraic model (48), namely, saddle-node bifurcation, Hopf bifurcation, and singularity induced bifurcation (see [42]).The singularity induced bifurcation is firstly introduced and analyzed in [42, 44]. It is a new type of bifurcation and does not occur in usual ordinary differential equation system, which has been characterized for differential-algebraic system, and later improved in [45, 46]. Roughly speaking, the singularity induced bifurcation refers to a stability change of the differential-algebraic hybrid model (48) owing to some eigenvalues of related linearization hx-hygy-1gx diverging to infinity when Jacobian gy is singular.One of the important consequences of the singularity induced bifurcation is that it leads to an impulse phenomenon, which may result in the collapse of the differential-algebraic system (see [45]). More detailed introductions of the singularity induced bifurcation can be found in [42, 44–46].Remark 13.
It follows from Theorem11 that there is a phenomenon of singularity induced bifurcation around the interior equilibrium when economic interest increases through zero, which can cause local stability switch of model system (4). As stated in Remark 12, the singularity induced bifurcation can result in impulse phenomenon, which may lead to the collapse of the proposed model. In the harvested epidemiological-economic system, the impulse phenomenon is vividly reflected with the outbreak of infectious disease during a short period in the real world. Under this climate, the infected population will be beyond the carrying capacity of environment, which is disastrous for sustainable development of the harvested ecosystem as well as prosperous yield on recovered host individuals.


## 4.1.2. State Feedback Controller

In order to maintain the sustainable yield on recovered host individuals biological resource as well as economic interest of commercial harvesting at an ideal level, some corresponding control strategies should be taken to eliminate the impulse phenomenon caused by singularity induced bifurcation and stabilize model (4) without time delay. In this subsection, state feedback controllers are designed to stabilize model system (4) without time delay around corresponding interior equilibria in the case of v=0 and v>0, respectively.According to the leading matrixΞ(t) in model system (4) and JP* in (45) (the Jacobian of model system (4) without time delay around the interior equilibrium P*), it can be calculated that rank⁡(JP*,ΞJP*,Ξ2JP*,Ξ3JP*)=4. By using Theorem 2-2.1 in [47], it is easy to show that the model system (4) without time delay is locally controllable around the interior equilibrium P* in the case of v=0. Consequently, a state feedback controller can be applied to stabilize the model system (4) without time delay around P*. By using Theorem 3-1.2 in [47], a state feedback controller u(t)=l(E(t)-E*) (l is a feedback gain and E* is the component of the interior equilibrium P*) can be applied to stabilize model system (4) without time delay around P*.Furthermore, the controlled model system (4) without time delay takes the following form:
(49)S˙(t)=r(1-S(t)k)-βS(t)I(t),I˙(t)=βS(t)I(t)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-E(t)R(t),0=E(t)(wR(t)-c)-v+u(t).Theorem 14.
When economic interest of harvesting is zero,v=0, if the feedback gain l satisfies the following inequality:
(50)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mn+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmn+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)},
then singularity induced bifurcation is eliminated and model system (49) is stable around P*.Proof.
The Jacobian of the model system (49) evaluated at the interior equilibrium P* takes the form
(51)J~P*=[-rS*k-βS*00βI*0000m-(μ2+E*)-R*00wE*l].
According to the leading matrixΞ(t) in the model system (4) and J~P*, the characteristic equation of model system (49) around P* is det⁡(λΞ-J~P*)=0, which can be expressed as follows:
(52)λ3+B1λ2+B2λ+B3=0,
where B1=μ2+E*+rS*/k-wE*R*/l, B2=β2I*S*+(rS*/k)(μ2+E*-wE*R*/l), and B3=β2I*S*(μ2+E*-wE*R*/l).
By using the Routh-Hurwitz criteria [43], model (49) is locally stable around P* if and only if l satisfies
(53)l>max⁡{kwE*R*rS*+k(μ2+E*),wE*R*μ2+E*,rwE*R*kβ2I*+r(μ2+E*),m(4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*24k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2wE*R*[2kr(μ2+E*)+2k2β2I*+r2S*]mm+4k4β4I*2+8(μ2+E*)rk3β2I*+r4S*2)m×(2r[r(μ2+E*)S*+kβ2I*S*mmmmm+k(μ2+E*)2])-1kwE*R*rS*+k(μ2+E*)}.
Consequently, if the feedback gain satisfies the above inequality, then model system (49) is stable around P* in the case of zero interest of harvesting.LetP~*(S~*,I~*,R~*,E~*) denote interior equilibrium of model (4) in the case of positive economic interest of harvesting (v>0), where S~*=(μ1+m)/β, I~*=r(kβ-μ1-m)/kβ2, E~*=v/(wR*-c), and R~* satisfies the following equation:
(54)R~*2+B~2R~*+B~3=0,
where B~2=(kβ2(v-cμ2)+wmr(m+μ1-kβ))/kwμ2β2, B~3=cmr(kβ-mu1-m)/kwμ2β2.Based on Routh-Hurwitz criteria [43], (54) has two positive roots if economic interest v satisfies the following inequalities:
(55)0<v<min⁡{cμ2+wmr(kβ-m-μ1)kβ2,mmm(1-c)cμ2+wmr(c+1)(kβ-μ1-m)kβ2}∶=v~.As analyzed above, there are two interior equilibria (denoted byP~1* and P~2*) when 0<v*<v~. In this subsection, we only design the controller for the model (4) around the interior equilibrium P~1*. Some symmetric results about P~2* can be also obtained, and P~1* is denoted as P~* for simplicity in the following part.Theorem 15.
When economic interest of harvesting is positive,0<v*<v~, if feedback gain l~ of controller u(t)=l~(E(t)-E~*) satisfies following inequality:
(56)l~>max⁡{kwE~*R~*rS~*+k(μ2+E~*),wE~*R~*μ2+E~*,rwE~*R~*kβ2I~*+r(μ2+E~*),m(+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2wE~*R~*[2kr(μ2+E~*)+2k2β2I~*+r2S~*]mm+4k4β4I~*2+8(μ2+E~*)rk3β2I~*+r4S~*2)m×(+k(μ2+E~*)2]2r[(μ2+E~*)2r(μ2+E~*)S~*+kβ2I~*S~*mmmmm+k(μ2+E~*)2])kwE~*R~*rS~*+k(μ2+E~*)},
then model system (49) is stable around the interior equilibrium P~*(S~*,I~*,R~*,E~*).Proof.
The proof is similar to the proof of Theorem14 of this paper.Remark 16.
It follows from (55) and Theorem 15 that economic interest of commercial harvesting should be regulated within certain interval v∈(0,v~), which guarantees the existence of interior equilibrium in the case of positive economic interest. After applying the state feedback controller into model system (4) without time delay, model system can be stabilized around the corresponding interior equilibrium, respectively. The elimination of the singularity induced bifurcation means the harvested epidemiological-economic system restores to ecological balance and avoidance of infectious disease outbreak.


## 4.2. Model System with Time Delay

By analyzing corresponding characteristic equation of model system with time delay, local stability analysis around the interior equilibrium due to variation of time delay is discussed. Conditions for existence of Hopf bifurcation are studied. Furthermore, directions of Hopf bifurcation and stability of periodic solutions are investigated.


### 4.2.1. Local Stability and Hopf Bifurcation

As analyzed in the above subsection, in the case of time delay and positive economic interest of harvesting0<v*<v~ where v~ is defined in (55), there are two interior equilibria P~1* and P~2* for model system (4) with respect to the positive economic interest v*.In this subsection, we only investigate dynamical behavior of model system (4) around the interior equilibrium P~1*. Some symmetric results about the interior equilibrium P~2* can be also obtained, and P~1* is denoted as P~* for simplicity. According to Jacobian evaluated at the interior equilibrium P~* and the leading matrix Ξ(t) in model system (4), we can obtain the characteristic equation of model system (4) around P~*, which can be expressed as follows:
(57)|λ+rS~*kβS~*e-λτ00-βI~*λ-βS~*e-λτ+μ1+m000-mλ+μ1+E~*R~*00-wE~*-v*E~*|=0.
⟹(58)M(λ)+N(λ)e-λτ=0,
where
(59)M(λ)=λ3+m1λ2+m2λ+m3,N(λ)=n1λ2+n2λ+n3,m1=μ2+(μ1+m)(1+rkβ)-cv*(wR~*-c)2,m2=μ1+mkβ[r(μ1+m+μ2-cv*(wR~*-c)2)mmmmm+kβ(μ2-cv*(wR~*-c)2)],m3=r(μ1+m)2kβ[μ2-cv*(wR~*-c)2],n1=-μ1-m,n2=-(μ1+m)×[μ2+2r(μ1+m)-krβkβ-cv*(wR~*-c)2],n3=-r(μ1+m)[2(μ1+m)kβ-1]×[μ2-cv*(wR~*-c)2].Now substitutingλ=iσ (σ is a positive real number) into (58) and separating the real and imaginary parts, two transcendental equations can be obtained as follows:
(60)σ3-m2σ=n2σcos⁡(στ)-(n3-n1σ2)sin(στ),(61)m1σ2-m3=(n3-n1σ2)cos⁡(στ)+n2σsin(στ).By squaring and adding (60) and (61), it can be calculated that
(62)(n3-n1σ2)2+n22σ2=(m1σ2-m3)2+(σ3-m2σ)2,
⟹(63)σ6+C1σ4+C2σ2+C3=0,
where C1=m12-2m2-n12, C2=m22-2m1m3+2n1n3-n22, C3=m32-n32.According to the values ofCj,(j=1,2,3) and the Routh-Hurwitz criteria [43], a simple assumption that (58) has at least one positive real root σ0 is C3<0, which derives that kβ>3(μ1+m). Hence, under this assumption, (58) will have a pair of purely imaginary roots of the form ±iσ0.By eliminatingsin(στ) from (60) and (61), it can be calculated that the τj* corresponding to σ0 is as follows:
(64)τj*=1σ0×arccos[n2σ02(σ02-m2)+(n3-n1σ02)(m1σ02-m3)(n3-n1σ02)2+(n2σ0)2]+2jπσ0,
where j=0,1,2,….By using Butler’s lemma [48], model system (4) is locally stable around P~* for τ<τ0*. Subsequently, conditions for existence of Hopf bifurcation in [39] are utilized to investigate whether Hopf bifurcation occurs as τ increases through τj*.Theorem 17.
Ifkβ>3(μ1+m), then model system (4) undergoes Hopf bifurcation around the interior equilibrium P~* when τ=τj*,  j=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Proof.
As mentioned above, letλ=iσ0 represent the purely imaginary root of (58). It follows from (58) that |M(iσ0)|=|N(iσ0)|, which determines a set of possible values of σ0.
In the following part, we determine the direction of motion ofλ=iσ0 as τ is varied; namely, we determine
(65)Θ=sign⁡[d(Reλ)dτ]λ=iσ0=sign⁡[Re(dλdτ)-1]λ=iσ0.
By differentiating (58) with respect to τ, it can be obtained that
(66)(dλdτ)-1=3λ2+2m1λ+m2λe-λτ(n1λ2+n2λ+n3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=3λ2+2m1λ+m2-λ(λ3+m1λ2+m2λ+m3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=2λ3+m1λ2-m3-λ2(λ3+m1λ2+m2λ+m3)+n1λ2-n3λ2(n1λ2+n2λ+n3)-τλ.
From (62) and the above equation, it can be obtained that
(67)Θ=sign⁡[Re(dλdτ)-1]λ=iσ0=1σ02sign⁡[(m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)(m1σ02-m3)2+(σ03-m2σ0)2mmmml+(n1σ02+n3)(n3-n1σ02)(n3-n1σ02)2+(n2σ0)2]=sign⁡[((m1σ02-m3)2+(σ03-m2σ0)2)-1((m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)mmmm+(n1σ02+n3)(n3-n1σ02))mmm×((m1σ02-m3)2+(σ03-m2σ0)2)-1]=sign⁡[2σ06+(m12-2m2-n12)σ04+n32-m32(m1σ02-m3)2+(σ03-m2σ0)2],=sign⁡[2σ06+C1σ04-C3(m1σ02-m3)2+(σ03-m2σ0)2].
According to the values ofCj,(j=1,2,3) given in (58) of this paper, it is easy to show that C1=m12-2m2-n12=r2(μ1+m)2/k2β2+[μ2-cv*/(wR~-c)2]2>0.
Furthermore, ifkβ>3(μ1+m), then it can be shown that C3>0. Hence, it can be concluded that 2σ06+C1σ04-C3>0, which derives sign⁡[d(Reλ)/dτ]τ=τj*,σ=σ0>0. Consequently, the transversality condition holds and Hopf bifurcation occurs at σ=σ0, τ=τj*. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Remark 18.
sign
⁡
[
d
(
Re
λ
)
/
d
τ
]
τ
=
τ
j
*
>
0 signifies that there exists at least one eigenvalue with positive real part for τ=τj*, and the conditions for Hopf bifurcation in [39] are also satisfied yielding the required periodic solution.


### 4.2.2. Properties of Hopf Bifurcation

By using normal theory and center manifold theorem [49], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 4.1.2, when economic interest of harvesting 0<v*<v~ (v¯ is defined in (55)), it follows from implicit function theorem [41] and the fourth equation of model system (4) that E(t)=v*/(wR(t)-c). Furthermore, model system (4) can be transformed into the following form:
(68)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-v*R(t)wR(t)-c.Firstly, some transformations associated with component(S~*,I~*,R~*) of interior equilibrium P~* are given as follows:
(69)y1=S-S~*,y2=I-I~*,y3=R-R~*,y¯i(t)=yi(τt),τ=ρ+τj,ρ∈ℝ=(-∞,+∞).
Then ρ=0 is the Hopf bifurcation value of model system (4). Bars of variables are dropped for simplicity of notations; model system (4) is transformed to a functional differential equation in C=C([-1,0],ℝ3) as
(70)y˙(t)=Lρ(yt)+f(ρ,yt),
where C=C([-1,0],ℝ3) is the Banach space of continuous functions mapping the interval [-τ,0] into ℝ3, y(t)=(y1(t),y2(t),y3(t))T∈ℝ3, yt(θ)=y(t+θ) for θ∈[-τ,0] and Lρ:C→ℝ3, f:ℝ×C→ℝ3 are defined as follows, respectively:
(71)Lρ(ϕ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×(ϕ1(0)ϕ2(0)ϕ3(0))+(τj+ρ)(0-βS~*00βS~*0000)×(ϕ1(-1)ϕ2(-1)ϕ3(-1)),(72)f(ρ,ϕ)=(τj+ρ)(-rkϕ12(0)-βϕ1(0)ϕ2(-1)βϕ1(0)ϕ2(-1)∑j=1∞fj(3)ϕ3j(0)),
where fj(3)=(di(-v*R/(wR-c))/dRi)|(S~*,I~*,R~*), j=1,2,3.It is easy to show thatLρ is a continuous linear function mapping C into ℝ3. According to Riesz representation theorem [40], there exists a 3×3 matrix function η(θ,ρ) of bounded variation for θ∈[-1,0] such that
(73)Lρ(ϕ)=∫-10‍dη(θ,ρ)ϕ(θ),
where ϕ∈C([-1,0],ℝ3).In fact, we can choose(74)η(θ,ρ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×δ(θ)-(τj+ρ)(0-βS~*00βS~*0000)δ(θ+1),
where δ denotes the Dirac delta function.Ifϕ is any given function in C([-1,0],ℝ3) and y(ϕ) is the unique solution of the linearized equation y˙(t)=Lρ(yt) of (70) with initial function ϕ at zero, then the solution operator T~(t):C→C is defined by
(75)T~(t)ϕ=yt(ϕ),t≥0.It follows from Lemma  7.1.1 in [39] that T~(t),  t≥0 is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator Aρ of T~(t), t≥0 is as follows:
(76)Aρ(ϕ)={dϕ(θ)dθ,θ∈[-1,0)∫-10‍dη(ρ,s)ϕ(s),θ=0,
for ϕ∈C1([-1,0],ℝ3), the space of functions mapping the interval [-1,0] into ℝ3 which have continuous first derivative and also define
(77)R(ρ)(ϕ)={0,θ∈[-1,0)f(ρ,ϕ),θ=0
then model system (70) is equivalent to
(78)y˙t=A(ρ)yt+R(ρ)yt.Forψ∈C1([0,1],(ℝ3)*), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative, define
(79)A*ϕ(s)={-dψ(s)ds,s∈(0,1]∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(80)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). It follows from the above analysis A(0) and A* are adjoint operators.By virtue of discussion in Section4.2.1, ±iω0τj are eigenvalues of A(0). Hence, they are also eigenvalues of A*. In the following, eigenvectors of A(0) and A* are corresponding to iω0τj and -iω0τj, respectively.Supposeq(θ)=(1,a,b)Teiω0τjθ is the eigenvectors of A(0) corresponding to iω0τj, which derives that A(0)q(θ)=iω0τjq(θ). By using the definition of A(0), (71) and (72), it gives that
(81)(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)q(0)+(0-βS~*00βS~*0000)q(-1)=iω0q(0).Forq(-1)=q(0)e-iω0τj, then it can be obtained that
(82)a=-r(μ1+m)+ikβω0kβ(μ1+m)e-iω0τj,b=m(wR~*-c)2[r(m+μ1)+ikβω0]kβ(μ1+m)[cv*-(μ2+iω0)(wR~*-c)2]e-iω0τj.Similarly, it follows from simple computation that eigenvectorq*(s)=J(1,a*,b*)eiω0τjs of A* is corresponding to -iω0τj, where
(83)a*=-r(μ1+m)-ikβω0kβ(μ1+m)eiω0τj,b*=m(wR~*-c)2[r(μ1+m)-ikβω0]kβ(μ1+m)[cv*-(μ2-iω0)(wR~*-c)2]eiω0τj.In order to assume〈q*(s),q(θ)〉=1, we need to determine the value of J in the following part.By virtue of (80), it derives that
(84)〈q*(s),q(θ)〉=J¯(1,a¯*,b¯*)(1,a,b)T-∫-10‍∫ξ=0θ‍J¯(1,a*,b*)e-iω0τj(ξ-θ)dη(θ)mmmmmmm×(1,a,b)Teiω0τjξdξ=J¯[∫-10‍1+aa¯*+bb¯*mn-∫-10‍(1,a¯*,b¯*)θeiω0τjθdη(θ)(1,a,b)T]=J¯[1+aa¯*+bb¯*+kb¯*τj(R~*+bS~*)eiω0τj].Hence, we can chooseJ as follows:
(85)J=11+a¯a*+b¯b*+kb*τj(R~*+bS~*)eiω0τj.Next, we will compute the coordinate to describe the centre manifoldC0 at ρ=0. Let yt be the solution of (78) when ρ=0.Define(86)z(t)=〈q*,yt〉,W(t,θ)=yt(θ)-2Re{z(t)q(θ)}.On the center manifoldC0, it derives that
(87)W(t,θ)=W(z(t),z¯(t),θ),
where
(88)W(z(t),z¯(t),θ)=W20(θ)z22+W11(θ)zz¯+W02z¯22+⋯,z and z¯ are local coordinates for center manifold C0 in the direction of q* and q¯*.It is noted thatW is real if yt is real, and we only consider real solutions. For solution yt∈C0 of (78), since ρ=0, it derives that
(89)z˙(t)=iω0τjz+q¯*(0)f(0,W(z,z¯,0)+2Re{zq(θ)})≜iω0τjz+q¯*(0)f0(z,z¯).The above equation can be rewritten as follows:(90)z˙(t)=iω0τjz(t)+g(z,z¯),
where
(91)g(z,z¯)=q¯*(0)f0(z,z¯)=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.It follows from (86) and (88) that
(92)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W20(θ)z22+W11(θ)zz¯+W02(θ)z¯22+(1,a,b)Teiω0τjθz+(1,a¯,b¯)Te-iω0τjθz¯+⋯.By virtue of (72), (91), and (92), it derives that(93)g(z,z¯)=q¯*(0)f0(z,z¯)=q¯*(0)f(0,yt)=τjJ¯(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0))=τjJ¯[∑j=1∞‍-rky1t2(0)-βy1t(0)y2t(-1)mml+a¯*βy1t(0)y2t(-1)+∑j=1∞‍b¯*fj(3)y3tj(0)]=τjJ¯[∑j=1∞‍-rky1t2(0)+(a¯*-1)βy1t(0)y2t(-1)mml+∑j=1∞‍b¯*fj(3)y3tj(0)]=-rkτjJ¯[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mmmmll+W02(1)(0)z¯22+o(|(z,z¯)|3)]2+(a¯*-1)βτjJ¯×[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mm+W02(1)(0)z¯22+o(|z,z¯|3)]×[e-iω0τjz+eiω0τjz¯+W20(2)(-1)z22mm+W11(2)(-1)zz¯+W02(2)(-1)z¯22+o(|z,z¯|3)]+τjJ¯b¯*[f1(3)(z+z¯+W20(3)(0)z22+W11(3)(0)zz¯mmmmmmmml+W02(1)(0)z¯22+o(|z,z¯|3))mmmmn+f2(3)(z+z¯+W20(3)(0)z22mmmmmmmmmm+W11(3)(0)zz¯+W02(1)(0)z¯22mmmmmmmmmm+W02(1)(0)z¯22+o(|z,z¯|3))+⋯].By comparing the coefficients with (91), it gives that
(94)g20=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W20(3)(0)2+f2(3))],g11=τjJ¯[(f1(3)W20(3)(0)2+f2(3))-2rk+2(a¯*-1)βcos⁡ω0τjmml+f1(3)W11(3)(0)+2f2(3)(f1(3)W20(3)(0)2+f2(3))],g02=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W02(3)(0)2+f2(3))],g21=2τjJ¯[W20(3)(-1)+W20(1)(0)eiω0τj2β(a¯*-1)mmll×(W11(2)(-1)+W20(3)(-1)+W20(1)(0)eiω0τj2mmmmml+W11(1)(0)e-iω0τjW20(3)(-1)+W20(1)(0)eiω0τj2)mmll-rk(W20(1)(0)+2W11(1)(0))mmll+f2(3)(W20(3)(0)+2W11(3)(0))W20(3)(-1)+W20(1)(0)eiω0τj2].Sinceg21 is associated with W20(θ) and W11(θ), further attempts should be carried out to compute W20(θ) and W11(θ).By virtue of (78) and (86), we have
(95)W˙=y˙t-z˙q-z¯˙q¯={AW-2Re{q¯*(0)f0q(θ)},θ∈[-1,0)AW-2Re{q¯*(0)f0q(0)+f0},θ=0.≜AW+H(z,z¯,θ),
where
(96)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯.By substituting the corresponding series into (95) and comparing the coefficients, we have
(97)(A-2iω0τj)W20(θ)=-H20(θ),AW11(θ)=-H11(θ),….It follows from (95) that for θ∈[-1,0)(98)H(z,z¯,θ)=-q¯*(0)f0q(θ)-q*(0)f¯0q¯(θ)=-g(z,z¯)q(θ)-g¯(z,z¯)q¯(θ).By comparing coefficients in (96) with those in (94), it derives that
(99)H20(θ)=-g20q(θ)-g¯02q¯(θ),(100)H11(θ)=-g11q(θ)-g¯11q¯(θ).Based on the definition ofA and (97) and (99), it can be obtained that
(101)W˙20(θ)=2iω0τjW20(θ)+g20q(θ)+g¯02q¯(θ).Forq(θ)=(1,a,b)Teiω0τjθ,
(102)W20(θ)=ig20ω0τjq(0)eiω0τjθ+ig¯023ω0τjq¯(0)e-iω0τjθ+G1e2iω0τjθ,
where G1=(G1(1),G1(2),G1(3)) is a constant vector.Similarly, it follows from (97) and (100) that
(103)W11(θ)=-ig11ω0τjq(0)eiω0τjθ+ig¯11ω0τjq¯(0)e-iω0τjθ+G2,
where G2=(G2(1),G2(2),G2(3)) is a constant vector.Subsequently, values ofG1 and G2 should be computed. By using the definition of A and (95), we have
(104)∫-10‍dηw20(θ)=2iω0τjW20(0)-H20(0),(105)∫-10‍dη(θ)W11(θ)=-H11(0),
where η(θ)=η(0,θ). Based on (95), it derives that in the case of θ=0,
(106)H(z,z¯,0)=-2Re{q¯*(0)f0q(0)}+f(0)=-q¯*(0)f0q(0)-q*(0)f¯0q¯(0)+f0=-g(z,z¯)q(0)-g¯(z,z¯)q¯(0)+f0,
which follows that
(107)H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯=-q(0)(g20z22+g11zz¯+g02z¯22+⋯)-q¯(0)(g¯20z¯22+g¯11zz¯+g¯02z22+⋯)+f0.By virtue of (72), it gives that
(108)f0=τk(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0)).By virtue of (86), it can be obtained that
(109)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W(t,θ)+z(t)q(θ)+z¯(t)q¯(θ)=W20(θ)z22+W21(θ)zz¯+z(t)q(θ)+z¯(t)q¯(θ)+⋯.
Then we have
(110)f0=τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3))z2+τk(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0))zz¯+⋯.According to (107) and (110), we have
(111)H20(0)=-g20q(0)-g¯02q¯(0)+2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),(112)H11(0)=-g11q(0)-g¯11q¯(0)+τj(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Sinceiω0τj is the eigenvalue of A(0) and q(0) is the corresponding eigenvector, we obtain that
(113)(iω0τjI-∫-10‍eiω0τjθdη(θ))q(0)=0,(-iω0τjI-∫-10‍e-iω0τjθdη(θ))q¯(0)=0,
where I is identity matrix.By substituting (102) and (111) into (104), it can be obtained that
(114)(2iω0τjI-∫-10‍e2iω0τjθdη(θ))G1=2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),
which can be rewritten as follows:
(115)(2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2)G1=2(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)).Based on Grammar Law [43], G1(1), G1(2), and G1(3) can be obtained as follows:(116)G1(1)=2U1|-rk-βe-iω0τjβS~*e-2iω0τj0βe-iω0τj2iω0+(μ1+m)-βS~*e-2iω0τj0f1(3)W20(3)(0)2+f2(3)-m2iω0+μ2-cv*(wR~*-c)2|,G1(2)=2U1|2iω0+rS~*k-rk-βe-iω0τj0-βI~*βe-iω0τj00f1(3)W20(3)(0)2+f2(3)2iω0+μ2-cv*(wR~*-c)2|,G1(3)=2U1|2iω0+rS~*kβS~*e-2iω0τj-rk-βe-iω0τj-βI~*2iω0+(μ1+m)-βS~*e-2iω0τjβe-iω0τj0-mf1(3)W20(3)(0)2+f2(3)|,
where
(117)U1=|2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2|.Similarly, substituting (103) and (112) into (105), it can be obtained that
(118)(rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2)G2=2(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Based on Grammar Law [43], G2(1), G2(2), and G2(3) can be obtained as follows:
(119)G2(1)=1U2|-2rk-2βcos⁡ω0τjβS~*02βcos⁡ω0τjμ1+m-βS~*02f2(3)+f1(3)W11(3)(0)-mμ2-cv*(wR~*-c)2|,G2(2)=1U2|rS~*k-2rk-2βcos⁡ω0τj0-βI~*2βcos⁡ω0τj002f2(3)+f1(3)W11(3)(0)μ2-cv*(wR~*-c)2|,G2(3)=1U2|rS~*kβS~*-2rk-2βcos⁡ω0τj-βI~*μ1+m-βS~*2βcos⁡ω0τj0-m2f2(3)+f1(3)W11(3)(0)|,
where
(120)U2=|rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2|.It follows from the above computation and analysis thatW20(θ) and W11(θ) can be determined based on (102) and (103).Furthermore, we can computeg21 based on (94). Hence, the following values can be computed as follows:
(121)d1(0)=i2ω0τj(g20g11-2|g11|2-|g02|23)+g212,δ2=-Re{d1(0)}Re{λ′(τj)},γ2=2Re(d1(0)),T2=Im⁡{d1(0)}+δ2Im⁡{λ′(τj)}ω0τj.By using similar arguments in [49], some properties of bifurcating periodic solutions of model (4) in the center manifold at the critical values are discussed in this paper. Based on the analysis in Section 4.2.2 of this paper, the following theorem can be concluded.Theorem 19.
The properties of Hopf bifurcation are determined by values in (121).(i)
δ2 determines directions of Hopf bifurcation: if δ2>0 (δ2<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τj (τ<τj).(ii)
γ2 determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if γ2<0 (γ2>0).(iii)
T2 determines the period of bifurcating periodic solutions: period increases (decreases) if T2>0 (T2<0).


## 4.2.1. Local Stability and Hopf Bifurcation

As analyzed in the above subsection, in the case of time delay and positive economic interest of harvesting0<v*<v~ where v~ is defined in (55), there are two interior equilibria P~1* and P~2* for model system (4) with respect to the positive economic interest v*.In this subsection, we only investigate dynamical behavior of model system (4) around the interior equilibrium P~1*. Some symmetric results about the interior equilibrium P~2* can be also obtained, and P~1* is denoted as P~* for simplicity. According to Jacobian evaluated at the interior equilibrium P~* and the leading matrix Ξ(t) in model system (4), we can obtain the characteristic equation of model system (4) around P~*, which can be expressed as follows:
(57)|λ+rS~*kβS~*e-λτ00-βI~*λ-βS~*e-λτ+μ1+m000-mλ+μ1+E~*R~*00-wE~*-v*E~*|=0.
⟹(58)M(λ)+N(λ)e-λτ=0,
where
(59)M(λ)=λ3+m1λ2+m2λ+m3,N(λ)=n1λ2+n2λ+n3,m1=μ2+(μ1+m)(1+rkβ)-cv*(wR~*-c)2,m2=μ1+mkβ[r(μ1+m+μ2-cv*(wR~*-c)2)mmmmm+kβ(μ2-cv*(wR~*-c)2)],m3=r(μ1+m)2kβ[μ2-cv*(wR~*-c)2],n1=-μ1-m,n2=-(μ1+m)×[μ2+2r(μ1+m)-krβkβ-cv*(wR~*-c)2],n3=-r(μ1+m)[2(μ1+m)kβ-1]×[μ2-cv*(wR~*-c)2].Now substitutingλ=iσ (σ is a positive real number) into (58) and separating the real and imaginary parts, two transcendental equations can be obtained as follows:
(60)σ3-m2σ=n2σcos⁡(στ)-(n3-n1σ2)sin(στ),(61)m1σ2-m3=(n3-n1σ2)cos⁡(στ)+n2σsin(στ).By squaring and adding (60) and (61), it can be calculated that
(62)(n3-n1σ2)2+n22σ2=(m1σ2-m3)2+(σ3-m2σ)2,
⟹(63)σ6+C1σ4+C2σ2+C3=0,
where C1=m12-2m2-n12, C2=m22-2m1m3+2n1n3-n22, C3=m32-n32.According to the values ofCj,(j=1,2,3) and the Routh-Hurwitz criteria [43], a simple assumption that (58) has at least one positive real root σ0 is C3<0, which derives that kβ>3(μ1+m). Hence, under this assumption, (58) will have a pair of purely imaginary roots of the form ±iσ0.By eliminatingsin(στ) from (60) and (61), it can be calculated that the τj* corresponding to σ0 is as follows:
(64)τj*=1σ0×arccos[n2σ02(σ02-m2)+(n3-n1σ02)(m1σ02-m3)(n3-n1σ02)2+(n2σ0)2]+2jπσ0,
where j=0,1,2,….By using Butler’s lemma [48], model system (4) is locally stable around P~* for τ<τ0*. Subsequently, conditions for existence of Hopf bifurcation in [39] are utilized to investigate whether Hopf bifurcation occurs as τ increases through τj*.Theorem 17.
Ifkβ>3(μ1+m), then model system (4) undergoes Hopf bifurcation around the interior equilibrium P~* when τ=τj*,  j=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Proof.
As mentioned above, letλ=iσ0 represent the purely imaginary root of (58). It follows from (58) that |M(iσ0)|=|N(iσ0)|, which determines a set of possible values of σ0.
In the following part, we determine the direction of motion ofλ=iσ0 as τ is varied; namely, we determine
(65)Θ=sign⁡[d(Reλ)dτ]λ=iσ0=sign⁡[Re(dλdτ)-1]λ=iσ0.
By differentiating (58) with respect to τ, it can be obtained that
(66)(dλdτ)-1=3λ2+2m1λ+m2λe-λτ(n1λ2+n2λ+n3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=3λ2+2m1λ+m2-λ(λ3+m1λ2+m2λ+m3)+2n1λ+n2λ(n1λ2+n2λ+n3)-τλ=2λ3+m1λ2-m3-λ2(λ3+m1λ2+m2λ+m3)+n1λ2-n3λ2(n1λ2+n2λ+n3)-τλ.
From (62) and the above equation, it can be obtained that
(67)Θ=sign⁡[Re(dλdτ)-1]λ=iσ0=1σ02sign⁡[(m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)(m1σ02-m3)2+(σ03-m2σ0)2mmmml+(n1σ02+n3)(n3-n1σ02)(n3-n1σ02)2+(n2σ0)2]=sign⁡[((m1σ02-m3)2+(σ03-m2σ0)2)-1((m3+m1σ02)(m1σ02-m3)+2σ04(σ02-m2)mmmm+(n1σ02+n3)(n3-n1σ02))mmm×((m1σ02-m3)2+(σ03-m2σ0)2)-1]=sign⁡[2σ06+(m12-2m2-n12)σ04+n32-m32(m1σ02-m3)2+(σ03-m2σ0)2],=sign⁡[2σ06+C1σ04-C3(m1σ02-m3)2+(σ03-m2σ0)2].
According to the values ofCj,(j=1,2,3) given in (58) of this paper, it is easy to show that C1=m12-2m2-n12=r2(μ1+m)2/k2β2+[μ2-cv*/(wR~-c)2]2>0.
Furthermore, ifkβ>3(μ1+m), then it can be shown that C3>0. Hence, it can be concluded that 2σ06+C1σ04-C3>0, which derives sign⁡[d(Reλ)/dτ]τ=τj*,σ=σ0>0. Consequently, the transversality condition holds and Hopf bifurcation occurs at σ=σ0, τ=τj*. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium P~* when τ>τ0* and ∥τ-τ0*∥≪1.Remark 18.
sign
⁡
[
d
(
Re
λ
)
/
d
τ
]
τ
=
τ
j
*
>
0 signifies that there exists at least one eigenvalue with positive real part for τ=τj*, and the conditions for Hopf bifurcation in [39] are also satisfied yielding the required periodic solution.


## 4.2.2. Properties of Hopf Bifurcation

By using normal theory and center manifold theorem [49], directions of Hopf bifurcation and stability of the bifurcating periodic solutions are discussed in this section. As analyzed in Section 4.1.2, when economic interest of harvesting 0<v*<v~ (v¯ is defined in (55)), it follows from implicit function theorem [41] and the fourth equation of model system (4) that E(t)=v*/(wR(t)-c). Furthermore, model system (4) can be transformed into the following form:
(68)S˙(t)=r(1-S(t)k)-βS(t)I(t-τ),I˙(t)=βS(t)I(t-τ)-μ1I(t)-mI(t),R˙(t)=mI(t)-μ2R(t)-v*R(t)wR(t)-c.Firstly, some transformations associated with component(S~*,I~*,R~*) of interior equilibrium P~* are given as follows:
(69)y1=S-S~*,y2=I-I~*,y3=R-R~*,y¯i(t)=yi(τt),τ=ρ+τj,ρ∈ℝ=(-∞,+∞).
Then ρ=0 is the Hopf bifurcation value of model system (4). Bars of variables are dropped for simplicity of notations; model system (4) is transformed to a functional differential equation in C=C([-1,0],ℝ3) as
(70)y˙(t)=Lρ(yt)+f(ρ,yt),
where C=C([-1,0],ℝ3) is the Banach space of continuous functions mapping the interval [-τ,0] into ℝ3, y(t)=(y1(t),y2(t),y3(t))T∈ℝ3, yt(θ)=y(t+θ) for θ∈[-τ,0] and Lρ:C→ℝ3, f:ℝ×C→ℝ3 are defined as follows, respectively:
(71)Lρ(ϕ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×(ϕ1(0)ϕ2(0)ϕ3(0))+(τj+ρ)(0-βS~*00βS~*0000)×(ϕ1(-1)ϕ2(-1)ϕ3(-1)),(72)f(ρ,ϕ)=(τj+ρ)(-rkϕ12(0)-βϕ1(0)ϕ2(-1)βϕ1(0)ϕ2(-1)∑j=1∞fj(3)ϕ3j(0)),
where fj(3)=(di(-v*R/(wR-c))/dRi)|(S~*,I~*,R~*), j=1,2,3.It is easy to show thatLρ is a continuous linear function mapping C into ℝ3. According to Riesz representation theorem [40], there exists a 3×3 matrix function η(θ,ρ) of bounded variation for θ∈[-1,0] such that
(73)Lρ(ϕ)=∫-10‍dη(θ,ρ)ϕ(θ),
where ϕ∈C([-1,0],ℝ3).In fact, we can choose(74)η(θ,ρ)=(τj+ρ)×(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)×δ(θ)-(τj+ρ)(0-βS~*00βS~*0000)δ(θ+1),
where δ denotes the Dirac delta function.Ifϕ is any given function in C([-1,0],ℝ3) and y(ϕ) is the unique solution of the linearized equation y˙(t)=Lρ(yt) of (70) with initial function ϕ at zero, then the solution operator T~(t):C→C is defined by
(75)T~(t)ϕ=yt(ϕ),t≥0.It follows from Lemma  7.1.1 in [39] that T~(t),  t≥0 is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator Aρ of T~(t), t≥0 is as follows:
(76)Aρ(ϕ)={dϕ(θ)dθ,θ∈[-1,0)∫-10‍dη(ρ,s)ϕ(s),θ=0,
for ϕ∈C1([-1,0],ℝ3), the space of functions mapping the interval [-1,0] into ℝ3 which have continuous first derivative and also define
(77)R(ρ)(ϕ)={0,θ∈[-1,0)f(ρ,ϕ),θ=0
then model system (70) is equivalent to
(78)y˙t=A(ρ)yt+R(ρ)yt.Forψ∈C1([0,1],(ℝ3)*), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative, define
(79)A*ϕ(s)={-dψ(s)ds,s∈(0,1]∫-10‍dηT(t,0)ψ(-t),s=0,
and a bilinear inner product
(80)〈ψ(s),ϕ(θ)〉=ψ¯(0)ϕ(0)-∫-10‍∫ξ=0θ‍ψ¯(ξ-θ)dη(θ)ϕ(ξ)dξ,
where η(θ)=η(θ,0). It follows from the above analysis A(0) and A* are adjoint operators.By virtue of discussion in Section4.2.1, ±iω0τj are eigenvalues of A(0). Hence, they are also eigenvalues of A*. In the following, eigenvectors of A(0) and A* are corresponding to iω0τj and -iω0τj, respectively.Supposeq(θ)=(1,a,b)Teiω0τjθ is the eigenvectors of A(0) corresponding to iω0τj, which derives that A(0)q(θ)=iω0τjq(θ). By using the definition of A(0), (71) and (72), it gives that
(81)(-rS~*k00βI~*-(μ1+m)00m-μ2+cv*(wR~*-c)2)q(0)+(0-βS~*00βS~*0000)q(-1)=iω0q(0).Forq(-1)=q(0)e-iω0τj, then it can be obtained that
(82)a=-r(μ1+m)+ikβω0kβ(μ1+m)e-iω0τj,b=m(wR~*-c)2[r(m+μ1)+ikβω0]kβ(μ1+m)[cv*-(μ2+iω0)(wR~*-c)2]e-iω0τj.Similarly, it follows from simple computation that eigenvectorq*(s)=J(1,a*,b*)eiω0τjs of A* is corresponding to -iω0τj, where
(83)a*=-r(μ1+m)-ikβω0kβ(μ1+m)eiω0τj,b*=m(wR~*-c)2[r(μ1+m)-ikβω0]kβ(μ1+m)[cv*-(μ2-iω0)(wR~*-c)2]eiω0τj.In order to assume〈q*(s),q(θ)〉=1, we need to determine the value of J in the following part.By virtue of (80), it derives that
(84)〈q*(s),q(θ)〉=J¯(1,a¯*,b¯*)(1,a,b)T-∫-10‍∫ξ=0θ‍J¯(1,a*,b*)e-iω0τj(ξ-θ)dη(θ)mmmmmmm×(1,a,b)Teiω0τjξdξ=J¯[∫-10‍1+aa¯*+bb¯*mn-∫-10‍(1,a¯*,b¯*)θeiω0τjθdη(θ)(1,a,b)T]=J¯[1+aa¯*+bb¯*+kb¯*τj(R~*+bS~*)eiω0τj].Hence, we can chooseJ as follows:
(85)J=11+a¯a*+b¯b*+kb*τj(R~*+bS~*)eiω0τj.Next, we will compute the coordinate to describe the centre manifoldC0 at ρ=0. Let yt be the solution of (78) when ρ=0.Define(86)z(t)=〈q*,yt〉,W(t,θ)=yt(θ)-2Re{z(t)q(θ)}.On the center manifoldC0, it derives that
(87)W(t,θ)=W(z(t),z¯(t),θ),
where
(88)W(z(t),z¯(t),θ)=W20(θ)z22+W11(θ)zz¯+W02z¯22+⋯,z and z¯ are local coordinates for center manifold C0 in the direction of q* and q¯*.It is noted thatW is real if yt is real, and we only consider real solutions. For solution yt∈C0 of (78), since ρ=0, it derives that
(89)z˙(t)=iω0τjz+q¯*(0)f(0,W(z,z¯,0)+2Re{zq(θ)})≜iω0τjz+q¯*(0)f0(z,z¯).The above equation can be rewritten as follows:(90)z˙(t)=iω0τjz(t)+g(z,z¯),
where
(91)g(z,z¯)=q¯*(0)f0(z,z¯)=g20z22+g11zz¯+g02z¯22+g21z2z¯2+⋯.It follows from (86) and (88) that
(92)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W20(θ)z22+W11(θ)zz¯+W02(θ)z¯22+(1,a,b)Teiω0τjθz+(1,a¯,b¯)Te-iω0τjθz¯+⋯.By virtue of (72), (91), and (92), it derives that(93)g(z,z¯)=q¯*(0)f0(z,z¯)=q¯*(0)f(0,yt)=τjJ¯(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0))=τjJ¯[∑j=1∞‍-rky1t2(0)-βy1t(0)y2t(-1)mml+a¯*βy1t(0)y2t(-1)+∑j=1∞‍b¯*fj(3)y3tj(0)]=τjJ¯[∑j=1∞‍-rky1t2(0)+(a¯*-1)βy1t(0)y2t(-1)mml+∑j=1∞‍b¯*fj(3)y3tj(0)]=-rkτjJ¯[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mmmmll+W02(1)(0)z¯22+o(|(z,z¯)|3)]2+(a¯*-1)βτjJ¯×[z+z¯+W20(1)(0)z22+W11(1)(0)zz¯mm+W02(1)(0)z¯22+o(|z,z¯|3)]×[e-iω0τjz+eiω0τjz¯+W20(2)(-1)z22mm+W11(2)(-1)zz¯+W02(2)(-1)z¯22+o(|z,z¯|3)]+τjJ¯b¯*[f1(3)(z+z¯+W20(3)(0)z22+W11(3)(0)zz¯mmmmmmmml+W02(1)(0)z¯22+o(|z,z¯|3))mmmmn+f2(3)(z+z¯+W20(3)(0)z22mmmmmmmmmm+W11(3)(0)zz¯+W02(1)(0)z¯22mmmmmmmmmm+W02(1)(0)z¯22+o(|z,z¯|3))+⋯].By comparing the coefficients with (91), it gives that
(94)g20=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W20(3)(0)2+f2(3))],g11=τjJ¯[(f1(3)W20(3)(0)2+f2(3))-2rk+2(a¯*-1)βcos⁡ω0τjmml+f1(3)W11(3)(0)+2f2(3)(f1(3)W20(3)(0)2+f2(3))],g02=2τjJ¯[(f1(3)W20(3)(0)2+f2(3))-rk+(a¯*-1)βe-iω0τjmmm+b¯*(f1(3)W02(3)(0)2+f2(3))],g21=2τjJ¯[W20(3)(-1)+W20(1)(0)eiω0τj2β(a¯*-1)mmll×(W11(2)(-1)+W20(3)(-1)+W20(1)(0)eiω0τj2mmmmml+W11(1)(0)e-iω0τjW20(3)(-1)+W20(1)(0)eiω0τj2)mmll-rk(W20(1)(0)+2W11(1)(0))mmll+f2(3)(W20(3)(0)+2W11(3)(0))W20(3)(-1)+W20(1)(0)eiω0τj2].Sinceg21 is associated with W20(θ) and W11(θ), further attempts should be carried out to compute W20(θ) and W11(θ).By virtue of (78) and (86), we have
(95)W˙=y˙t-z˙q-z¯˙q¯={AW-2Re{q¯*(0)f0q(θ)},θ∈[-1,0)AW-2Re{q¯*(0)f0q(0)+f0},θ=0.≜AW+H(z,z¯,θ),
where
(96)H(z,z¯,θ)=H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯.By substituting the corresponding series into (95) and comparing the coefficients, we have
(97)(A-2iω0τj)W20(θ)=-H20(θ),AW11(θ)=-H11(θ),….It follows from (95) that for θ∈[-1,0)(98)H(z,z¯,θ)=-q¯*(0)f0q(θ)-q*(0)f¯0q¯(θ)=-g(z,z¯)q(θ)-g¯(z,z¯)q¯(θ).By comparing coefficients in (96) with those in (94), it derives that
(99)H20(θ)=-g20q(θ)-g¯02q¯(θ),(100)H11(θ)=-g11q(θ)-g¯11q¯(θ).Based on the definition ofA and (97) and (99), it can be obtained that
(101)W˙20(θ)=2iω0τjW20(θ)+g20q(θ)+g¯02q¯(θ).Forq(θ)=(1,a,b)Teiω0τjθ,
(102)W20(θ)=ig20ω0τjq(0)eiω0τjθ+ig¯023ω0τjq¯(0)e-iω0τjθ+G1e2iω0τjθ,
where G1=(G1(1),G1(2),G1(3)) is a constant vector.Similarly, it follows from (97) and (100) that
(103)W11(θ)=-ig11ω0τjq(0)eiω0τjθ+ig¯11ω0τjq¯(0)e-iω0τjθ+G2,
where G2=(G2(1),G2(2),G2(3)) is a constant vector.Subsequently, values ofG1 and G2 should be computed. By using the definition of A and (95), we have
(104)∫-10‍dηw20(θ)=2iω0τjW20(0)-H20(0),(105)∫-10‍dη(θ)W11(θ)=-H11(0),
where η(θ)=η(0,θ). Based on (95), it derives that in the case of θ=0,
(106)H(z,z¯,0)=-2Re{q¯*(0)f0q(0)}+f(0)=-q¯*(0)f0q(0)-q*(0)f¯0q¯(0)+f0=-g(z,z¯)q(0)-g¯(z,z¯)q¯(0)+f0,
which follows that
(107)H20(θ)z22+H11(θ)zz¯+H02(θ)z¯22+⋯=-q(0)(g20z22+g11zz¯+g02z¯22+⋯)-q¯(0)(g¯20z¯22+g¯11zz¯+g¯02z22+⋯)+f0.By virtue of (72), it gives that
(108)f0=τk(-rky1t2(0)-βy1t(0)y2t(-1)βy1t(0)y2t(-1)∑j=1∞‍fj(3)y3tj(0)).By virtue of (86), it can be obtained that
(109)yt(θ)=W(t,θ)+2Re{z(t)q(θ)}=W(t,θ)+z(t)q(θ)+z¯(t)q¯(θ)=W20(θ)z22+W21(θ)zz¯+z(t)q(θ)+z¯(t)q¯(θ)+⋯.
Then we have
(110)f0=τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3))z2+τk(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0))zz¯+⋯.According to (107) and (110), we have
(111)H20(0)=-g20q(0)-g¯02q¯(0)+2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),(112)H11(0)=-g11q(0)-g¯11q¯(0)+τj(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Sinceiω0τj is the eigenvalue of A(0) and q(0) is the corresponding eigenvector, we obtain that
(113)(iω0τjI-∫-10‍eiω0τjθdη(θ))q(0)=0,(-iω0τjI-∫-10‍e-iω0τjθdη(θ))q¯(0)=0,
where I is identity matrix.By substituting (102) and (111) into (104), it can be obtained that
(114)(2iω0τjI-∫-10‍e2iω0τjθdη(θ))G1=2τj(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)),
which can be rewritten as follows:
(115)(2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2)G1=2(-rk-βe-iω0τjβe-iω0τjf1(3)W20(3)(0)2+f2(3)).Based on Grammar Law [43], G1(1), G1(2), and G1(3) can be obtained as follows:(116)G1(1)=2U1|-rk-βe-iω0τjβS~*e-2iω0τj0βe-iω0τj2iω0+(μ1+m)-βS~*e-2iω0τj0f1(3)W20(3)(0)2+f2(3)-m2iω0+μ2-cv*(wR~*-c)2|,G1(2)=2U1|2iω0+rS~*k-rk-βe-iω0τj0-βI~*βe-iω0τj00f1(3)W20(3)(0)2+f2(3)2iω0+μ2-cv*(wR~*-c)2|,G1(3)=2U1|2iω0+rS~*kβS~*e-2iω0τj-rk-βe-iω0τj-βI~*2iω0+(μ1+m)-βS~*e-2iω0τjβe-iω0τj0-mf1(3)W20(3)(0)2+f2(3)|,
where
(117)U1=|2iω0+rS~*kβS~*e-2iω0τj0-βI~*2iω0+(μ1+m)-βS~*e-2iω0τj00-m2iω0+μ2-cv*(wR~*-c)2|.Similarly, substituting (103) and (112) into (105), it can be obtained that
(118)(rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2)G2=2(-2rk-2βcos⁡ω0τj2βcos⁡ω0τj2f2(3)+f1(3)W11(3)(0)).Based on Grammar Law [43], G2(1), G2(2), and G2(3) can be obtained as follows:
(119)G2(1)=1U2|-2rk-2βcos⁡ω0τjβS~*02βcos⁡ω0τjμ1+m-βS~*02f2(3)+f1(3)W11(3)(0)-mμ2-cv*(wR~*-c)2|,G2(2)=1U2|rS~*k-2rk-2βcos⁡ω0τj0-βI~*2βcos⁡ω0τj002f2(3)+f1(3)W11(3)(0)μ2-cv*(wR~*-c)2|,G2(3)=1U2|rS~*kβS~*-2rk-2βcos⁡ω0τj-βI~*μ1+m-βS~*2βcos⁡ω0τj0-m2f2(3)+f1(3)W11(3)(0)|,
where
(120)U2=|rS~*kβS~*0-βI~*μ1+m-βS~*00-mμ2-cv*(wR~*-c)2|.It follows from the above computation and analysis thatW20(θ) and W11(θ) can be determined based on (102) and (103).Furthermore, we can computeg21 based on (94). Hence, the following values can be computed as follows:
(121)d1(0)=i2ω0τj(g20g11-2|g11|2-|g02|23)+g212,δ2=-Re{d1(0)}Re{λ′(τj)},γ2=2Re(d1(0)),T2=Im⁡{d1(0)}+δ2Im⁡{λ′(τj)}ω0τj.By using similar arguments in [49], some properties of bifurcating periodic solutions of model (4) in the center manifold at the critical values are discussed in this paper. Based on the analysis in Section 4.2.2 of this paper, the following theorem can be concluded.Theorem 19.
The properties of Hopf bifurcation are determined by values in (121).(i)
δ2 determines directions of Hopf bifurcation: if δ2>0 (δ2<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for τ>τj (τ<τj).(ii)
γ2 determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if γ2<0 (γ2>0).(iii)
T2 determines the period of bifurcating periodic solutions: period increases (decreases) if T2>0 (T2<0).


## 5. Numerical Simulation

In this section, some numerical simulations are provided to substantiate the theoretical results obtained in Section4 of this paper.


### 5.1. Numerical Simulation of State Feedback Controller for Singularity Induced Bifurcation and Local Stability Switch

In this subsection, values of parameters are partially taken from Section 5 of [11] and set in appropriate units, which can be found in Table 2. Numerical simulations are provided to illustrate the effectiveness of the state feedback controllers designed in Section 4.1 in the case of zero economic interest and positive economic interest, respectively.Table 2
Values of parameters for numerical simulation.
Parameter
Value
r
0.1
k
1
β
1
m
0.1
μ
1
0.1
μ
2
0.1
w
20
c
1In the case of zero economic interest, model (4) without time delay takes the following form:
(122)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1).By using Theorem11 of this paper, it can be shown that the model (4) without time delay has a singularity induced bifurcation around the interior equilibrium P*(0.2,0.08,0.05,0.06), and local stability switch occurs as v increases through 0.Based on the analysis in Section4.1.2 and Theorem 14 of this paper, a state feedback controller u(t)=l(E(t)-0.06) can be applied to stabilize model system (122) around P*, and then the model system (122) with the state feedback controller takes the following form:
(123)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)+l(E(t)-0.06).By using Theorem14 of this paper, if the feedback gain l satisfies l>2.2429, then model system (122) is stable around P* and singularity induced bifurcation of the model system (122) is also eliminated. The dynamical responses of model system (122) can be shown in Figure 1.Figure 1
Dynamical responses of model system (123) with state feedback controller, which shows that model system (123) is stable around (0.2,0.08,0.05,0.06).Furthermore, based on the analysis and inequality (55) in Section 4.1.2, there are two interior equilibria (denoted by P~1* and P~2*) when 0<v<0.00668. In the following part, we focus on the case of 0<v<0.00668, and the economic interest is set as v*=0.005∈(0,0.00668) in appropriate unit, which is arbitrarily selected within the interval (0,0.00668) and is enough to merit the theoretical analysis obtained in Section 4.1.2. By virtue of the given values of parameters in Table 2 and (54), two interior equilibria can be obtained as follows: P~1*(0.2,0.08,0.0557,0.0439) and P~2*(0.2,0.08,0.07175,0.0115). By using Theorem 11 of this paper, it can be shown that model system (124) is unstable around P~1* and P~2*:
(124)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005.Based on the analysis in Section4.1.2 of this paper, state feedback controllers u(t)=l~(E(t)-0.0439) and u(t)=l~(E(t)-0.0115) can be applied to stabilize the model system (124) around P~1* and P~2*, respectively. The model system (124) with respective state feedback controller takes the following form:
(125)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005+l~(E(t)-0.0439),(126)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005+l~(E(t)-0.0115).By using Theorem15 of this paper, if the feedback gain l~ satisfies l~>3.753, then model system (125) is stable around P~1* and model system (126) is stable around P~2*. The dynamical responses of model system (125) and (126) can be shown in Figures 2 and 3, respectively.Figure 2
Dynamical responses of model system (125) with state feedback controller, which shows that model system (125) is stable around (0.2,0.08,0.0557,0.0439).Figure 3
Dynamical responses of model system (126) with state feedback controller, which shows that model system (126) is stable around (0.2,0.08,0.07175,0.0115).


### 5.2. Numerical Simulation for Hopf Bifurcation and Local Stability Switch

In this subsection,values of parameters are partially taken from Section 5 of [11] and set in appropriate units, which can be found in Table 3. Numerical simulations are provided to support the theoretical findings obtained in Section 4.2 of this paper.Table 3
Values of parameters for numerical simulation.
Parameter
Value
r
0.1
k
1
β
1
m
0.13
μ
1
0.1
μ
2
0.1
w
20
c
1Based on the analysis and inequality (55) in Section 4.1.2, there are two interior equilibria (denoted by P~1* and P~2*) when 0<v<0.01831. In the following part, we focus on the case of 0<v<0.01831, and the economic interest is set as v*=0.012∈(0,0.01831) in appropriate unit, which is arbitrarily selected within the interval (0,0.01831) and is enough to merit the theoretical analysis obtained in Section 4.2. By virtue of the given values of parameters in Table 3 and (54), two interior equilibria can be obtained as follows: P~1*(0.2,0.08,0.0907,0.0147) and P~2*(0.2,0.08,0.0573,0.0821). Furthermore, it can be computed that kβ>3(μ1+m). Based on the analysis in Section 4.2.1, it satisfies the assumption that (58) has a positive root, and then the corresponding τ0*=2.7814 can be calculated by solving (64). It follows from (121) that δ2=1.2963>0, γ2=-0.3012<0, and T2=1.2467>0. Consequently, the interior equilibrium P~1* remains stable for τ<τ0*, and dynamical responses of model system (4) with τ=1.54 are plotted in Figure 4. It should be noted that τ=1.54 in Figure 4 is randomly selected in the interval (0,2.7814), which is enough to merit the above mathematical study. Only the dynamical responses and corresponding phase portrait of model (4) around P~1* are plotted; some symmetric results about P~2* can be also obtained. As τ increases through τ0*, a periodic solution caused by the phenomenon of Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from the interior equilibrium P~1*. Since δ2>0 and γ2<0, the Hopf bifurcation is supercritical, the directions of the Hopf bifurcation is τ>τ0*, and these bifurcating periodic solutions from the interior equilibrium P~1* at τ0* are stable. Dynamical responses of model (4) with τ=3>τ0* are plotted in Figure 5. Figures 6 and 7 show a limit cycle corresponding to the periodic solution in Figure 5 in theS-I-R andS-I-E space, respectively.Figure 4
Dynamical responses of model system (4) with time delay τ=1.54, which shows that model system (4) is stable around (0.2,0.08,0.0907,0.0147).Figure 5
Dynamical responses of model system (4) with time delay τ=3, which shows that model system (4) is unstable around (0.2,0.08,0.0907,0.0147).Figure 6
A limit cycle corresponding to the periodic solution in Figure5 in theS-I-R space.Figure 7
Limit cycle corresponding to the periodic solution in Figure5 in theS-I-E space.


## 5.1. Numerical Simulation of State Feedback Controller for Singularity Induced Bifurcation and Local Stability Switch

In this subsection, values of parameters are partially taken from Section 5 of [11] and set in appropriate units, which can be found in Table 2. Numerical simulations are provided to illustrate the effectiveness of the state feedback controllers designed in Section 4.1 in the case of zero economic interest and positive economic interest, respectively.Table 2
Values of parameters for numerical simulation.
Parameter
Value
r
0.1
k
1
β
1
m
0.1
μ
1
0.1
μ
2
0.1
w
20
c
1In the case of zero economic interest, model (4) without time delay takes the following form:
(122)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1).By using Theorem11 of this paper, it can be shown that the model (4) without time delay has a singularity induced bifurcation around the interior equilibrium P*(0.2,0.08,0.05,0.06), and local stability switch occurs as v increases through 0.Based on the analysis in Section4.1.2 and Theorem 14 of this paper, a state feedback controller u(t)=l(E(t)-0.06) can be applied to stabilize model system (122) around P*, and then the model system (122) with the state feedback controller takes the following form:
(123)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)+l(E(t)-0.06).By using Theorem14 of this paper, if the feedback gain l satisfies l>2.2429, then model system (122) is stable around P* and singularity induced bifurcation of the model system (122) is also eliminated. The dynamical responses of model system (122) can be shown in Figure 1.Figure 1
Dynamical responses of model system (123) with state feedback controller, which shows that model system (123) is stable around (0.2,0.08,0.05,0.06).Furthermore, based on the analysis and inequality (55) in Section 4.1.2, there are two interior equilibria (denoted by P~1* and P~2*) when 0<v<0.00668. In the following part, we focus on the case of 0<v<0.00668, and the economic interest is set as v*=0.005∈(0,0.00668) in appropriate unit, which is arbitrarily selected within the interval (0,0.00668) and is enough to merit the theoretical analysis obtained in Section 4.1.2. By virtue of the given values of parameters in Table 2 and (54), two interior equilibria can be obtained as follows: P~1*(0.2,0.08,0.0557,0.0439) and P~2*(0.2,0.08,0.07175,0.0115). By using Theorem 11 of this paper, it can be shown that model system (124) is unstable around P~1* and P~2*:
(124)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005.Based on the analysis in Section4.1.2 of this paper, state feedback controllers u(t)=l~(E(t)-0.0439) and u(t)=l~(E(t)-0.0115) can be applied to stabilize the model system (124) around P~1* and P~2*, respectively. The model system (124) with respective state feedback controller takes the following form:
(125)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005+l~(E(t)-0.0439),(126)S˙(t)=0.1(1-S(t))-S(t)I(t),I˙(t)=S(t)I(t)-0.1I(t)-0.1I(t),R˙(t)=0.1I(t)-0.1R(t)-E(t)R(t),0=E(t)(20R(t)-1)-0.005+l~(E(t)-0.0115).By using Theorem15 of this paper, if the feedback gain l~ satisfies l~>3.753, then model system (125) is stable around P~1* and model system (126) is stable around P~2*. The dynamical responses of model system (125) and (126) can be shown in Figures 2 and 3, respectively.Figure 2
Dynamical responses of model system (125) with state feedback controller, which shows that model system (125) is stable around (0.2,0.08,0.0557,0.0439).Figure 3
Dynamical responses of model system (126) with state feedback controller, which shows that model system (126) is stable around (0.2,0.08,0.07175,0.0115).


## 5.2. Numerical Simulation for Hopf Bifurcation and Local Stability Switch

In this subsection,values of parameters are partially taken from Section 5 of [11] and set in appropriate units, which can be found in Table 3. Numerical simulations are provided to support the theoretical findings obtained in Section 4.2 of this paper.Table 3
Values of parameters for numerical simulation.
Parameter
Value
r
0.1
k
1
β
1
m
0.13
μ
1
0.1
μ
2
0.1
w
20
c
1Based on the analysis and inequality (55) in Section 4.1.2, there are two interior equilibria (denoted by P~1* and P~2*) when 0<v<0.01831. In the following part, we focus on the case of 0<v<0.01831, and the economic interest is set as v*=0.012∈(0,0.01831) in appropriate unit, which is arbitrarily selected within the interval (0,0.01831) and is enough to merit the theoretical analysis obtained in Section 4.2. By virtue of the given values of parameters in Table 3 and (54), two interior equilibria can be obtained as follows: P~1*(0.2,0.08,0.0907,0.0147) and P~2*(0.2,0.08,0.0573,0.0821). Furthermore, it can be computed that kβ>3(μ1+m). Based on the analysis in Section 4.2.1, it satisfies the assumption that (58) has a positive root, and then the corresponding τ0*=2.7814 can be calculated by solving (64). It follows from (121) that δ2=1.2963>0, γ2=-0.3012<0, and T2=1.2467>0. Consequently, the interior equilibrium P~1* remains stable for τ<τ0*, and dynamical responses of model system (4) with τ=1.54 are plotted in Figure 4. It should be noted that τ=1.54 in Figure 4 is randomly selected in the interval (0,2.7814), which is enough to merit the above mathematical study. Only the dynamical responses and corresponding phase portrait of model (4) around P~1* are plotted; some symmetric results about P~2* can be also obtained. As τ increases through τ0*, a periodic solution caused by the phenomenon of Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from the interior equilibrium P~1*. Since δ2>0 and γ2<0, the Hopf bifurcation is supercritical, the directions of the Hopf bifurcation is τ>τ0*, and these bifurcating periodic solutions from the interior equilibrium P~1* at τ0* are stable. Dynamical responses of model (4) with τ=3>τ0* are plotted in Figure 5. Figures 6 and 7 show a limit cycle corresponding to the periodic solution in Figure 5 in theS-I-R andS-I-E space, respectively.Figure 4
Dynamical responses of model system (4) with time delay τ=1.54, which shows that model system (4) is stable around (0.2,0.08,0.0907,0.0147).Figure 5
Dynamical responses of model system (4) with time delay τ=3, which shows that model system (4) is unstable around (0.2,0.08,0.0907,0.0147).Figure 6
A limit cycle corresponding to the periodic solution in Figure5 in theS-I-R space.Figure 7
Limit cycle corresponding to the periodic solution in Figure5 in theS-I-E space.


## 6. Conclusion

It is well known that the recovered host individuals are naturally immune to vector disease [1], and its potential economic interest can be commercially exploited. Furthermore, harvest effort is usually influenced by variation of economic interest under market economy. Consequently, it is necessary to discuss the coexistence and interaction mechanism of population within harvested epidemiological ecosystem as well as dynamical effect of harvest effort due to variation of economic interest.By introducing commercial harvest effort into model proposed in [11], a delayed hybrid mathematical model is established to investigate the dynamical effect of commercial harvesting and incubation time delay on epidemiological economic system, which extends the work done in [11] from a bioeconomic perspective. Positivity and persistence of solutions of model system are discussed in Theorems 3 and 10, respectively. The economic interest of commercial harvesting should be restricted within certain interval that guarantees the existence of interior equilibrium, which can be found in Remark 16. Since the interior equilibrium biologically interprets that susceptible, infective, and recovered host individuals survive as well as harvest on recovered host individuals exist, the bifurcation phenomena around the interior equilibria can reveal the instability mechanism of model system, which are theoretically relevant to infectious disease control and sustainable yield on recovered host individuals. Consequently, we will mainly concentrate on dynamical behavior and local stability switch around interior equilibrium of model system (4) in this paper. As analyzed in Theorem 11 of this paper, a singularity induced bifurcation occurs which leads to local stability switch in the case of positive economic interest of harvesting. In the perspective of practical viewpoint, a direct damage done by the singularity induced bifurcation to the proposed model is impulse phenomenon, which may lead to outbreak of infectious disease and hamper prosperous harvesting on recovered host individual population resource in the harvested ecosystem, which can be found in Remark 13. With the purpose of maintaining the economic interest at an ideal level, state feedback controllers are designed to stabilize model system around the desirable interior equilibria in the case of zero economic interest and positive economic interest, respectively. The design of the state feedback controller can be found in Theorems 14 and 15 of this paper. The theoretical results and numerical simulations obtained in this paper suggest that incorporating harvest effort on recovered host individuals can not only prevent the stability switch of model system,but also drive model system to stable equilibrium, which will contribute to the persistence and sustainable yield of the harvested ecosystem.Further attempts are made to understand the dynamical effect of incubation time delay and economic interest on local stability of model system around interior equilibrium. Local stability analysis reveals that incubation delay is responsible for local stability switch of the proposed model, and a family of periodic solutions bifurcate from the interior equilibrium which occurs as incubation delay increases through a critical threshold, which can be found in Theorem17. The direction and stability of Hopf bifurcation are also discussed in Theorem 19 of this paper, which reveals that Hopf bifurcation is supercritical, the directions of Hopf bifurcation is τ>τ0*, and these bifurcating periodic solutions from the interior equilibrium are stable. The model proposed in [11] does not discuss the harvest effort on economic population. For the model proposed in [11], the threshold value of incubation delay where Hopf bifurcation occurs in [11] is τ^0=2.0842. However, the harvest effort on recovered host individuals is considered in this paper. As calculated in Section 5.2 of this paper, the Hopf bifurcation occurs at τ0*=2.7814 in the case of positive economic interest. It is obvious that τ0*>τ^0, which implies that the harvesting has a stabilizing impact on the dynamical behavior of population dynamics; cyclic behavior caused by incubation delay can be deferred by introduction of commercial harvesting effort.It should be noted that some hybrid dynamical models are proposed in [23–28], which are utilized to discuss the interaction mechanism of harvested ecosystem from an economic perspective in recent years. Complex dynamical behavior and stability analysis in prey-predator ecosystems with stage-structured population and gestation delay are considered. However, as far as knowledge goes, nobody has explicitly proposed a mathematical model to discuss the dynamic effect of commercial harvest on epidemiological system under the market economy environment. The main objective of this paper is to investigate the transmission mechanism of infectious disease and dynamical effect of commercial harvest on population dynamics, especially the complex dynamical behavior and stability switch due to variation of incubation and commercial harvest economic interest, which makes the work studied in this paper has some new and positive features.



---

*Source: 639405-2014-05-12.xml*