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

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.

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−𝑆(𝑡)𝑘(𝑡)𝐼(),𝐼̇(𝑡)(𝑡)𝐼()−𝛽𝑆𝑡−𝜏=𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)−𝑚𝐼=𝑚𝐼−𝜇₂𝑅(𝑡), (1) where 𝑆(𝑡), 𝐼(𝑡), and 𝑅(𝑡) represent the population density of susceptible, infective, and recovered host individuals at time 𝑡, respectively. It is assumed that the population growth of susceptible host individuals is governed by the logistic growth with a carrying capacity 𝑘>0 as well as intrinsic birth rate constant 𝑟>0. 𝛽>0 is the average number of constants per infective per unit time and 𝜏≥0 denotes the incubation time, and 𝜇₁>0 and 𝜇₂>0 stand for the death rate of infective and recovered host individuals, respectively. 𝑚>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: ()()().NetEconomicRevenueNER=TotalRevenueTR−TotalCostTC (2)

Associated with model (1), an algebraic equation, which considers the economic interest 𝑣 of the harvest effort on recovered host individuals in epidemiological system, that is, 𝑅(𝑡), is established as follows: 𝐸(𝑡)(𝑤𝑅(𝑡)−𝑐)=𝑣, (3) where 𝐸(𝑡) represents the harvest effort on recovered host individuals at time 𝑡. 𝑣 represents the economic interest of harvest effort on the recovered host individuals. 𝑤 and 𝑐 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: 𝑆̇(𝑡)=𝑟1−𝑆(𝑡)𝑘(𝑡)𝐼(),𝐼̇(𝑡)(𝑡)𝐼()−𝛽𝑆𝑡−𝜏=𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)−𝑚𝐼=𝑚𝐼−𝜇₂𝑅(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))−𝐸0=𝐸)(𝑤𝑅−𝑐−𝑣, (4) where 𝑆(𝑡), 𝐼(𝑡), 𝑅(𝑡), 𝐸(𝑡), and other parameters share the same interpretations mentioned in (1) and (3), and initial conditions 𝜓=(𝜓₁,𝜓₂,𝜓₃,𝜓₄) for model system (4) are defined in the Banach space: 𝜓∈𝐶⁺[−𝜏,0],ℝ⁴₊∣𝜓₁(𝜃)=𝑆(𝜃),𝜓₂(𝜃)(𝜃),𝜓=𝐼₃(𝜃)=𝑅(𝜃),𝜓₄(𝜃)=𝐸(𝜃), (5) where ℝ⁴₊={(𝑆,𝐼,𝑅,𝐸)∈ℝ⁴∶𝑆≥0,𝐼≥0,𝑅≥0,𝐸≥0}. It is also assumed that 𝜓_𝑖(0)>0 (𝑖=1,2,3,4) for a biological reason.

Model system (4) can be expressed in the following form: Ξ(𝑡)𝑋̇(𝑡)=𝐹(𝑋(𝑡)), (6) where 𝑋(𝑡)=(𝑆(𝑡),𝐼(𝑡),𝑅(𝑡),𝐸(𝑡))^𝑇,Ξ(𝑡)=,𝐹(𝑋(𝑡=𝐹1000010000100000))¹(𝑋(𝑡𝐹))₂(𝑋(𝑡𝐹))₃(𝑋(𝑡𝐹))₄(𝑋(𝑡=𝑟𝑆(𝑡)))1−𝑘(𝑡)𝐼()(𝑡)𝐼()−𝛽𝑆𝑡−𝜏𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)(𝑡)(𝑡)−𝑚𝐼𝑚𝐼−𝜇₂𝑅(𝑡)(𝑡)𝑅(𝑡)𝐸(𝑡(𝑡)).−𝐸)(𝑤𝑅−𝑐−𝑣 (7)

Remark 1. The algebraic equation in model system (6) contains no differentiated variables; hence, the leading matrix Ξ(𝑡) 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 𝐹_𝑖∶ℝ₊⁴⁺¹→ℝ⁴ is locally Lipschitz and satisfies the condition, 𝐹_𝑖(𝑋(𝑡))|_𝑋∈ℝ⁴>0, where 𝐹_𝑖(𝑋(𝑡)) (𝑖=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,𝐴₀) for some 𝐴₀>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 𝑣≥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 Ω₀∈intΩ such that every solution of model system (4) with initial conditions will eventually enter and remain in region Ω₀.

Definition 5 (see [39]). Consider a metric space 𝑄 with metric 𝑑. The distance 𝑑(𝑥,𝑦) of a point 𝑥∈𝑄 from a subset 𝑌 of 𝑄 is defined by 𝑑𝑥,𝑦=inf_𝑦∈𝑌𝑑𝑥,𝑦. (8)

It is further assumed that 𝑄 is the closure of an open set 𝑄⁰, and 𝑄⁰=𝜕𝑄⁰ is nonempty and is the boundary of 𝑄⁰. Consequently, 𝑄⁰∪𝑄₀=𝑄,𝑄⁰∩𝑄₀=∅. We will also suppose that 𝑇(𝑡) is a 𝑄⁰ semigroup on 𝑄 satisfying 𝑇(𝑡)∶𝑄⁰⟶𝑄⁰,𝑇(𝑡)∶𝑄₀⟶𝑄₀. (9)

Let 𝑇_𝜕(𝑡)=𝑇(𝑡)|_𝑄0 and 𝐴_𝜕 be the global attractor for 𝑇_𝜕(𝑡).

Lemma 6 (see [39]). Suppose that 𝑇(𝑡) satisfies (9) and the following conditions hold.(i)There is a 𝑡₀≥0 such that 𝑇(𝑡) is compact for 𝑡>𝑡₀;(ii)𝑇(𝑡) is point dissipative in 𝑄;(iii)𝐴_𝜕=⋃_{𝑥∈𝐴𝜕}𝜔(𝑥) is isolated and has an acyclic covering 𝑍.Then 𝑇(𝑡) is uniformly persistent if and only if for each 𝑍_𝑖∈𝑍, 𝑊^𝑠(𝑍_𝑖)∩𝑄⁰=∅ for 𝑖=1,2,…,𝑛.

Lemma 7 (see [40]). Consider the following equation: 𝑢̇(𝑡)=𝑎𝑢(𝑡−𝜏)−𝑏𝑢(𝑡), (10) where 𝑎,𝑏,𝜏>0 and 𝑢(𝑡)>0 for all −𝜏≤𝑡≤0; it derives the following:(i)If 𝑎<𝑏, then lim_𝑡→+∞𝑢(𝑡)=0,(ii)If 𝑎>𝑏, then lim_𝑡→+∞𝑢(𝑡)=+∞.

Lemma 8. For any solutions of model system (4), we have limsup_𝑡→+∞𝑆(𝑡)≤𝑘,limsup_𝑡→+∞(𝑆(𝑡)(𝑡≤𝑘+𝐼))𝑟+𝜇₁+𝑚²4𝑟𝜇₁.+𝑚 (11)

Proof. By using Theorem 3, it follows from the first equation of model system (4) that 𝑆̇(𝑡)≤𝑟1−𝑆(𝑡)𝑘𝑆(𝑡), (12) which derives that limsup_𝑡→+∞𝑆(𝑡)≤𝑘.
According to Theorem 3 and the first and second equation of model system (4), it gives that 𝑆̇(𝑡)+𝐼̇(𝑡)≤𝑟1−𝑆(𝑡)𝑘𝑆(𝑡)−𝜇₁+𝑚𝐼(𝑡), (13) which derives that limsup_𝑡→+∞(𝑆(𝑡)+𝐼(𝑡))≤𝑘(𝑟+𝜇₁+𝑚)²/4𝑟(𝜇₁+𝑚).

Lemma 9. If 𝜇₁+𝑚<1, then (𝑆(𝑡),𝐼(𝑡)) of solution of model system (4) with initial conditions satisfies liminf_𝑡→+∞𝑆(𝑡)≥𝑆_𝜂,liminf_𝑡→+∞𝐼(𝑡)≥𝐼_𝜂, (14) where 𝑆_𝜂>0 and 𝐼_𝜂>0 are independent of corresponding initial values of model system (4).

Proof. Firstly, let 𝐶⁺([−𝜏,0],ℝ²₊) denote space of continuous functions mapping [−𝜏,0] into ℝ²₊, where ℝ²₊={(𝑥,𝑦)∣𝑥≥0,𝑦≥0}: 𝑄₁=𝜙₁,𝜙₂∈𝐶⁺[−𝜏,0],ℝ²₊∣𝜙₁(𝜃)[,𝑄=0,𝜃∈−𝜏,0]₂=𝜙₁,𝜙₂∈𝐶⁺[−𝜏,0],ℝ²₊∣𝜙₁(𝜃)𝜙>0,₂(𝜃)=0,𝜃∈[−𝜏,0]. (15)
Denote 𝑄₀=𝑄₁∪𝑄₂,  𝑄=𝐶⁺([−𝜏,0],ℝ²₊) and 𝑄⁰=int𝐶⁺([−𝜏,0],ℝ²₊).
Next, all conditions in Lemma 6 will be checked. In order to facilitate the proof, we consider the following subsystem of model system (4): 𝑆̇(𝑡)=𝑟1−𝑆(𝑡)𝑘(𝑡)𝐼(),𝐼̇(𝑡)(𝑡)𝐼()−𝛽𝑆𝑡−𝜏=𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)−𝑚𝐼(𝑡), (16) where 𝑆(𝜃)≥0,  𝐼(𝜃)≥0 are continuous on 𝜃∈[−𝜏,0] and 𝑆(0)>0,  𝐼(0)>0.
By Definition 5 and model system (16), it is easy to see that 𝑄⁰ and 𝑄₀ are positively invariant, and conditions (i) and (ii) of Lemma 6 clearly hold.
Since model system (16) possesses two constant solutions in 𝑄₀: 𝑃₀∈𝑄₁, 𝑃₁∈𝑄₂ with the following form: 𝑃₀=𝜙₁,𝜙₂∈𝐶⁺[−𝜏,0],ℝ²₊∣𝜙₁(𝜃)=𝜙₂(𝜃)[,𝑃=0,𝜃∈−𝜏,0]₁=𝜙₁,𝜙₂∈𝐶⁺[−𝜏,0],ℝ²₊∣𝜙₁(𝜃)=1,𝜙₂(𝜃)[.=0,𝜃∈−𝜏,0] (17)
It follows from simple computation that 𝑆̇(𝑡)_{(𝜙1,𝜙2)∈𝑄1}=0,𝑆(𝑡)_{(𝜙1,𝜙2)∈𝑄1}=0for𝑡≥0. (18)
Moreover, it follows from the second equation of model system (16) that 𝐼̇(𝑡)_{(𝜙1,𝜙2)∈𝑄1}=−𝜇₁+𝑚𝐼(𝑡)≤0, (19) which reveals that all points in 𝑄₁ approach to 𝑃₀; that is, 𝑄₁=𝑊^𝑠(𝑃₀). By using the similar analysis mentioned above, it can be also concluded that all points in 𝑄₂ approach to 𝑃₁; that is, 𝑄₂=𝑊^𝑠(𝑃₁). Based on the above analysis, it shows that invariant sets 𝑃₀ and 𝑃₁ are isolated invariant, and {𝑃₀,𝑃₁} is isolated and an acyclic covering. It can be concluded that condition (iii) of Lemma 6 holds.
Finally, we will show that 𝑊^𝑠(𝑃_𝑖)∩𝑄⁰=∅ for 𝑖=0,1. Based on the definition of 𝑃₀, it is easy to show that 𝑊^𝑠(𝑃₀)∩𝑄⁰=∅. We will show 𝑊^𝑠(𝑃₁)∩𝑄⁰=∅ in the following part.
If 𝑊^𝑠(𝑃₁)∩𝑄⁰≠∅, then there exists a positive solution (𝑆(𝑡),𝐼(𝑡)) to model system (16) with lim_𝑡→+∞(𝑆(𝑡),𝐼(𝑡))=(1,0). If 𝜇₁+𝑚<1, then 𝜇₁+𝑚<1−𝜖 holds for sufficiently small 𝜖>0 and there exists a positive constant 𝑇=𝑇(𝜖) such that 𝑆(𝑡)>1−𝜖>0, and 0<𝐼(𝑡)<𝜖 for all 𝑡≥𝑇.
By the second equation of model system (16), it derives that 𝐼̇(𝑡)≥(1−𝜖)𝐼(𝑡−𝜏)−𝜇₁+𝑚𝐼(𝑡) (20) holds for all 𝑡≥𝑇+𝜏.
Consider the following equation: 𝑥̇(𝑡)=(1−𝜖)𝑥(𝑡−𝜏)−𝜇₁𝑥(𝑡)𝑥(𝑡)(𝑡)+𝑚,𝑡≥𝑇+𝜏,=𝐼,𝑇≤𝑡≤𝑇+𝜏. (21)
Based on (21) and the comparison principle, it derives that 𝐼(𝑡)≥𝑥(𝑡) for all 𝑡>𝑇.
On the other hand, if 𝜇₁+𝑚<1, then it follows from Lemma 7 that lim_𝑡→+∞𝑥(𝑡)=+∞ for all solutions of (21). It can be concluded that lim_𝑡→+∞𝐼(𝑡)=∞, which is a contradiction to 𝐼(𝑡)<𝜖. Consequently, it can be derived that 𝑊^𝑠(𝑃₁)∩𝑄⁰=∅.
According to the above analysis, all conditions of Lemma 6 hold. By using Lemma 6, it can be obtained that liminf_𝑡→+∞𝑆(𝑡)≥𝑆_𝜂,liminf_𝑡→+∞𝐼(𝑡)≥𝐼_𝜂, (22) where 𝑆_𝜂>0 and 𝐼_𝜂>0 are independent of the corresponding initial values of model system (4).

Theorem 10. If 𝜇₁+𝑚<1, 𝑐𝜇₂<𝑚𝑤𝐼_𝜂, and 0≤𝑣<𝑐𝜇₂+𝑚𝑤𝐼_𝜂, then all solutions of model system (4) with initial conditions are persistent.

Proof. According to Lemmas 8 and 9, it can be obtained that 𝑆_𝜂(𝑡)(𝑡)(𝑡)≤𝑘≤𝑆≤𝑘,𝑆+𝐼𝑟+𝜇₁+𝑚²4𝑟𝜇₁,+𝑚 (23) hold for all 𝑡>0, which derive that 𝐼_𝜂(𝑡)≤𝑘≤𝐼𝑟+𝜇₁+𝑚²4𝑟𝜇₁+𝑚−𝑆_𝜂. (24)
When the economic interest 𝑣=0, it follows from Theorem 3 and the fourth equation of model system (4) that 𝑅(𝑡)=𝑐𝑤 (25) and 𝑅̇(𝑡)=0. Based on the third equation of model system (4), it can be computed that 𝐸(𝑡)=(𝑚𝑤/𝑐)𝐼(𝑡)−𝜇₂. According to (24), it derives that 0<𝑚𝑤𝐼_𝜂−𝑐𝜇₂𝑐(𝑡)≤𝑘≤𝐸𝑚𝑤𝑟+𝜇₁+𝑚²−4𝑆_𝜂𝑟𝜇₁+𝑚4𝑐𝑟𝜇₁+𝑚−𝜇₂ (26) provided that 𝑐𝜇₂<𝑚𝑤𝐼_𝜂.
In the case of 𝑣>0, it derives that 𝐸(𝑡)=𝑣/(𝑤𝑅(𝑡)−𝑐) based on implicit function theory [41]. According to the third equation of model system (4), it can be obtained that 𝑅̇(𝑡)≥𝑚𝐼_𝜂+𝑐𝜇₂−𝑣𝑤−𝜇₂𝑅(𝑡), (27) which derives that liminf_𝑡→+∞𝑅(𝑡)≥𝑚𝑤𝐼_𝜂+𝑐𝜇₂−𝑣𝑤𝜇₂∶=𝑅>0, (28) provided that 0<𝑣<𝑐𝜇₂+𝑚𝑤𝐼_𝜂.
It follows from Theorem 3 and the third equation of model system (4) that 𝑅̇(𝑡)≤𝑚𝐼(𝑡)−𝜇₂𝑅(𝑡), (29) which derives that limsup_𝑡→+∞𝑅(𝑡)≤𝑚𝑘𝑟+𝜇₁+𝑚²−4𝑆_𝜂𝑟𝜇₁+𝑚4𝑟𝜇₂𝜇₁+𝑚∶=𝑅. (30)
Hence, it gives that 𝑣/(𝑤𝑅−𝑐)≤𝐸(𝑡)≤𝑣/(𝑤𝑅−𝑐), and it can be rewritten as follows: 4𝑟𝜇₂𝑣𝜇₁+𝑚𝑤𝑚𝑘𝑟+𝜇₁+𝑚²−4𝑆_𝜂𝑟𝜇₁+𝑚−4𝑟𝑐𝜇₂𝜇₁(𝑡)≤𝜇+𝑚≤𝐸²𝑣𝑚𝑤𝐼_𝜂.−𝑣 (31)
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 𝑣=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<𝑣<𝑐𝜇₂+𝑚𝑤𝐼_𝜂.

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 𝑣=0 is a bifurcation value. Furthermore, local stability switch occurs as 𝑣 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, 𝑣=0. An interior equilibrium can be obtained as follows: 𝑃^∗(𝑆^∗,𝐼^∗,𝑅^∗,𝐸^∗), where 𝑆^∗=(𝜇₁+𝑚)/𝛽, 𝐼^∗=𝑟(𝑘𝛽−𝜇₁−𝑚)/𝑘𝛽², 𝑅^∗=𝑐/𝑤, and 𝐸^∗=(𝑤𝑚𝑟(𝑘𝛽−𝜇₁−𝑚)−𝑐𝑘𝜇₂𝛽²)/𝑐𝑘𝛽².
According to biological interpretation of the interior equilibrium, it follows that 𝑆^∗>0,  𝐼^∗>0,𝑅^∗>0 and 𝐸^∗>0. In order to guarantee the existence of interior equilibrium, some inequalities are satisfied: 𝑘𝛽−𝜇₁−𝑚>0,𝑤𝑚𝑟𝑘𝛽−𝜇₁−𝑚−𝑐𝑘𝜇₂𝛽²>0. (32)
Let 𝑣 be a bifurcation parameter, 𝐻(𝑡)=(𝑆(𝑡),𝐼(𝑡),𝑅(𝑡))^𝑇, ℎ₁(𝐻(𝑡)(𝑡))=𝑟𝑆(𝑡),𝐸,𝑣1−𝑘(𝑡)𝐼(),(𝑡)𝐼()−𝛽𝑆𝑡−𝜏𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)(𝑡),(𝑡)−𝑚𝐼𝑚𝐼−𝜇₂𝑅(𝑡)(𝑡)𝑅(𝑡),ℎ−𝐸₂(𝐻(𝑡),𝐸(𝑡),𝑣)=𝐸(𝑡)(𝑤𝑅(𝑡)−𝑐)−𝑣. (33)
It can be calculated that trace𝐷_𝐸ℎ₁adj𝐷_𝐸ℎ₂𝐷_𝐻ℎ₂,𝐷_𝐸ℎ_2𝑃∗=−𝑤𝑚𝑟𝑘𝛽−𝜇₁−𝑚−𝑐𝑘𝜇₂𝛽²𝑘𝛽². (34) By virtue of (32), it can be obtained that trace𝐷_𝐸ℎ₁adj𝐷_𝐸ℎ₂𝐷_𝐻ℎ₂,𝐷_𝐸ℎ_2𝑃∗≠0. (35)
Furthermore, it can be also calculated that 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ_2𝑃∗=𝑟𝜇₁+𝑚𝑘𝛽−𝜇₁−𝑚𝑤𝑚𝑟𝑘𝛽−𝜇₁−𝑚−𝑐𝑘𝜇₂𝛽²𝑘²𝛽³. (36)
It follows from (32) that 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ_2𝑃∗≠0. (37)
Based on Section IV(A) in [42], ℎ₃(𝐻(𝑡),𝐸(𝑡),𝑣) can be defined as follows: ℎ₃(𝐻(𝑡),𝐸(𝑡),𝑣)=det𝐷_𝐸𝑔=𝑤𝑅(𝑡)−𝑐. (38)
By simple computing, 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝑣ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂𝐷_𝑣ℎ₂𝐷_𝐻ℎ₃𝐷_𝐸ℎ₃𝐷_𝑣ℎ_3𝑃∗=𝑐𝑟𝜇₁+𝑚𝑘𝛽−𝜇₁+𝑚.𝑘𝛽 (39)
According to (32), it derives that 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝑣ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂𝐷_𝑣ℎ₂𝐷_𝐻ℎ₃𝐷_𝐸ℎ₃𝐷_𝑣ℎ_3𝑃∗≠0. (40)
Based on the above analysis, four items (i–iv) can be obtained as follows.(i)It is easy to show that 𝐷_𝐸ℎ₂ has a simple zero eigenvalue: ℎ₁(𝐻(𝑡),𝐸(𝑡),𝑣)_𝑃∗=0,ℎ₂(𝐻(𝑡),𝐸(𝑡),𝑣)_𝑃∗=0, (41) and trace(𝐷_𝐸ℎ₁adj(𝐷_𝐸ℎ₂)(𝐷_𝐻ℎ₂,𝐷_𝐸ℎ₂))|_𝑃^∗≠0 based on (35).(ii)It follows from (37) that 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂ is nonsingular around 𝑃^∗.(iii)By virtue of (40), it can be shown that 𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝑣ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂𝐷_𝑣ℎ₂𝐷_𝐻ℎ₃𝐷_𝐸ℎ₃𝐷_𝑣ℎ₃ is nonsingular around 𝑃^∗; hence rank^{𝐷𝐻ℎ1𝐷𝐸ℎ1𝐷𝑣ℎ1𝐷}_𝐻^ℎ₂^𝐷_𝐸^ℎ₂^𝐷_𝑣^ℎ_{2𝐷𝐻ℎ3𝐷𝐸ℎ3𝐷𝑣ℎ3}=5.(iv)It is easy to show rank(ℎ₁(𝐻(𝑡),𝐸(𝑡),𝑣))=3 and rank(ℎ₂(𝐻(𝑡),𝐸(𝑡),𝑣))=1, which follows 𝐷rank^𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝑣ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂𝐷_𝑣ℎ₂𝐷_𝐻ℎ₃𝐷_𝐸ℎ₃𝐷_𝑣ℎ₃=rankℎ₁(𝐻(𝑡)(𝑡))ℎ,𝐸,𝑣+rank₂(𝐻(𝑡),𝐸(𝑡),𝑣)+1. (42)
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 𝑃^∗ and the bifurcation value is 𝑣=0.
Along with the line of the above proof, for model (4) without time delay, it follows from simple computing that 𝑀=−trace𝐷_𝐸ℎ₁adj𝐷_𝐸ℎ₂𝐷_𝐻ℎ₂,𝐷_𝐸ℎ_2𝑃∗=𝑤𝑚𝑟𝑘𝛽−𝜇₁−𝑚−𝑐𝑘𝜇₂𝛽²𝑘𝛽²,𝑁=𝐷_𝑣ℎ₃−𝐷_𝐻ℎ₃,𝐷_𝐸ℎ₃𝐷_𝐻ℎ₁𝐷_𝐸ℎ₁𝐷_𝐻ℎ₂𝐷_𝐸ℎ₂⁻¹𝐷_𝑣ℎ₁𝐷_𝑣ℎ_2𝑃∗=𝑐𝑟𝜇₁+𝑚𝑘𝛽−𝜇₁−𝑚.𝑘𝛽 (43) It follows from (32) that 𝑀𝑁=𝑚𝑟𝑤𝑘𝛽−𝜇₁−𝑚−𝑐𝑘𝜇₂𝛽²𝑐𝑟𝛽𝜇₁+𝑚𝑘𝛽−𝜇₁−𝑚>0. (44)
Inequality (44) satisfies Theorem 3 of [42]. According to Theorem 3 of [42], when 𝑣 increases through 0, one eigenvalue (denoted by 𝜆₁) 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 𝑃^∗ takes the following form: 𝐽_𝑃^∗=−𝑟𝑆^∗𝑘−𝛽𝑆^∗00𝛽𝐼^∗𝜇0000𝑚−₂+𝐸^∗−𝑅^∗00𝑤𝐸^∗0, (45) according to the leading matrix Ξ(𝑡) in model system (4) and 𝐽_𝑃^∗, the characteristic equation of the model system (4) without time delay around 𝑃^∗ is det𝜆Ξ−𝐽_𝑃^∗=0. (46)
By virtue of simple computation, the characteristic equation is as follows: 𝜆²+𝑟𝑆^∗𝑘𝜆+𝛽²𝑆^∗𝐼^∗=0. (47)
It can be concluded that the rest eigenvalues of model system (4) without time delay (denoted by 𝜆₂ and 𝜆₃) 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 𝑣 increases through 0, and the rest eigenvalues are continuous and nonzero and cannot jump from one half open complex plane to another one as 𝑣 increases through 0. It has been shown that 𝜆₁ moves from ℂ⁻ to ℂ⁺ along the real axis by diverging through infinity. Therefore, 𝜆₂ and 𝜆₃ are continuous and bounded in the ℂ⁻ half plane as 𝑣 increases through 0 and their movement behaviors have no influence on the stability of model system (4) without time delay around the interior equilibrium 𝑃^∗.

According to Table 1 and the stability theory, it can be concluded that model system (4) without time delay is stable around 𝑃^∗ as 𝑣<0 and model system (4) without time delay is unstable around 𝑃^∗ as 𝑣>0. Consequently, a stability switch occurs as 𝑣 increases through 0.

Table 1 

Signs of real parts of eigen values of model (4) without time delay around interior equilibrium 𝑃^∗.

Remark 12. Some preliminaries of singularity induced bifurcation are introduced below. Parameter dependent differential-algebraic hybrid system of the form 𝑥̇(𝑡)=ℎ𝑥(𝑡),𝑦(𝑡),𝜆,ℎ∶𝑅^𝑛×𝑅^𝑚×𝑅^𝑝⟶𝑅^𝑛,0=𝑔𝑥(𝑡),𝑦(𝑡),𝜆,𝑔∶𝑅^𝑛×𝑅^𝑚×𝑅^𝑝⟶𝑅^𝑚, (48) where 𝑥(𝑡), 𝑦(𝑡), 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 ℎ_𝑥−ℎ_𝑦𝑔_𝑦⁻¹𝑔_𝑥 diverging to infinity when Jacobian 𝑔_𝑦 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 Theorem 11 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 𝑣=0 and 𝑣>0, respectively.

According to the leading matrix Ξ(𝑡) in model system (4) and 𝐽_𝑃^∗ in (45) (the Jacobian of model system (4) without time delay around the interior equilibrium 𝑃^∗), it can be calculated that rank(𝐽_𝑃^∗,Ξ𝐽_𝑃^∗,Ξ²𝐽_𝑃^∗,Ξ³𝐽_𝑃^∗)=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 𝑃^∗ in the case of 𝑣=0. Consequently, a state feedback controller can be applied to stabilize the model system (4) without time delay around 𝑃^∗. By using Theorem 3-1.2 in [47], a state feedback controller 𝑢(𝑡)=𝑙(𝐸(𝑡)−𝐸^∗) (𝑙 is a feedback gain and 𝐸^∗ is the component of the interior equilibrium 𝑃^∗) can be applied to stabilize model system (4) without time delay around 𝑃^∗.

Furthermore, the controlled model system (4) without time delay takes the following form: 𝑆̇(𝑡)=𝑟1−𝑆(𝑡)𝑘(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)−𝛽𝑆=𝛽𝑆−𝜇₁𝐼(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)−𝑚𝐼=𝑚𝐼−𝜇₂𝑅(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))(𝑡).−𝐸0=𝐸)(𝑤𝑅−𝑐−𝑣+𝑢 (49)

Theorem 14. When economic interest of harvesting is zero, 𝑣=0, if the feedback gain 𝑙 satisfies the following inequality: 𝑙>max𝑘𝑤𝐸^∗𝑅^∗𝑟𝑆^∗+𝑘𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗𝜇₂+𝐸^∗,𝑟𝑤𝐸^∗𝑅^∗𝑘𝛽²𝐼^∗+𝑟𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗2𝑘𝑟𝜇₂+𝐸^∗+2𝑘²𝛽²𝐼^∗+𝑟²𝑆^∗+4𝑘⁴𝛽⁴𝐼^∗2+8𝜇₂+𝐸^∗𝑟𝑘³𝛽²𝐼^∗+𝑟⁴𝑆^∗2×2𝑟𝑟𝜇₂+𝐸^∗𝑆^∗+𝑘𝛽²𝐼^∗𝑆^∗+𝑘𝜇₂+𝐸^∗2−1, (50) then singularity induced bifurcation is eliminated and model system (49) is stable around 𝑃^∗.

Proof. The Jacobian of the model system (49) evaluated at the interior equilibrium 𝑃^∗ takes the form 𝐽_𝑃^∗=−𝑟𝑆^∗𝑘−𝛽𝑆^∗00𝛽𝐼^∗𝜇0000𝑚−₂+𝐸^∗−𝑅^∗00𝑤𝐸^∗𝑙. (51)
According to the leading matrix Ξ(𝑡) in the model system (4) and 𝐽_𝑃^∗, the characteristic equation of model system (49) around 𝑃^∗ is det(𝜆Ξ−𝐽_𝑃^∗)=0, which can be expressed as follows: 𝜆³+𝐵₁𝜆²+𝐵₂𝜆+𝐵₃=0, (52) where 𝐵₁=𝜇₂+𝐸^∗+𝑟𝑆^∗/𝑘−𝑤𝐸^∗𝑅^∗/𝑙, 𝐵₂=𝛽²𝐼^∗𝑆^∗+(𝑟𝑆^∗/𝑘)(𝜇₂+𝐸^∗−𝑤𝐸^∗𝑅^∗/𝑙), and 𝐵₃=𝛽²𝐼^∗𝑆^∗(𝜇₂+𝐸^∗−𝑤𝐸^∗𝑅^∗/𝑙).
By using the Routh-Hurwitz criteria [43], model (49) is locally stable around 𝑃^∗ if and only if 𝑙 satisfies 𝑙>max𝑘𝑤𝐸^∗𝑅^∗𝑟𝑆^∗+𝑘𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗𝜇₂+𝐸^∗,𝑟𝑤𝐸^∗𝑅^∗𝑘𝛽²𝐼^∗+𝑟𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗2𝑘𝑟𝜇₂+𝐸^∗+2𝑘²𝛽²𝐼^∗+𝑟²𝑆^∗+4𝑘⁴𝛽⁴𝐼^∗2+8𝜇₂+𝐸^∗𝑟𝑘³𝛽²𝐼^∗+𝑟⁴𝑆^∗2×2𝑟𝑟𝜇₂+𝐸^∗𝑆^∗+𝑘𝛽²𝐼^∗𝑆^∗+𝑘(𝜇₂+𝐸^∗)²⁻¹. (53)
Consequently, if the feedback gain satisfies the above inequality, then model system (49) is stable around 𝑃^∗ in the case of zero interest of harvesting.

Let 𝑃^∗(𝑆^∗,𝐼^∗,𝑅^∗,𝐸^∗) denote interior equilibrium of model (4) in the case of positive economic interest of harvesting (𝑣>0), where 𝑆^∗=(𝜇₁+𝑚)/𝛽, 𝐼^∗=𝑟(𝑘𝛽−𝜇₁−𝑚)/𝑘𝛽², 𝐸^∗=𝑣/(𝑤𝑅^∗−𝑐), and 𝑅^∗ satisfies the following equation: 𝑅^∗2+𝐵₂𝑅^∗+𝐵₃=0, (54) where 𝐵₂=(𝑘𝛽²(𝑣−𝑐𝜇₂)+𝑤𝑚𝑟(𝑚+𝜇₁−𝑘𝛽))/𝑘𝑤𝜇₂𝛽², 𝐵₃=𝑐𝑚𝑟(𝑘𝛽−𝑚𝑢₁−𝑚)/𝑘𝑤𝜇₂𝛽².

Based on Routh-Hurwitz criteria [43], (54) has two positive roots if economic interest 𝑣 satisfies the following inequalities: 0<𝑣<min𝑐𝜇₂+𝑤𝑚𝑟𝑘𝛽−𝑚−𝜇₁𝑘𝛽²,(1−𝑐)𝑐𝜇₂+𝑤𝑚𝑟(𝑐+1)𝑘𝛽−𝜇₁−𝑚𝑘𝛽²∶=𝑣.̃ (55)

As analyzed above, there are two interior equilibria (denoted by 𝑃^∗₁ and 𝑃^∗₂) when 0<𝑣^∗<𝑣̃. In this subsection, we only design the controller for the model (4) around the interior equilibrium 𝑃^∗₁. Some symmetric results about 𝑃^∗₂ can be also obtained, and 𝑃^∗₁ is denoted as 𝑃^∗ for simplicity in the following part.

Theorem 15. When economic interest of harvesting is positive, 0<𝑣^∗<𝑣̃, if feedback gain 𝑙̃ of controller 𝑢(𝑡)=𝑙(̃𝐸(𝑡)−𝐸^∗) satisfies following inequality: 𝑙̃>max𝑘𝑤𝐸^∗𝑅^∗𝑟𝑆^∗+𝑘𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗𝜇₂+𝐸^∗,𝑟𝑤𝐸^∗𝑅^∗𝑘𝛽²𝐼^∗+𝑟𝜇₂+𝐸^∗,𝑤𝐸^∗𝑅^∗2𝑘𝑟𝜇₂+𝐸^∗+2𝑘²𝛽²𝐼^∗+𝑟²𝑆^∗+4𝑘⁴𝛽⁴𝐼^∗2+8𝜇₂+𝐸^∗𝑟𝑘³𝛽²𝐼^∗+𝑟⁴𝑆^∗2×2𝑟𝑟𝜇₂+𝐸^∗𝑆^∗+𝑘𝛽²𝐼^∗𝑆^∗+𝑘𝜇₂+𝐸^∗2, (56) then model system (49) is stable around the interior equilibrium 𝑃^∗(𝑆^∗,𝐼^∗,𝑅^∗,𝐸^∗).

Proof. The proof is similar to the proof of Theorem 14 of this paper.

Remark 16. It follows from (55) and Theorem 15 that economic interest of commercial harvesting should be regulated within certain interval 𝑣∈(0,𝑣)̃, 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 harvesting 0<𝑣^∗<𝑣̃ where 𝑣̃ is defined in (55), there are two interior equilibria 𝑃^∗₁ and 𝑃^∗₂ for model system (4) with respect to the positive economic interest 𝑣^∗.

In this subsection, we only investigate dynamical behavior of model system (4) around the interior equilibrium 𝑃^∗₁. Some symmetric results about the interior equilibrium 𝑃^∗₂ can be also obtained, and 𝑃^∗₁ is denoted as 𝑃^∗ for simplicity. According to Jacobian evaluated at the interior equilibrium 𝑃^∗ and the leading matrix Ξ(𝑡) in model system (4), we can obtain the characteristic equation of model system (4) around 𝑃^∗, which can be expressed as follows: 𝑟𝑆𝜆+^∗𝑘𝛽𝑆^∗𝑒^−𝜆𝜏𝐼00−𝛽^∗𝜆−𝛽𝑆^∗𝑒^−𝜆𝜏+𝜇₁+𝑚000−𝑚𝜆+𝜇₁+𝐸^∗𝑅^∗00−𝑤𝐸^∗−𝑣^∗𝐸^∗=0. (57)  ⟹𝑀(𝜆)+𝑁(𝜆)𝑒^−𝜆𝜏=0, (58) where 𝑀(𝜆)=𝜆³+𝑚₁𝜆²+𝑚₂𝜆+𝑚₃,𝑁(𝜆)=𝑛₁𝜆²+𝑛₂𝜆+𝑛₃,𝑚₁=𝜇₂+𝜇₁𝑟+𝑚1+−𝑘𝛽𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝑚₂=𝜇₁+𝑚𝑟𝜇𝑘𝛽¹+𝑚+𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²+𝑘𝛽𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝑚₃=𝑟𝜇₁+𝑚²𝜇𝑘𝛽²−𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝑛₁=−𝜇₁𝑛−𝑚,₂=−𝜇₁×𝜇+𝑚₂+2𝑟𝜇₁+𝑚−𝑘𝑟𝛽−𝑘𝛽𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝑛₃=−𝑟𝜇₁2𝜇+𝑚¹+𝑚×𝜇𝑘𝛽−1₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐². (59)

Now substituting 𝜆=𝑖𝜎 (𝜎 is a positive real number) into (58) and separating the real and imaginary parts, two transcendental equations can be obtained as follows: 𝜎³−𝑚₂𝜎=𝑛₂𝜎cos(𝜎𝜏)−𝑛₃−𝑛₁𝜎²sin(𝜎𝜏),𝑚₁𝜎²−𝑚₃=𝑛₃−𝑛₁𝜎²cos(𝜎𝜏)+𝑛₂𝜎sin(𝜎𝜏). (60)(61)

By squaring and adding (60) and (61), it can be calculated that 𝑛₃−𝑛₁𝜎²²+𝑛²₂𝜎²=𝑚₁𝜎²−𝑚₃²+𝜎³−𝑚₂𝜎², (62)  ⟹𝜎⁶+𝐶₁𝜎⁴+𝐶₂𝜎²+𝐶₃=0, (63) where 𝐶₁=𝑚²₁−2𝑚₂−𝑛²₁, 𝐶₂=𝑚²₂−2𝑚₁𝑚₃+2𝑛₁𝑛₃−𝑛²₂, 𝐶₃=𝑚²₃−𝑛²₃.

According to the values of 𝐶_𝑗,(𝑗=1,2,3) and the Routh-Hurwitz criteria [43], a simple assumption that (58) has at least one positive real root 𝜎₀ is 𝐶₃<0, which derives that 𝑘𝛽>3(𝜇₁+𝑚). Hence, under this assumption, (58) will have a pair of purely imaginary roots of the form ±𝑖𝜎₀.

By eliminating sin(𝜎𝜏) from (60) and (61), it can be calculated that the 𝜏^∗_𝑗 corresponding to 𝜎₀ is as follows: 𝜏^∗_𝑗=1𝜎₀𝑛×arccos²𝜎²₀𝜎²₀−𝑚₂+𝑛₃−𝑛₁𝜎²₀𝑚₁𝜎²₀−𝑚₃𝑛₃−𝑛₁𝜎²₀²+𝑛₂𝜎₀²+2𝑗𝜋𝜎₀, (64) where 𝑗=0,1,2,….

By using Butler’s lemma [48], model system (4) is locally stable around 𝑃^∗ for 𝜏<𝜏^∗₀. Subsequently, conditions for existence of Hopf bifurcation in [39] are utilized to investigate whether Hopf bifurcation occurs as 𝜏 increases through 𝜏^∗_𝑗.

Theorem 17. If 𝑘𝛽>3(𝜇₁+𝑚), then model system (4) undergoes Hopf bifurcation around the interior equilibrium 𝑃^∗ when 𝜏=𝜏^∗_𝑗,  𝑗=0,1,2,…. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium 𝑃^∗ when 𝜏>𝜏^∗₀ and ‖𝜏−𝜏^∗₀‖≪1.

Proof. As mentioned above, let 𝜆=𝑖𝜎₀ represent the purely imaginary root of (58). It follows from (58) that |𝑀(𝑖𝜎₀)|=|𝑁(𝑖𝜎₀)|, which determines a set of possible values of 𝜎₀.
In the following part, we determine the direction of motion of 𝜆=𝑖𝜎₀ as 𝜏 is varied; namely, we determine d()Θ=signRe𝜆d𝜏_{𝜆=𝑖𝜎0}=signRed𝜆d𝜏⁻¹_{𝜆=𝑖𝜎0}. (65)
By differentiating (58) with respect to 𝜏, it can be obtained that d𝜆d𝜏⁻¹=3𝜆²+2𝑚₁𝜆+𝑚₂𝜆𝑒^−𝜆𝜏𝑛₁𝜆²+𝑛₂𝜆+𝑛₃+2𝑛₁𝜆+𝑛₂𝜆𝑛₁𝜆²+𝑛₂𝜆+𝑛₃−𝜏𝜆=3𝜆²+2𝑚₁𝜆+𝑚₂−𝜆𝜆³+𝑚₁𝜆²+𝑚₂𝜆+𝑚₃+2𝑛₁𝜆+𝑛₂𝜆𝑛₁𝜆²+𝑛₂𝜆+𝑛₃−𝜏𝜆=2𝜆³+𝑚₁𝜆²−𝑚₃−𝜆²𝜆³+𝑚₁𝜆²+𝑚₂𝜆+𝑚₃+𝑛₁𝜆²−𝑛₃𝜆²𝑛₁𝜆²+𝑛₂𝜆+𝑛₃−𝜏𝜆. (66)
From (62) and the above equation, it can be obtained that Θ=signRed𝜆d𝜏⁻¹_{𝜆=𝑖𝜎0}=1𝜎²₀𝑚sign³+𝑚₁𝜎²₀𝑚₁𝜎²₀−𝑚₃+2𝜎⁴₀𝜎²₀−𝑚₂𝑚₁𝜎²₀−𝑚₃²+𝜎³₀−𝑚₂𝜎₀²+𝑛₁𝜎²₀+𝑛₃𝑛₃−𝑛₁𝜎²₀𝑛₃−𝑛₁𝜎²₀²+𝑛₂𝜎₀²=sign𝑚₃+𝑚₁𝜎²₀𝑚₁𝜎²₀−𝑚₃+2𝜎⁴₀𝜎²₀−𝑚₂+𝑛₁𝜎²₀+𝑛₃𝑛₃−𝑛₁𝜎²₀×𝑚₁𝜎²₀−𝑚₃²+𝜎³₀−𝑚₂𝜎₀²⁻¹=sign2𝜎⁶₀+𝑚²₁−2𝑚₂−𝑛²₁𝜎⁴₀+𝑛²₃−𝑚²₃𝑚₁𝜎²₀−𝑚₃²+𝜎³₀−𝑚₂𝜎₀²,=sign2𝜎⁶₀+𝐶₁𝜎⁴₀−𝐶₃𝑚₁𝜎²₀−𝑚₃²+𝜎³₀−𝑚₂𝜎₀². (67)
According to the values of 𝐶_𝑗,(𝑗=1,2,3) given in (58) of this paper, it is easy to show that 𝐶₁=𝑚²₁−2𝑚₂−𝑛²₁=𝑟²(𝜇₁+𝑚)²/𝑘²𝛽²+[𝜇₂−𝑐𝑣^∗/(𝑤𝑅−𝑐)²]²>0.
Furthermore, if 𝑘𝛽>3(𝜇₁+𝑚), then it can be shown that 𝐶₃>0. Hence, it can be concluded that 2𝜎⁶₀+𝐶₁𝜎⁴₀−𝐶₃>0, which derives sign[d(Re𝜆)/d𝜏]_𝜏=𝜏^∗_{𝑗,𝜎=𝜎0}>0. Consequently, the transversality condition holds and Hopf bifurcation occurs at 𝜎=𝜎₀, 𝜏=𝜏^∗_𝑗. Furthermore, an attracting invariant closed curve bifurcates from interior equilibrium 𝑃^∗ when 𝜏>𝜏^∗₀ and ‖𝜏−𝜏^∗₀‖≪1.

Remark 18. sign[d(Re𝜆)/d𝜏]_𝜏=𝜏^∗_𝑗>0 signifies that there exists at least one eigenvalue with positive real part for 𝜏=𝜏^∗_𝑗, 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<𝑣^∗<𝑣̃ (𝑣 is defined in (55)), it follows from implicit function theorem [41] and the fourth equation of model system (4) that 𝐸(𝑡)=𝑣^∗/(𝑤𝑅(𝑡)−𝑐). Furthermore, model system (4) can be transformed into the following form: 𝑆̇(𝑡)=𝑟1−𝑆(𝑡)𝑘(𝑡)𝐼(),𝐼̇(𝑡)(𝑡)𝐼()−𝛽𝑆𝑡−𝜏=𝛽𝑆𝑡−𝜏−𝜇₁𝐼(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)−𝑚𝐼=𝑚𝐼−𝜇₂𝑅(𝑡)−𝑣^∗𝑅(𝑡)(𝑡).𝑤𝑅−𝑐 (68)

Firstly, some transformations associated with component (𝑆^∗,𝐼^∗,𝑅^∗) of interior equilibrium 𝑃^∗ are given as follows: 𝑦₁=𝑆−𝑆^∗,𝑦₂=𝐼−𝐼^∗,𝑦₃=𝑅−𝑅^∗,𝑦_𝑖(𝑡)=𝑦_𝑖(𝜏𝑡),𝜏=𝜌+𝜏_𝑗,𝜌∈ℝ=(−∞,+∞). (69) 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 𝐶=𝐶([−1,0],ℝ³) as 𝑦̇(𝑡)=𝐿_𝜌𝑦_𝑡+𝑓𝜌,𝑦_𝑡, (70) where 𝐶=𝐶([−1,0],ℝ³) is the Banach space of continuous functions mapping the interval [−𝜏,0] into ℝ³, 𝑦(𝑡)=(𝑦₁(𝑡),𝑦₂(𝑡),𝑦₃(𝑡))^𝑇∈ℝ³, 𝑦_𝑡(𝜃)=𝑦(𝑡+𝜃) for 𝜃∈[−𝜏,0] and 𝐿_𝜌∶𝐶→ℝ³, 𝑓∶ℝ×𝐶→ℝ³ are defined as follows, respectively: 𝐿_𝜌𝜙=𝜏_𝑗×−𝑟𝑆+𝜌_∗𝑘𝛽𝐼00_∗−𝜇₁+𝑚00𝑚−𝜇₂+𝑐𝑣^∗𝑤𝑅^∗−𝑐²×𝜙₁(0)𝜙₂(0)𝜙₃(0)+𝜏_𝑗𝑆+𝜌0−𝛽^∗00𝛽𝑆^∗0×𝜙000₁()𝜙−1₂()𝜙−1₃(),−1𝑓𝜌,𝜙=𝜏_𝑗−𝑟+𝜌𝑘𝜙²₁(0)−𝛽𝜙₁(0)𝜙₂(−1)𝛽𝜙₁(0)𝜙₂(−1_∞𝑗=1)𝑓_𝑗⁽³⁾𝜙^𝑗₃(0), (71)(72) where 𝑓_𝑗⁽³⁾=(d^𝑖(−𝑣^∗𝑅/(𝑤𝑅−𝑐))/d𝑅^𝑖)|_(𝑆^∗_,𝐼^∗_,𝑅^∗₎, 𝑗=1,2,3.

It is easy to show that 𝐿_𝜌 is a continuous linear function mapping 𝐶 into ℝ³. According to Riesz representation theorem [40], there exists a 3×3 matrix function 𝜂(𝜃,𝜌) of bounded variation for 𝜃∈[−1,0] such that 𝐿_𝜌𝜙=⁰₋₁d𝜂𝜃,𝜌𝜙(𝜃), (73) where 𝜙∈𝐶([−1,0],ℝ³).

In fact, we can choose 𝜂𝜃,𝜌=𝜏_𝑗×−𝑟𝑆+𝜌_∗𝑘𝛽𝐼00_∗−𝜇₁+𝑚00𝑚−𝜇₂+𝑐𝑣^∗𝑤𝑅^∗−𝑐²(𝜃)−𝜏×𝛿_𝑗𝑆+𝜌0−𝛽^∗00𝛽𝑆^∗0𝛿(),000𝜃+1 (74) where 𝛿 denotes the Dirac delta function.

If 𝜙 is any given function in 𝐶([−1,0],ℝ³) and 𝑦(𝜙) is the unique solution of the linearized equation 𝑦(̇𝑡)=𝐿_𝜌(𝑦_𝑡) of (70) with initial function 𝜙 at zero, then the solution operator 𝑇(𝑡)∶𝐶→𝐶 is defined by 𝑇(𝑡)𝜙=𝑦_𝑡𝜙,𝑡≥0. (75)

It follows from Lemma  7.1.1 in [39] that 𝑇(𝑡),  𝑡≥0 is a strongly continuous semigroup of linear transformation on [0,+∞) and the infinitesimal generator 𝐴_𝜌 of 𝑇(𝑡), 𝑡≥0 is as follows: 𝐴_𝜌𝜙=d𝜙(𝜃)[)d𝜃,𝜃∈−1,0₀₋₁d𝜂𝜌,𝑠𝜙(𝑠),𝜃=0, (76) for 𝜙∈𝐶¹([−1,0],ℝ³), the space of functions mapping the interval [−1,0] into ℝ³ which have continuous first derivative and also define 𝑅𝜌𝜙=[)𝑓0,𝜃∈−1,0𝜌,𝜙,𝜃=0 (77) then model system (70) is equivalent to 𝑦̇_𝑡=𝐴𝜌𝑦_𝑡+𝑅𝜌𝑦_𝑡. (78)

For 𝜓∈𝐶¹([0,1],(ℝ³)^∗), the space of functions mapping interval [0,1] into the three-dimensional row vectors which have continuous first derivative, define 𝐴^∗𝜙(𝑠)=−d𝜓(𝑠)(]d𝑠,𝑠∈0,1₀₋₁d𝜂^𝑇(𝑡,0)𝜓(−𝑡),𝑠=0, (79) and a bilinear inner product 𝜓(𝑠),𝜙(𝜃)=𝜓(0)𝜙(0)−⁰₋₁^𝜃_𝜉=0𝜓(𝜉−𝜃)d𝜂(𝜃)𝜙(𝜉)d𝜉, (80) where 𝜂(𝜃)=𝜂(𝜃,0). It follows from the above analysis 𝐴(0) and 𝐴^∗ are adjoint operators.

By virtue of discussion in Section 4.2.1, ±𝑖𝜔₀𝜏_𝑗 are eigenvalues of 𝐴(0). Hence, they are also eigenvalues of 𝐴^∗. In the following, eigenvectors of 𝐴(0) and 𝐴^∗ are corresponding to 𝑖𝜔₀𝜏_𝑗 and −𝑖𝜔₀𝜏_𝑗, respectively.

Suppose 𝑞(𝜃)=(1,𝑎,𝑏)^𝑇𝑒^{𝑖𝜔0𝜏𝑗𝜃} is the eigenvectors of 𝐴(0) corresponding to 𝑖𝜔₀𝜏_𝑗, which derives that 𝐴(0)𝑞(𝜃)=𝑖𝜔₀𝜏_𝑗𝑞(𝜃). By using the definition of 𝐴(0), (71) and (72), it gives that −𝑟𝑆^∗𝑘𝛽𝐼00_∗−𝜇₁+𝑚00𝑚−𝜇₂+𝑐𝑣^∗𝑤𝑅^∗−𝑐²𝑞(0)+0−𝛽𝑆^∗00𝛽𝑆^∗0𝑞()000−1=𝑖𝜔₀𝑞(0). (81)

For 𝑞(−1)=𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑗}, then it can be obtained that 𝑟𝜇𝑎=−¹+𝑚+𝑖𝑘𝛽𝜔₀𝑘𝛽𝜇₁+𝑚𝑒^{−𝑖𝜔0𝜏𝑗},𝑚𝑤𝑅𝑏=^∗−𝑐²𝑟𝑚+𝜇₁+𝑖𝑘𝛽𝜔₀𝑘𝛽𝜇₁+𝑚𝑐𝑣^∗−𝜇₂+𝑖𝜔₀𝑤𝑅^∗−𝑐²𝑒^{−𝑖𝜔0𝜏𝑗}. (82)

Similarly, it follows from simple computation that eigenvector 𝑞^∗(𝑠)=𝐽(1,𝑎^∗,𝑏^∗)𝑒^{𝑖𝜔0𝜏𝑗𝑠} of 𝐴^∗ is corresponding to −𝑖𝜔₀𝜏_𝑗, where 𝑎^∗𝑟𝜇=−¹+𝑚−𝑖𝑘𝛽𝜔₀𝑘𝛽𝜇₁+𝑚𝑒^{𝑖𝜔0𝜏𝑗},𝑏^∗=𝑚𝑤𝑅^∗−𝑐²𝑟𝜇₁+𝑚−𝑖𝑘𝛽𝜔₀𝑘𝛽𝜇₁+𝑚𝑐𝑣^∗−𝜇₂−𝑖𝜔₀𝑤𝑅^∗−𝑐²𝑒^{𝑖𝜔0𝜏𝑗}. (83)

In order to assume ⟨𝑞^∗(𝑠),𝑞(𝜃)⟩=1, we need to determine the value of 𝐽 in the following part.

By virtue of (80), it derives that ⟨𝑞^∗(𝑠)(𝜃)⟩=,𝑞𝐽1,𝑎^∗,𝑏^∗(1,𝑎,𝑏)^𝑇−⁰₋₁^𝜃_𝜉=0𝐽1,𝑎^∗,𝑏^∗𝑒^{−𝑖𝜔0𝜏𝑗(𝜉−𝜃)}(𝜃)×()d𝜂1,𝑎,𝑏^𝑇𝑒^{𝑖𝜔0𝜏𝑗𝜉}=d𝜉𝐽1+𝑎𝑎^∗+𝑏𝑏^∗−⁰₋₁1,𝑎^∗,𝑏^∗𝜃𝑒^{𝑖𝜔0𝜏𝑗𝜃}d𝜂(𝜃)(1,𝑎,𝑏)^𝑇=𝐽1+𝑎𝑎^∗+𝑏𝑏^∗+𝑘𝑏^∗𝜏_𝑗𝑅^∗+𝑏𝑆^∗𝑒^{𝑖𝜔0𝜏𝑗}. (84)

Hence, we can choose 𝐽 as follows: 1𝐽=1+𝑎𝑎^∗+𝑏𝑏^∗+𝑘𝑏^∗𝜏_𝑗𝑅^∗+𝑏𝑆^∗𝑒^{𝑖𝜔0𝜏𝑗}. (85)

Next, we will compute the coordinate to describe the centre manifold 𝐶₀ at 𝜌=0. Let 𝑦_𝑡 be the solution of (78) when 𝜌=0.

Define 𝑧(𝑡)=⟨𝑞^∗,𝑦_𝑡⟩,𝑊(𝑡,𝜃)=𝑦_𝑡(𝜃)−2Re𝑧(𝑡)𝑞(𝜃). (86)

On the center manifold 𝐶₀, it derives that 𝑊(𝑡,𝜃)=𝑊(𝑧(𝑡),𝑧(𝑡),𝜃), (87) where 𝑊(𝑧(𝑡),𝑧(𝑡),𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧+𝑊₀₂𝑧²2+⋯, (88)𝑧 and 𝑧 are local coordinates for center manifold 𝐶₀ in the direction of 𝑞^∗ and 𝑞^∗.

It is noted that 𝑊 is real if 𝑦_𝑡 is real, and we only consider real solutions. For solution 𝑦_𝑡∈𝐶₀ of (78), since 𝜌=0, it derives that 𝑧̇(𝑡)=𝑖𝜔₀𝜏_𝑗𝑧+𝑞^∗(0)𝑓0,𝑊(𝑧,𝑧,0)+2Re𝑧𝑞(𝜃)≜𝑖𝜔₀𝜏_𝑗𝑧+𝑞^∗(0)𝑓₀(𝑧,𝑧). (89)

The above equation can be rewritten as follows: 𝑧̇(𝑡)=𝑖𝜔₀𝜏_𝑗𝑧(𝑡)+𝑔(𝑧,𝑧), (90) where 𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+𝑔₂₁𝑧²𝑧2+⋯. (91)

It follows from (86) and (88) that 𝑦_𝑡(𝜃)=𝑊(𝑡,𝜃)+2Re𝑧(𝑡)𝑞(𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₁₁(𝜃)𝑧𝑧+𝑊₀₂(𝜃)𝑧²2+(1,𝑎,𝑏)^𝑇𝑒^{𝑖𝜔0𝜏𝑗𝜃}𝑧+1,𝑎,𝑏^𝑇𝑒^{−𝑖𝜔0𝜏𝑗𝜃}𝑧+⋯. (92)

By virtue of (72), (91), and (92), it derives that𝑔(𝑧,𝑧)=𝑞^∗(0)𝑓₀(𝑧,𝑧)=𝑞^∗(0)𝑓0,𝑦_𝑡=𝜏_𝑗𝐽−𝑟𝑘𝑦²_1𝑡(0)−𝛽𝑦_1𝑡(0)𝑦_2𝑡(−1)𝛽𝑦_1𝑡(0)𝑦_2𝑡(−1_∞𝑗=1)𝑓_𝑗⁽³⁾𝑦^𝑗_3𝑡(0)=𝜏_𝑗𝐽−𝑟𝑘𝑦²_1𝑡(0)−𝛽𝑦_1𝑡(0)𝑦_2𝑡()+−1𝑎^∗𝛽𝑦_1𝑡(0)𝑦_2𝑡(−1)+^∞_𝑗=1𝑏^∗𝑓_𝑗⁽³⁾𝑦^𝑗_3𝑡(0)=𝜏_𝑗𝐽−𝑟𝑘𝑦²_1𝑡(0)+𝑎^∗−1𝛽𝑦_1𝑡(0)𝑦_2𝑡()+−1_∞𝑗=1𝑏^∗𝑓_𝑗⁽³⁾𝑦^𝑗_3𝑡(0)𝑟=−𝑘𝜏_𝑗𝐽𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2+𝑜|(𝑧,𝑧)|³²+𝑎^∗−1𝛽𝜏_𝑗𝐽×𝑧+𝑧+𝑊⁽¹⁾₂₀(0)𝑧²2+𝑊⁽¹⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2+𝑜|𝑧,𝑧|³×𝑒^{−𝑖𝜔0𝜏𝑗}𝑧+𝑒^{𝑖𝜔0𝜏𝑗}𝑧+𝑊⁽²⁾₂₀()𝑧−1²2+𝑊⁽²⁾₁₁(−1)𝑧𝑧+𝑊⁽²⁾₀₂(−1)𝑧²2+𝑜|𝑧,z|³+𝜏_𝑗𝐽𝑏^∗𝑓₁⁽³⁾𝑧+𝑧+𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2+𝑜|𝑧,𝑧|³+𝑓₂⁽³⁾𝑧+𝑧+𝑊⁽³⁾₂₀(0)𝑧²2+𝑊⁽³⁾₁₁(0)𝑧𝑧+𝑊⁽¹⁾₀₂(0)𝑧²2+𝑜|𝑧,𝑧|³+⋯. (93)

By comparing the coefficients with (91), it gives that 𝑔₂₀=2𝜏_𝑗𝐽−𝑟𝑘+𝑎^∗−1𝛽𝑒^{−𝑖𝜔0𝜏𝑗}+𝑏^∗𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾,𝑔₁₁=𝜏_𝑗𝐽−2𝑟𝑘+2𝑎^∗−1𝛽cos𝜔₀𝜏_𝑗+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0)+2𝑓₂⁽³⁾,𝑔₀₂=2𝜏_𝑗𝐽−𝑟𝑘+𝑎^∗−1𝛽𝑒^{−𝑖𝜔0𝜏𝑗}+𝑏^∗𝑓₁⁽³⁾𝑊⁽³⁾₀₂(0)2+𝑓₂⁽³⁾,𝑔₂₁=2𝜏_𝑗𝐽𝛽𝑎^∗×𝑊−1₍₂₎₁₁()+𝑊−1⁽³⁾²⁰(−1)+𝑊⁽¹⁾₂₀(0)𝑒^{𝑖𝜔0𝜏𝑗}2+𝑊⁽¹⁾₁₁(0)𝑒^{−𝑖𝜔0𝜏𝑗}−𝑟𝑘𝑊⁽¹⁾₂₀(0)+2𝑊⁽¹⁾₁₁(0)+𝑓₂⁽³⁾𝑊⁽³⁾₂₀(0)+2𝑊⁽³⁾₁₁(0). (94)

Since 𝑔₂₁ is associated with 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃), further attempts should be carried out to compute 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃).

By virtue of (78) and (86), we have 𝑊̇=𝑦̇_𝑡−𝑧𝑞̇−𝑧̇𝑞=𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(𝜃),𝜃∈[−1,0)𝐴𝑊−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀(,𝜃=0.≜𝐴𝑊+𝐻𝑧,𝑧,𝜃), (95) where 𝐻(𝑧,𝑧,𝜃)=𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2+⋯. (96)

By substituting the corresponding series into (95) and comparing the coefficients, we have 𝐴−2𝑖𝜔₀𝜏_𝑗𝑊₂₀(𝜃)=−𝐻₂₀(𝜃),𝐴𝑊₁₁(𝜃)=−𝐻₁₁(𝜃),…. (97)

It follows from (95) that for 𝜃∈[−1,0)𝐻(𝑧,𝑧,𝜃)=−𝑞^∗(0)𝑓₀𝑞(𝜃)−𝑞^∗(0)𝑓₀𝑞(𝜃)=−𝑔(𝑧,𝑧)𝑞(𝜃)−𝑔(𝑧,𝑧)𝑞(𝜃). (98)

By comparing coefficients in (96) with those in (94), it derives that 𝐻₂₀(𝜃)=−𝑔₂₀𝑞(𝜃)−𝑔₀₂𝑞(𝜃),𝐻₁₁(𝜃)=−𝑔₁₁𝑞(𝜃)−𝑔₁₁𝑞(𝜃). (99)(100)

Based on the definition of 𝐴 and (97) and (99), it can be obtained that 𝑊̇₂₀(𝜃)=2𝑖𝜔₀𝜏_𝑗𝑊₂₀(𝜃)+𝑔₂₀𝑞(𝜃)+𝑔₀₂𝑞(𝜃). (101)

For 𝑞(𝜃)=(1,𝑎,𝑏)^𝑇𝑒^{𝑖𝜔0𝜏𝑗𝜃}, 𝑊₂₀(𝜃)=𝑖𝑔₂₀𝜔₀𝜏_𝑗𝑞(0)𝑒^{𝑖𝜔0𝜏𝑗𝜃}+𝑖𝑔₀₂3𝜔₀𝜏_𝑗𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑗𝜃}+𝐺₁𝑒^{2𝑖𝜔0𝜏𝑗𝜃}, (102) where 𝐺₁=(𝐺₁⁽¹⁾,𝐺₁⁽²⁾,𝐺₁⁽³⁾) is a constant vector.

Similarly, it follows from (97) and (100) that 𝑊₁₁(𝜃)=−𝑖𝑔₁₁𝜔₀𝜏_𝑗𝑞(0)𝑒^{𝑖𝜔0𝜏𝑗𝜃}+𝑖𝑔₁₁𝜔₀𝜏_𝑗𝑞(0)𝑒^{−𝑖𝜔0𝜏𝑗𝜃}+𝐺₂, (103) where 𝐺₂=(𝐺₂⁽¹⁾,𝐺₂⁽²⁾,𝐺₂⁽³⁾) is a constant vector.

Subsequently, values of 𝐺₁ and 𝐺₂ should be computed. By using the definition of 𝐴 and (95), we have ⁰₋₁d𝜂𝑤₂₀(𝜃)=2𝑖𝜔₀𝜏_𝑗𝑊₂₀(0)−𝐻₂₀(0),₀₋₁d𝜂(𝜃)𝑊₁₁(𝜃)=−𝐻₁₁(0), (104)(105) where 𝜂(𝜃)=𝜂(0,𝜃). Based on (95), it derives that in the case of 𝜃=0, 𝐻(𝑧,𝑧,0)=−2Re𝑞^∗(0)𝑓₀𝑞(0)+𝑓(0)=−𝑞^∗(0)𝑓₀𝑞(0)−𝑞^∗(0)𝑓₀𝑞(0)+𝑓₀=−𝑔(𝑧,𝑧)𝑞(0)−𝑔(𝑧,𝑧)𝑞(0)+𝑓₀, (106) which follows that 𝐻₂₀(𝜃)𝑧²2+𝐻₁₁(𝜃)𝑧𝑧+𝐻₀₂(𝜃)𝑧²2(0)𝑔+⋯=−𝑞₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2−+⋯𝑞(0)𝑔₂₀𝑧²2+𝑔₁₁𝑧𝑧+𝑔₀₂𝑧²2+⋯+𝑓₀. (107)

By virtue of (72), it gives that 𝑓₀=𝜏_𝑘−𝑟𝑘𝑦²_1𝑡(0)−𝛽𝑦_1𝑡(0)𝑦_2𝑡(−1)𝛽𝑦_1𝑡(0)𝑦_2𝑡(−1_∞𝑗=1)𝑓_𝑗⁽³⁾𝑦^𝑗_3𝑡(0). (108)

By virtue of (86), it can be obtained that 𝑦_𝑡(𝜃)()𝑧(𝑡)𝑞(𝜃()(𝑡)𝑞(𝜃)+=𝑊𝑡,𝜃+2Re)=𝑊𝑡,𝜃+𝑧𝑧(𝑡)𝑞(𝜃)=𝑊₂₀(𝜃)𝑧²2+𝑊₂₁(𝜃)𝑧(𝑡)𝑞(𝜃)+𝑧+𝑧𝑧(𝑡)𝑞(𝜃)+⋯. (109) Then we have 𝑓₀=𝜏_𝑗−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾𝑧²+𝜏_𝑘−2𝑟𝑘−2𝛽cos𝜔₀𝜏_𝑗2𝛽cos𝜔₀𝜏_𝑗2𝑓₂⁽³⁾+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0)𝑧𝑧+⋯. (110)

According to (107) and (110), we have 𝐻₂₀(0)=−𝑔₂₀𝑞(0)−𝑔₀₂𝑞(0)+2𝜏_𝑗−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾,𝐻₁₁(0)=−𝑔₁₁𝑞(0)−𝑔₁₁𝑞(0)+𝜏_𝑗−2𝑟𝑘−2𝛽cos𝜔₀𝜏_𝑗2𝛽cos𝜔₀𝜏_𝑗2𝑓₂⁽³⁾+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0). (111)(112)

Since 𝑖𝜔₀𝜏_𝑗 is the eigenvalue of 𝐴(0) and 𝑞(0) is the corresponding eigenvector, we obtain that 𝑖𝜔₀𝜏_𝑗𝐼−⁰₋₁𝑒^{𝑖𝜔0𝜏𝑗𝜃}d𝜂(𝜃)𝑞(0)=0,−𝑖𝜔₀𝜏_𝑗𝐼−⁰₋₁𝑒^{−𝑖𝜔0𝜏𝑗𝜃}d𝜂(𝜃)𝑞(0)=0, (113) where 𝐼 is identity matrix.

By substituting (102) and (111) into (104), it can be obtained that 2𝑖𝜔₀𝜏_𝑗𝐼−⁰₋₁𝑒^{2𝑖𝜔0𝜏𝑗𝜃}d𝜂(𝜃)𝐺₁=2𝜏_𝑗−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾, (114) which can be rewritten as follows: 2𝑖𝜔₀+𝑟𝑆^∗𝑘𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}0−𝛽𝐼^∗2𝑖𝜔₀+𝜇₁+𝑚−𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}00−𝑚2𝑖𝜔₀+𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²𝐺₁−𝑟=2𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾. (115)

Based on Grammar Law [43], 𝐺₁⁽¹⁾, 𝐺₁⁽²⁾, and 𝐺₁⁽³⁾ can be obtained as follows:𝐺₁⁽¹⁾=2𝑈₁−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}0𝛽𝑒^{−𝑖𝜔0𝜏𝑗}2𝑖𝜔₀+𝜇₁+𝑚−𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}0𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾−𝑚2𝑖𝜔₀+𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝐺₁⁽²⁾=2𝑈₁2𝑖𝜔₀+𝑟𝑆^∗𝑘−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}0−𝛽𝐼^∗𝛽𝑒^{−𝑖𝜔0𝜏𝑗}00𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾2𝑖𝜔₀+𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝐺₁⁽³⁾=2𝑈₁2𝑖𝜔₀+𝑟𝑆^∗𝑘𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}−𝑟𝑘−𝛽𝑒^{−𝑖𝜔0𝜏𝑗}−𝛽𝐼^∗2𝑖𝜔₀+𝜇₁+𝑚−𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}𝛽𝑒^{−𝑖𝜔0𝜏𝑗}0−𝑚𝑓₁⁽³⁾𝑊⁽³⁾₂₀(0)2+𝑓₂⁽³⁾, (116) where 𝑈₁=2𝑖𝜔₀+𝑟𝑆^∗𝑘𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}0−𝛽𝐼^∗2𝑖𝜔₀+𝜇₁+𝑚−𝛽𝑆^∗𝑒^{−2𝑖𝜔0𝜏𝑗}00−𝑚2𝑖𝜔₀+𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐². (117)

Similarly, substituting (103) and (112) into (105), it can be obtained that 𝑟𝑆^∗𝑘𝛽𝑆^∗0−𝛽𝐼^∗𝜇₁+𝑚−𝛽𝑆^∗00−𝑚𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²𝐺₂−=22𝑟𝑘−2𝛽cos𝜔₀𝜏_𝑗2𝛽cos𝜔₀𝜏_𝑗2𝑓₂⁽³⁾+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0). (118)

Based on Grammar Law [43], 𝐺₂⁽¹⁾, 𝐺₂⁽²⁾, and 𝐺₂⁽³⁾ can be obtained as follows: 𝐺₂⁽¹⁾=1𝑈₂^{−2𝑟𝑘−2𝛽cos𝜔0𝜏𝑗𝛽𝑆∗02𝛽cos𝜔0𝜏𝑗𝜇1𝑆+𝑚−𝛽∗0}_{2𝑓2(3)+𝑓1(3)𝑊(3)11(0)−𝑚𝜇2−}^𝑐𝑣∗_{𝑤𝑅∗−𝑐2},𝐺₂⁽²⁾=1𝑈₂𝑟𝑆^∗𝑘−2𝑟𝑘−2𝛽cos𝜔₀𝜏_𝑗0−𝛽𝐼^∗2𝛽cos𝜔₀𝜏_𝑗002𝑓₂⁽³⁾+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0)𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐²,𝐺₂⁽³⁾=1𝑈₂𝑟𝑆^∗𝑘𝛽𝑆^∗−2𝑟𝑘−2𝛽cos𝜔₀𝜏_𝑗−𝛽𝐼^∗𝜇₁+𝑚−𝛽𝑆^∗2𝛽cos𝜔₀𝜏_𝑗0−𝑚2𝑓₂⁽³⁾+𝑓₁⁽³⁾𝑊⁽³⁾₁₁(0), (119) where 𝑈₂=𝑟𝑆^∗𝑘𝛽𝑆^∗0−𝛽𝐼^∗𝜇₁+𝑚−𝛽𝑆^∗00−𝑚𝜇₂−𝑐𝑣^∗𝑤𝑅^∗−𝑐². (120)

It follows from the above computation and analysis that 𝑊₂₀(𝜃) and 𝑊₁₁(𝜃) can be determined based on (102) and (103).

Furthermore, we can compute 𝑔₂₁ based on (94). Hence, the following values can be computed as follows: 𝑑₁(0)=𝑖2𝜔₀𝜏_𝑗𝑔₂₀𝑔₁₁−2𝑔₁₁²−𝑔₀₂²3+𝑔₂₁2,𝛿₂𝑑=−Re₁(0)Re𝜆^′𝜏_𝑗,𝛾₂=2Re𝑑₁(0),𝑇₂=Im𝑑₁(0)+𝛿₂Im𝜆^′𝜏_𝑗𝜔₀𝜏_𝑗. (121)

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)𝛿₂ determines directions of Hopf bifurcation: if 𝛿₂>0 (𝛿₂<0), then Hopf bifurcation is supercritical (subcritical) and the bifurcating periodic solutions exist for 𝜏>𝜏_𝑗 (𝜏<𝜏_𝑗).(ii)𝛾₂ determines the stability of bifurcating periodic solutions: bifurcating periodic solutions are stable (unstable) if 𝛾₂<0 (𝛾₂>0).(iii)𝑇₂ determines the period of bifurcating periodic solutions: period increases (decreases) if 𝑇₂>0 (𝑇₂<0).

5. Numerical Simulation

In this section, some numerical simulations are provided to substantiate the theoretical results obtained in Section 4 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.

In the case of zero economic interest, model (4) without time delay takes the following form: 𝑆̇(𝑡)((𝑡(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡)).=0.11−𝑆))−𝑆=𝑆−0.1𝐼−0.1𝐼=0.1𝐼−0.1𝑅−𝐸0=𝐸)(20𝑅−1 (122)

By using Theorem 11 of this paper, it can be shown that the model (4) without time delay has a singularity induced bifurcation around the interior equilibrium 𝑃^∗(0.2,0.08,0.05,0.06), and local stability switch occurs as 𝑣 increases through 0.

Based on the analysis in Section 4.1.2 and Theorem 14 of this paper, a state feedback controller 𝑢(𝑡)=𝑙(𝐸(𝑡)−0.06) can be applied to stabilize model system (122) around 𝑃^∗, and then the model system (122) with the state feedback controller takes the following form: 𝑆̇(𝑡)((𝑡(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))(𝐸(𝑡)).=0.11−𝑆))−𝑆=𝑆−0.1𝐼−0.1𝐼=0.1𝐼−0.1𝑅−𝐸0=𝐸)(20𝑅−1+𝑙−0.06 (123)

By using Theorem 14 of this paper, if the feedback gain 𝑙 satisfies 𝑙>2.2429, then model system (122) is stable around 𝑃^∗ 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 𝑃^∗₁ and 𝑃^∗₂) when 0<𝑣<0.00668. In the following part, we focus on the case of 0<𝑣<0.00668, and the economic interest is set as 𝑣^∗=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: 𝑃^∗₁(0.2,0.08,0.0557,0.0439) and 𝑃^∗₂(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 𝑃^∗₁ and 𝑃^∗₂: 𝑆̇(𝑡)((𝑡(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))=0.11−𝑆))−𝑆=𝑆−0.1𝐼−0.1𝐼=0.1𝐼−0.1𝑅−𝐸0=𝐸)(20𝑅−1−0.005. (124)

Based on the analysis in Section 4.1.2 of this paper, state feedback controllers 𝑢(𝑡)=𝑙(̃𝐸(𝑡)−0.0439) and 𝑢(𝑡)=𝑙(̃𝐸(𝑡)−0.0115) can be applied to stabilize the model system (124) around 𝑃^∗₁ and 𝑃^∗₂, respectively. The model system (124) with respective state feedback controller takes the following form: 𝑆̇(𝑡)((𝑡(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))𝑙̃(𝐸(𝑡)),=0.11−𝑆))−𝑆=𝑆−0.1𝐼−0.1𝐼=0.1𝐼−0.1𝑅−𝐸0=𝐸)(20𝑅−1−0.005+−0.0439𝑆̇(𝑡)((𝑡(𝑡)𝐼(𝑡),𝐼̇(𝑡)(𝑡)𝐼(𝑡)(𝑡)(𝑡),𝑅̇(𝑡)(𝑡)(𝑡)(𝑡)𝑅(𝑡),(𝑡(𝑡))𝑙̃(𝐸(𝑡)).=0.11−𝑆))−𝑆=𝑆−0.1𝐼−0.1𝐼=0.1𝐼−0.1𝑅−𝐸0=𝐸)(20𝑅−1−0.005+−0.0115 (125)(126)

By using Theorem 15 of this paper, if the feedback gain 𝑙̃ satisfies 𝑙̃>3.753, then model system (125) is stable around 𝑃^∗₁ and model system (126) is stable around 𝑃^∗₂. 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.

Based on the analysis and inequality (55) in Section 4.1.2, there are two interior equilibria (denoted by 𝑃^∗₁ and 𝑃^∗₂) when 0<𝑣<0.01831. In the following part, we focus on the case of 0<𝑣<0.01831, and the economic interest is set as 𝑣^∗=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: 𝑃^∗₁(0.2,0.08,0.0907,0.0147) and 𝑃^∗₂(0.2,0.08,0.0573,0.0821). Furthermore, it can be computed that 𝑘𝛽>3(𝜇₁+𝑚). Based on the analysis in Section 4.2.1, it satisfies the assumption that (58) has a positive root, and then the corresponding 𝜏^∗₀=2.7814 can be calculated by solving (64). It follows from (121) that 𝛿₂=1.2963>0, 𝛾₂=−0.3012<0, and 𝑇₂=1.2467>0. Consequently, the interior equilibrium 𝑃^∗₁ remains stable for 𝜏<𝜏^∗₀, 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 𝑃^∗₁ are plotted; some symmetric results about 𝑃^∗₂ can be also obtained. As 𝜏 increases through 𝜏^∗₀, a periodic solution caused by the phenomenon of Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from the interior equilibrium 𝑃^∗₁. Since 𝛿₂>0 and 𝛾₂<0, the Hopf bifurcation is supercritical, the directions of the Hopf bifurcation is 𝜏>𝜏^∗₀, and these bifurcating periodic solutions from the interior equilibrium 𝑃^∗₁ at 𝜏^∗₀ are stable. Dynamical responses of model (4) with 𝜏=3>𝜏^∗₀ are plotted in Figure 5. Figures 6 and 7 show a limit cycle corresponding to the periodic solution in Figure 5 in the S-I-R and S-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 Figure 5 in the S-I-R space.

Figure 7 

Limit cycle corresponding to the periodic solution in Figure 5 in the S-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 Theorem 17. 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 𝜏>𝜏^∗₀, 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 𝜏₀=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 𝜏^∗₀=2.7814 in the case of positive economic interest. It is obvious that 𝜏^∗₀>𝜏₀, 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.

Conflict of Interests

All authors of this paper declare that there is no conflict of interests regarding the publication of this paper. They have no proprietary, financial, professional, or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, this paper.

Acknowledgments

This work is supported by National Natural Science Foundation of China, Grant no. 61104003, Grant no. 61273008, and Grant no. 61104093; Research Foundation for Doctoral Program of Higher Education of Education Ministry, Grant no. 20110042120016; Hebei Province Natural Science Foundation, Grant no. F2011501023; Fundamental Research Funds for the Central Universities, Grant no. N120423009; and Research Foundation for Science and Technology Pillar Program of Northeastern University at Qinhuangdao, Grant no. XNK201301. This work is supported by State Key Laboratory of Integrated Automation of Process Industry, Northeastern University and supported by Hong Kong Admission Scheme for Mainland Talents and Professionals, Hong Kong Special Administrative Region.