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

Abstract

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

1. Introduction

In the natural world, prey predator ecosystem with infectious disease is a well-studied ecological phenomenon. Generally, species within prey predator ecosystem does not exist alone. While species spreads the disease, it also competes with the other species for space or food, or is predated by other species. Kermack-Mc kendrick made the pioneering work on (susceptible-infective-removal-susceptible) SIRS [1]. Along with the line of this research, mathematical epidemiology has become an important subject of research in recent decades, where the evolution of a disease which gets transmitted upon contact is described. In the work done by Hadeler and Freedman [2], a prey predator model where the prey is infected by a parasite and in turn infects the predator with the parasite. Thereafter, plenty of papers were made to investigate prey predator system with infection in the prey. Arino et al. [3, 4] observed that the introduction of an infected population in the classical ratio-dependent predator-prey model may act as a biological control to save the population from extinction. Hethcote et al. [5, 6] considered the effects of disease-induced mortality or disease-reduced reproduction in regulating natural populations, decreasing their population sizes, reducing their natural fluctuations, or causing destabilization of equilibria into oscillations of the population states. In the work done by Haque and Venturino [7], the mixed terms of the underlying demographic model describing the mass action law of populations affect each other, and Holling-type functions are taken into account [8]. In the above mentioned work, by performing dynamic analysis with the help of linear approximation to model systems that consist of nonlinear equations, it concludes that it is of importance from ecological and mathematical points of view to study ecological systems subject to epidemiological factors.

Recently, an eco-epidemiological model consisting of susceptible prey, infected prey, and predator is proposed and analyzed in [9]. It is assumed that the predation functional response (trophic function) allows Holling type II (Michaelis-Menten kinetics). From the above assumptions model system can be expressed as follows: d𝑠(𝑡)(𝑡)𝑠(𝑡)(𝑡)d𝑡=𝑟𝑠1−+𝑖𝑘−𝑎₁𝑠(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)𝑖(𝑡),(𝑡)−𝑐𝑠d𝑖(𝑡)𝑖(𝑡)−𝑎d𝑡=𝑐𝑠²𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡),(𝑡)−𝑒𝑖d𝑝=𝑑d𝑡¹𝑎₁𝑠(𝑡)+𝑑₂𝑎₂𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)−𝑓𝑝(𝑡), (1) where 𝑠(𝑡) is the susceptible prey population densities per hectare and 𝑖(𝑡) is the infected prey population densities per hectare, so that the total prey population densities at any time 𝑡 per hectare is 𝑠(𝑡)+𝑖(𝑡), and it is assumed that the infected prey population contributes with the susceptible prey population to the prey population growth towards the environment carrying capacity. Consequently, the susceptible prey population grows according to a logistic law with carrying capacity 𝑘>0 and intrinsic growth rate constant 𝑟>0 involving the whole prey population (susceptible prey and infected prey population), which is best regarded as a purely descriptive equation. The transmission rate among the susceptible prey population and the infected prey population follows the simple law of mass action 𝑐𝑠(𝑡)𝑖(𝑡) for all 𝑡>0, where 𝑐>0 denotes the transmission coefficient. The infectious disease spreads among the prey population only and that disease is not genetically inherited. The infected prey population does not recover or become immune. 𝑝(𝑡) is the predator population densities per hectare which is assumed to feed on both the susceptible and the infected prey population. 𝑒 denotes the natural death rate of the infected prey populations; 𝑓 is the death rate of the predator population; 𝑑₁ and 𝑑₂ which are within the interval (0,1) are the conversion rate of the susceptible and the infected prey population to the predator population, respectively. 𝑎₁ and 𝑎₂ are the searching efficiency constants or equivalently the predation rates on the susceptible and the infected prey populations, respectively. 𝑏₁ and 𝑏₂ are the positive constants.

It is well known that harvesting has a strong impact on the dynamic evolution of a population. Evidence shows that many species have already become extinct and many others are at the verge of extinction due to several natural or man-made reasons like over exploitation, indiscriminate harvesting and mismanagement of natural resources, and so forth. The severity of this impact that may range from rapid depletion to complete preservation of a population depends on the implemented harvesting agency [10, 11]. Due to practical and economic utilization, biological resources in the prey predator system are extensively harvested nowadays. Furthermore, exploitation of biological resources has been increased by people’s multifarious material needs, which attracts a global concern to protect the limited biological resources. Therefore, regulation of exploitation of biological resources has become a problem of major concern in view of dwindling resource stocks and the deteriorating environment. It is necessary to establish a constructive management of commercial exploitation of the biological resources. The techniques and issues associated with bioeconomic exploitation have been discussed in details by Clark [12]. Taxation, license fees, lease of property rights, seasonal harvesting, and so forth are usually considered as possible governing instruments in regulation for harvesting as well as regulatory mechanisms to keep the damage to the ecosystem minimal. Out of these regulating options, taxation is considered to be superior because of its economic flexibility [13], and economists are particularly attracted to taxation because a competitive system can be better maintained under taxation rather than other regulatory methods. There has been considerable interest in the modeling of harvesting of biological resources in recent decades, where harvest effort is considered to be a dynamic variable, and optimal harvesting policies with taxation are discussed. Ganguly and Chaudhuri [14], Krishna et al. [15], and Dubey et al. [16, 17] investigated the optimal harvesting of a class of models of a single fishery species with taxation as a control. Kar [18, 19] studied the optimal tax policies for harvesting of the prey predator system.

Recently, stimulated by new analytical and numerical findings, scholars have found interest to discuss the applications of optimal control theory to eco-epidemiological model, which can be found in [20–22]. Mathematical modelling of infectious disease with application of optimal control theory is investigated in [23–25], where vaccination cost is considered to discuss associated optimal problems. It should be noted that harvesting issues are also considered in eco-epidemiological model [26, 27], where dynamical behaviour and impulsive harvesting control are investigated. However, to authors’ best knowledge, optimal control strategy associated with maximum net economic revenue of harvesting and vaccination cost in eco-epidemiological model has not been investigated in related research.

In this paper, a harvested eco-epidemiological prey predator system is established in the second section. We extend the work done in [9] by incorporating the harvest effort on predator population. Based on superiority of taxation control in regulating harvesting and effectiveness of vaccination in protecting population from infectious disease introduced above, vaccination is introduced to control number of susceptible prey population, and taxation is chosen to regulate exploitation of biological resource within the harvested eco-epidemiological prey predator system. In the third section, nonnegativity and boundedness of solutions of the proposed model are discussed. In the fourth section, global stability analysis of model system around disease-free equilibrium is studied based on robust Bendixson criterion. Furthermore, by taking taxation and vaccination control as time dependent, an optimal strategy associated with net harvesting revenue and vaccination cost is investigated by using Pontryagin’s maximum principle in the fifth section. It is assumed that vaccination cost is considered and exploitation of the predator species is regulated by an authority by imposing taxation per unit biomass of the harvested predator, and the imposition of taxation acts as a deterrent to control harvesting. We aim to find an optimal strategy which guarantees a lasting exploitation of the biological resource and maximize the total discounted net economic revenue to society as well as protect prey population from infectious disease. Numerical simulations are provided to support the analytical findings obtained in this paper. Finally, this paper ends with a conclusion.

2. Model Formulation

In 1954, Gordon [28] proposed the economic theory of a common property resource, which studies the effect of harvest effort on ecosystem from an economic perspective. It is well known that the harvested predator is usually sold as commodity in the market in order to achieve economic interest, which implies that the unit price of harvested population is influenced by the fluctuation of supply and demand for harvested population in the market, and the unit cost of harvest effort depends on fecundity of harvested population [12, 28]. Consequently, the unit price and unit cost are not always constant, but time varying. Assuming that there is always constant demand for harvested predator, then 𝑣₁(𝑡) the unit price of harvested predator and 𝑣₂(𝑡) the unit cost of harvest effort can be expressed as 𝑣₁(𝑡)=ℎ−𝜎𝛽+𝐸(𝑡)𝑝(𝑡),𝑣₂(𝑡)=𝑔𝑝(𝑡), (2) respectively, where 𝐸(𝑡) represents the harvest effort on predator population 𝑝(𝑡) at any time 𝑡>0, ℎ, 𝜎, 𝛽 and 𝑔 are all positive constants. To conserve the population in the harvested eco-epidemiological prey predator system, the regulatory agency imposes a tax 𝜎>0 per unit biomass of the predator, and 𝜎<0 denotes the subsidies given to the harvest effort. 𝐸(𝑡)𝑝(𝑡) denotes supply amount of harvested predator into the market.

By simple computation, it can be found that the unit price will decrease if the supply of harvested predator is larger than the demand for harvested predator; on the other hand, the unit price will increase if the supply of harvested predator cannot meet the demand for harvested predator, and the maximum unit price is (ℎ−𝜎)/𝛽. It is also obvious that lim_{𝑝(𝑡)→0}𝑣₂(𝑡)=lim_{𝑝(𝑡)→0}𝑔𝑝(𝑡)=∞,lim_{𝑝(𝑡)→∞}𝑣₂(𝑡)=lim_{𝑝(𝑡)→∞}𝑔𝑝(𝑡)=0, (3) which imply that the unit cost of harvest effort is inversely proportional to the population density of harvested predator. More richer the population density of harvested predator becomes, more lower the unit cost of harvest effort is; on the other hand, the unit cost of harvest effort will increase if the population density of harvested predator becomes relatively rare.

According to the above analysis, a harvested eco-epidemiological prey predator model with vaccination and taxation is proposed based on (1) and (2), which can be governed by the following differential equations: d𝑠(𝑡)(𝑡)𝑠(𝑡)(𝑡)d𝑡=𝑟𝑠1−+𝑖𝑘−𝑎₁𝑠(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)𝑖(𝑡)(𝑡),(𝑡)−𝑐𝑠−𝑢𝑠d𝑖(𝑡)𝑖(𝑡)−𝑎d𝑡=𝑐𝑠²𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡),(𝑡)−𝑒𝑖d𝑝=𝑑d𝑡¹𝑎₁𝑠(𝑡)+𝑑₂𝑎₂𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)(𝑡)𝑝(𝑡),(𝑡)−𝑓𝑝−𝐸d𝐸(𝑡)𝑝(𝑡)d𝑡=𝛼𝐸ℎ−𝜎(𝑡)𝑝(𝑡)−𝑔𝛽+𝐸𝑝(𝑡). (4) The initial conditions of model system (4) are given as follows: 𝑠(0)≥0,𝑖(0)≥0,𝑝(0)≥0,𝐸(0)≥0, (5) where 𝛼 describes stiffness parameter measuring the strength of reaction of harvest effort and susceptible prey population is vaccinated at vaccination rate 𝑢. Other parameters that share the same interpretations are mentioned in (1) and (2).

Remark 1. It should be noted that it is assumed that harvesting issue in model system (4) does not affect the growth of prey population (susceptible prey and infected prey) directly. Effect of harvest effort and harvesting regulation on prey population is biologically functioned through the predator population.

Remark 2. It is further assumed that vaccination rate 𝑢 such that constraint 0≤𝑢<1 for any time, which vividly reflects the fact that vaccination on the entire susceptible prey population is impossible.

3. Nonnegativity and Boundedness of Solution

Theorem 3. All solutions of model system (4) with initial conditions are nonnegative.

Proof. The model system (4) can be interpreted as the matrix form: d𝑋(𝑡)(𝑋(𝑡,d𝑡=𝐻)) (6) where 𝑋(𝑡)=(𝑠(𝑡),𝑖(𝑡),𝑝(𝑡),𝐸(𝑡))^𝑇∈ℝ⁴ and 𝐻(𝑋(𝑡)) is given as follows:𝐻(𝑋(𝑡=𝐻))¹(𝑋(𝑡𝐻))₂(𝑋(𝑡𝐻))₃(𝑋(𝑡𝐻))₄(𝑋(𝑡=(𝑡)𝑠(𝑡)(𝑡)))𝑟𝑠1−+𝑖𝑘−𝑎₁𝑠(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)𝑖(𝑡)(𝑡)(𝑡)𝑖(𝑡)−𝑎−𝑐𝑠−𝑢𝑠𝑐𝑠²𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)𝑑−𝑒𝑖₁𝑎₁𝑠(𝑡)+𝑑₂𝑎₂𝑖(𝑡)𝑝(𝑡)1+𝑏₁𝑠(𝑡)+𝑏₂𝑖(𝑡)(𝑡)(𝑡)𝑝(𝑡)(𝑡)𝑝(𝑡)−𝑓𝑝−𝐸𝛼𝐸ℎ−𝜎(𝑡)𝑝(𝑡)−𝑔𝛽+𝐸𝑝(𝑡). (7)Let ℝ⁴₊=[0,∞)⁴ be the nonnegative octant in ℝ⁴, then 𝐿∶ℝ₊⁴⁺¹→ℝ⁴ is locally Lipschitz and satisfy the condition 𝐻_𝑖(𝑋(𝑡))|_𝑋∈ℝ⁴≥0,𝑖=1,2,3,4.
Due to lemma in [29] and [30, Theorem A.4], any solution of the model system (4) with nonnegative initial conditions exist uniquely, and each component of the solution remains within the interval [0,𝐴₀) for some 𝐴₀>0. Standard and simple arguments show that solutions of model system (4) always exist and stay nonnegative. Hence, this completes the nonnegativity of the solutions of model system (4).

Theorem 4. If ℎ−𝜎>0 and 𝑑₁≥𝑑₂, then all solutions of model system (4) are bounded within a region Ω for all 𝑡>0, 𝑠(𝑡)(𝑡)(𝑡)(𝑡)∣𝑠(𝑡)(𝑡)≤Ω=,𝑖,𝑝,𝐸≤𝑘,𝐸ℎ−𝜎𝑔,𝑤₁(𝑡)≤𝑘(𝑟+𝑒)²4𝑒𝑟,𝑤₂(𝑡)≤𝑘𝑑₁𝑟4𝑚₁. (8)
On the other hand, if ℎ−𝜎>0 and 𝑑₁<𝑑₂, then all solutions of model system (4) are bounded within a region Ω for all 𝑡>0, 𝑠(𝑡)(𝑡)(𝑡)(𝑡)∣𝑠(𝑡)(𝑡)≤Ω=,𝑖,𝑝,𝐸≤𝑘,𝐸ℎ−𝜎𝑔,𝑤₁(𝑡)≤𝑘(𝑟+𝑒)²,𝑤4𝑒𝑟₂(𝑡)≤𝑘𝑒𝑑₁𝑟²+𝑐𝑘²𝑚₁𝑑₂−𝑑₁(𝑟+𝑒)²4𝑒𝑟𝑚₁, (9) where 𝑤₁(𝑡)=𝑠(𝑡)+𝑖(𝑡), 𝑤₂(𝑡)=𝑑₁𝑠(𝑡)+𝑑₂𝑖(𝑡)+𝑝(𝑡), and 𝑚₁=min{𝑒,𝑢,𝑓}.

Proof. Firstly, we consider the first equation of model system (4). We have d𝑠(𝑡)(𝑡)𝑠(𝑡)d𝑡≤𝑟𝑠1−𝑘, (10) which implies that lim_𝑡→∞sup𝑠(𝑡)≤𝑘. (11)
Now, defining the function 𝑤₁(𝑡)=𝑠(𝑡)+𝑖(𝑡) by calculating the time derivative of 𝑤₁(𝑡) along the solutions of model system (4), we get d𝑤₁(𝑡)≤()𝑠(𝑡)−d𝑡𝑟+𝑒𝑟𝑠²(𝑡)𝑘−𝑒𝑤₁(𝑡). (12) By utilizing the standard comparison principle, it follows that lim_𝑡→∞sup𝑤₁(𝑡)≤𝑘(𝑟+𝑒)².4𝑒𝑟 (13)
Furthermore, by using nonnegativity of solution of model system (4), if ℎ−𝜎>0, then it derives that d𝐸(𝑡)()(𝑡),d𝑡≤𝛼ℎ−𝜎−𝛼𝑔𝐸 (14) which implies that lim_𝑡→∞sup𝐸(𝑡)≤ℎ−𝜎𝑔. (15)
Let 𝑤₂(𝑡)=𝑑₁𝑠(𝑡)+𝑑₂𝑖(𝑡)+𝑝(𝑡), and by calculating the time derivative of 𝑤₂(𝑡) along the solutions of model system (4), it can be obtained that d𝑤₂(𝑡)d𝑡≤𝑟𝑑₁𝑠(𝑡)𝑠(𝑡)1−𝑘+𝑐𝑑₂−𝑑₁𝑠(𝑡)𝑖(𝑡)−𝑚₁𝑑₁𝑠(𝑡)+𝑑₂𝑖(𝑡)+𝑝(𝑡)=𝑟𝑑₁𝑠(𝑡)𝑠(𝑡)1−𝑘+𝑐𝑑₂−𝑑₁𝑠(𝑡)𝑖(𝑡)−𝑚₁𝑤₂(𝑡). (16)
Hence, it derives that lim_𝑡→∞sup𝑤₂(𝑡)=𝑘𝑑₁𝑟4𝑚₁,𝑑₁≥𝑑₂,𝑘𝑒𝑑₁𝑟²+𝑐𝑚₁𝑘²𝑑₂−𝑑₁(𝑟+𝑒)²4𝑒𝑟𝑚₁,𝑑₂>𝑑₁, (17) where 𝑚₁=min{𝑒,𝑢,𝑓}.
Consequently, If ℎ−𝜎>0 and 𝑑₁≥𝑑₂, then all solutions of model system (4) are bounded within a region Ω for all 𝑡>0, 𝑠(𝑡)(𝑡)(𝑡)(𝑡)∣𝑠(𝑡)(𝑡)≤Ω=,𝑖,𝑝,𝐸≤𝑘,𝐸ℎ−𝜎𝑔,𝑤₁(𝑡)≤𝑘(𝑟+𝑒)²4𝑒𝑟,𝑤₂(𝑡)≤𝑘𝑑₁𝑟4𝑚₁. (18)
On the other hand, if ℎ−𝜎>0 and 𝑑₁<𝑑₂, then all solutions of model system (4) are bounded within a region Ω for all 𝑡>0, 𝑠(𝑡)(𝑡)(𝑡)(𝑡)∣𝑠(𝑡)(𝑡)≤Ω=,𝑖,𝑝,𝐸≤𝑘,𝐸ℎ−𝜎𝑔,𝑤₁(𝑡)≤𝑘(𝑟+𝑒)²,𝑤4𝑒𝑟₂(𝑡)≤𝑘𝑒𝑑₁𝑟²+𝑐𝑘²𝑚₁𝑑₂−𝑑₁(𝑟+𝑒)²4𝑒𝑟𝑚₁. (19)

Remark 5. Since the components (𝑠(𝑡),  𝑖(𝑡),  and  𝑝(𝑡)) of solution of model system (4) represent the population in the harvested eco-epidemiological prey predator system, boundedness reveals a natural restriction to growth as a consequence of limited resources. Especially, boundedness of 𝑖(𝑡) epidemiologically interprets infectious disease in prey population, it can not spread without limitation, and prey population may not die out due to rampant spread of infectious disease. Furthermore, with the purpose of maintaining the sustainable development of prey predator system, the harvesting can not increase without any restriction, which is reflected by the boundedness of 𝐸(𝑡). It is of inspiration for people to regulate the harvesting effort by means of economic instrument.

4. Global Stability Analysis around Disease-Free Equilibrium

In this section, global stability analysis around disease-free equilibrium is discussed. It is utilized to discuss stability mechanism of model system due to variation of taxation and vaccination control. It should be noted that we only concentrate on the dynamical behavior and stability analysis around disease-free equilibrium, since biological interpretation of disease-free equilibrium implies that susceptible prey, predator population all exist under infectious disease-free environment, only susceptible prey population is predated and harvest effort on predator population exists, which are relevant to our main study in this paper.

For simplicity, model system (4) is nondimensionalized with the following scaling: 𝑆=𝑠/𝑘,  𝐼=𝑖/𝑘,  𝑃=𝑎₁𝑝/𝑟,  𝑌=𝐸/𝑟,  𝑇=𝑟𝑡. By using these quantities, model system (4) can be transformed into a dimensionless form as follows: d𝑆(𝑇)(𝑇(𝑇)(𝑇−𝑆(𝑇)𝑃(𝑇)d𝑇=𝑆)(1−𝑆−𝐼))1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑐₁𝑆(𝑇)𝐼(𝑇)−𝑢₁𝑆(𝑇),d𝐼(𝑇)d𝑇=𝑐₁𝑆(𝑇)𝐼(𝑇)−𝑎₃𝐼(𝑇)𝑃(𝑇)1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑒₁𝐼(𝑇),d𝑃(𝑇)=𝑑d𝑇³𝑆(𝑇)+𝑑₄𝐼(𝑇)𝑃(𝑇)1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑓₁𝑃(𝑇)(𝑇)𝑌(𝑇),(𝑇)−𝑃d𝑌d𝑇=𝛼₁𝑃(𝑇)𝑌(𝑇)ℎ₁−𝜎₁𝛽₁(𝑇)𝑌(𝑇)−𝑔+𝑃𝑃(𝑇¹), (20) where 𝑐₁=𝑐𝑘/𝑟, 𝑒₁=𝑒/𝑟, 𝑎₃=𝑎₂/𝑎₁, 𝑏₃=𝑏₁𝑘, 𝑏₄=𝑏₂𝑘, 𝑓₁=𝑓/𝑟, 𝑑₃=(𝑑₁𝑎₁𝑘)/𝑟, 𝑑₄=𝑑₂𝑎₂𝑘/𝑟, 𝛼₁=𝛼/𝑎₁, 𝛽₁=(𝑎₁𝛽)/𝑟², ℎ₁−𝜎₁=(𝑎₁(ℎ−𝜎))/𝑟², 𝑢₁=𝑢/𝑟, and 𝑔₁=(𝑎₁𝑔)/𝑟.

By simple computation, disease-free equilibrium 𝑀^∗(𝑆^∗,0,𝑃^∗,𝑌^∗) can be obtained as follows.

𝑃^∗=(1−𝑢₁−𝑆^∗)(1+𝑏₃𝑆^∗), 𝑌^∗=((𝑑₃−𝑏₃𝑓₁)𝑆^∗−𝑓₁)/(1+𝑏₃𝑆^∗), and 𝑆^∗ satisfies the following equation: 𝑔₁𝑑₃−𝑏₃𝑓₁−𝑏₃ℎ₁−𝜎₁𝑆^∗2+1−𝑢₁𝑏₃×ℎ−1₁−𝜎₁+𝑓₁𝑔₁𝑆^∗−1−𝑢₁𝑑₃𝑔₁𝑆^∗+ℎ₁−𝜎₁+𝑔₁𝑓₁1−𝑢₁−𝑔₁𝛽₁=0. (21)

Conditions for existence of unique positive real root of (21) are as follows: 𝑔₁𝛽₁−𝑓₁≤𝜎₁≤ℎ₁+𝑔₁𝑓₁−𝑔₁𝑑₃𝑏₃,1−𝑏₃𝛽₁+𝑓₁+𝑑₃≤𝑢₁<𝑟. (22)

Remark 6. The above inequalities associated with taxation and vaccination control in (22) also guarantee existence and uniqueness of disease-free equilibrium 𝑀^∗.

Theorem 7. If the following inequalities hold 𝑎₃+𝑒₁>𝑎₃+𝑐₁𝑆^∗,𝑆^∗−1𝑏₃ℎ₁−𝜎₁+2𝑔₁𝑏₃𝑓₁−𝑑₃𝑢>0,₁>𝑏₃−1−2𝑑₃𝑏₃, (23) then model system (20) is locally stable around disease-free equilibrium 𝑀^∗.

Proof. For disease-free equilibrium 𝑀^∗(𝑆^∗,0,𝑃^∗,𝑌^∗), the variational matrix of model system (20) around 𝑀^∗ takes the following form: 𝑉𝑆^∗,0,𝑃^∗,𝑌^∗=𝑉₁₁𝑉₁₂𝑉₁₃𝑉₁₄𝑉₂₁𝑉₂₂𝑉₂₃𝑉₂₄𝑉₃₁𝑉₃₂𝑉₃₃𝑉₃₄𝑉₄₁𝑉₄₂𝑉₄₃𝑉₄₄, (24) where 𝑉₁₁=𝑆^∗1−𝑢₁𝑏₃−1−2𝑏₃𝑆^∗1+𝑏₃𝑆^∗,𝑉₁₂=𝑆^∗𝑏₄−𝑐₁−1−𝑏₄+𝑏₃+𝑏₃𝑐₁𝑆^∗1+𝑏₃𝑆^∗,𝑉₁₃𝑆=−^∗1+𝑏₃𝑆^∗,𝑉₂₂=𝑐₁+𝑎₃𝑆^∗−𝑎₃−𝑒₁,𝑉₃₁=𝑑₃1−𝑆^∗1+𝑏₃𝑆^∗,𝑉₃₂=1−𝑆^∗𝑑₄+𝑏₃𝑑₄−𝑏₄𝑑₃𝑆^∗1+𝑏₃𝑆^∗,𝑉₃₃𝑑=−³𝑆^∗1+𝑏₃𝑆^∗,𝑉₃₄=𝑆^∗−11+𝑏₃𝑆^∗,𝑉₄₃=Δ^∗𝑑₃−𝑏₃𝑓₁𝑆^∗−𝑓₁1+𝑏₃𝑆^∗,𝑉₄₄=Δ^∗1−𝑆^∗1+𝑏₃𝑆^∗−𝛼₁𝑔₁,𝑉₁₄=𝑉₂₁=𝑉₂₃=𝑉₂₄=𝑉₄₁=𝑉₄₂=0. (25)
It follows from 𝑉(𝑆^∗,0,𝑃^∗,𝑌^∗) that the characteristic equation around 𝑀^∗ is as follows: 𝜆+𝑎₃+𝑒₁−𝑎₃+𝑐₁𝑆^∗×𝜆³+𝐴^∗₁𝜆²+𝐴^∗₂𝜆+𝐴^∗₃=0, (26)where 𝐴^∗₁=𝛼₁𝑆^∗−11+𝑏₃𝑆^∗Δ^∗+𝑔₁,Δ^∗=𝛽₁ℎ₁−𝜎₁𝛽₁+(1−𝑆^∗)𝑆^∗𝑑₃−𝑏₃𝑓₁−𝑓₁²,𝐴^∗₂=2𝑏₃𝑏₃+𝑑₃𝑆^∗3+3𝑏₃−𝑏²₃−2𝑑₃𝑆^∗2+2𝑑₃−𝑏₃+𝑢₁𝑏₃+1𝑆^∗×1+𝑏₃𝑆^∗2−1+𝛼₁Δ^∗𝑆^∗−11+𝑏₃𝑆^∗2𝑑₃−𝑏₃𝑓₁𝑆^∗−𝑓₁1+𝑏₃𝑆^∗2,𝐴^∗₃=𝑆^∗𝐴^∗₁2𝑏₃𝑏₃+𝑑₃𝑆^∗2+3𝑏₃−𝑏²₃−2𝑑₃𝑆^∗+2𝑑₃−𝑏₃+𝑢₁𝑏₃×+11+𝑏₃𝑆^∗2−1+𝛼₁𝑆^∗𝐴^∗₁𝐴^∗₁−𝑔₁2𝑏₃𝑆^∗−𝑏₃×𝑑+1₃−𝑏₃𝑓₁𝑆^∗−𝑓₁×1+𝑏₃𝑆^∗2−1. (27)
It follows from Routh-Hurwitz criteria [13] that there are four eigenvalues with negative real part provided that 𝐴^∗₁>0,𝐴^∗₃>0, and 𝐴^∗₁𝐴^∗₂>𝐴^∗₃, which derives that 𝑎₃+𝑒₁>𝑎₃+𝑐₁𝑆^∗,𝑆^∗−1𝑏₃ℎ₁−𝜎₁+2𝑔₁𝑏₃𝑓₁−𝑑₃𝑢>0,₁>𝑏₃−1−2𝑑₃𝑏₃. (28)

In the following part, global stability analysis of model system (20) around disease-free equilibrium 𝑀^∗ is investigated. Firstly, some preliminaries are discussed as follows.

Lemma 8 (see [31]). The unique interior equilibrium is globally stable in the domain 𝐺 if it is locally stable and all trajectories in 𝐺 converge to the unique interior equilibrium, where domain 𝐺 satisfies the following conditions:(i)the domain 𝐺 is simply connected;(ii)there is a compact absorbing set 𝐽⊂𝐺.

Let 𝑉^|2| be the second additive compound matrix [31, 32]. In the case of 𝑛=4, 𝑉^|2| can be written based on variational matrix of model system (20) around 𝑀^∗, which is as follows:𝑉^|2|=𝑉₁₁+𝑉₂₂𝑉₂₃𝑉₂₄−𝑉₁₃−𝑉₁₄0𝑉₃₂𝑉₁₁+𝑉₃₃𝑉₃₄𝑉₁₂0−𝑉₁₄𝑉₄₂𝑉₄₃𝑉₁₁+𝑉₄₄0𝑉₁₂𝑉₁₃−𝑉₃₁𝑉₂₁0𝑉₂₂+𝑉₃₃𝑉₃₄−𝑉₂₄−𝑉₄₁0𝑉₂₁𝑉₄₃𝑉₂₂+𝑉₄₄𝑉₂₃0−𝑉₄₁𝑉₃₁−𝑉₄₂𝑉₃₂𝑉₃₃+𝑉₄₄, (29)where 𝑉_𝑖𝑗,  𝑖,𝑗=1,2,3,4, can be found in Theorem 7. Since 𝑉₁₄=0, 𝑉₂₁=0, 𝑉₂₃=0, 𝑉₂₄=0, 𝑉₄₁=0, and 𝑉₄₂=0, 𝑉^|2| can be deduced as follows:𝑉^|2|=𝑉₁₁+𝑉₂₂00−𝑉₁₃𝑉00₃₂𝑉₁₁+𝑉₃₃𝑉₃₄𝑉₁₂000𝑉₄₃𝑉₁₁+𝑉₄₄0𝑉₁₂𝑉₁₃−𝑉₃₁00𝑉₂₂+𝑉₃₃𝑉₃₄0000𝑉₄₃𝑉₂₂+𝑉₄₄000𝑉₃₁0𝑉₃₂𝑉₃₃+𝑉₄₄. (30)

Furthermore, we have 𝑁⁻¹𝑌=diag𝑆,𝑌𝑆,𝑌𝑆,𝑌𝑆,𝑁_𝐹(𝑍)=d𝑁𝑆̇d𝑍=diag𝑌−𝑆𝑌̇𝑌²,𝑆̇𝑌−𝑆𝑌̇𝑌²,𝑆̇𝑌−𝑆𝑌̇𝑌², (31) where 𝐹(𝑍(𝑇)) and 𝑍(𝑇) are defined based on model system (20), and the model system (20) can be interpreted as the matrix formd𝑍(𝑇)(𝑍(𝑇,d𝑇=𝐹)) (32) where 𝑍(𝑇)=(𝑆(𝑇),𝐼(𝑇),𝑃(𝑇),𝑌(𝑇))^𝑇 and 𝐹(𝑍(𝑇)) is given as follows:𝐹(𝑍(𝑇=𝑆(𝑇(𝑇)(𝑇−𝑆(𝑇)𝑃(𝑇))))(1−𝑆−𝐼))1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑐₁𝑆(𝑇)𝐼(𝑇)−𝑢₁𝑆(𝑇)𝑐₁𝑆(𝑇)𝐼(𝑇)−𝑎₃𝐼(𝑇)𝑃(𝑇)1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑒₁𝐼(𝑇)𝑑₃𝑆(𝑇)+𝑑₄𝐼(𝑇)𝑃(𝑇)1+𝑏₃𝑆(𝑇)+𝑏₄𝐼(𝑇)−𝑓₁𝑃(𝑇)(𝑇)𝑌(𝑇)𝛼−𝑃₁𝑃(𝑇)𝑌(𝑇ℎ)¹−𝜎₁𝛽₁(𝑇)𝑌(𝑇)−𝑔+𝑃𝑃(𝑇¹). (33)

By simple computation, it can be obtained that 𝑁_𝐹𝑁⁻¹𝑆̇=diag𝑆−𝑌̇𝑌,𝑆̇𝑆−𝑌̇𝑌,𝑆̇𝑆−𝑌̇𝑌,𝑁𝑉^|2|𝑁⁻¹=𝑉^|2|. (34)

Let 𝐷=𝑁_𝐹𝑁⁻¹+𝑁𝑉^|2|𝑁⁻¹=𝐷₁₁𝐷₁₂𝐷₁₃𝐷₂₁𝐷₂₂𝐷₂₃𝐷₃₁𝐷₃₂𝐷₃₃, (35)where𝐷₁₁=𝑆̇𝑆−𝑌̇𝑌+𝑉₁₁+𝑉₂₂,𝐷₁₂=(),𝐷0,0₁₃=−𝑉₁₃,0,0,𝐷₂₁=𝑉₃₂,0^𝑇,𝐷₂₂=𝑆̇𝑆−𝑌̇𝑌+𝑉₁₁+𝑉₃₃𝑉₃₄𝑉₄₃𝑆̇𝑆−𝑌̇𝑌+𝑉₁₁+𝑉₄₄,𝐷₂₃=𝑉₁₂000𝑉₁₂𝑉₁₃,𝐷₃₁=−𝑉₃₁,0,0^𝑇,𝐷₃₂=00000𝑉₃₁,𝐷₃₃=𝑆̇𝑆−𝑌̇𝑌+𝑉₂₂+𝑉₃₃𝑉₃₄0𝑉₄₃𝑆̇𝑆−𝑌̇𝑌+𝑉₂₂+𝑉₄₄00𝑉₃₂𝑆̇𝑆−𝑌̇𝑌+𝑉₃₃+𝑉₄₄. (36)

Consider the Lozinskii measure Γ of 𝐷 with respect to a vector norm |⋅| in ℝ⁴ (see, [32]) and 𝐼 (identity matrix), then it can be obtained that Γ(𝐷)=lim_𝜖→0+|𝐼+𝜖𝐷|−1𝜖. (37)

Hence, it follows from simple computation that Γ(𝐷)≤sup𝑓₁,𝑓₂,𝑓₃, (38) where 𝑓_𝑖=Γ₁(𝐷_𝑖𝑖)+|𝐷_𝑖𝑗| for 𝑖,𝑗=1,2,3, and 𝑖≠𝑗. |𝐷₁₂|, |𝐷₁₃|, |𝐷₂₁|, |𝐷₂₃|, |𝐷₃₁|, and |𝐷₃₂| are matrix norms with respect to the 𝐿¹ vector norm, and Γ₁ is the Lozinskii measure with respect to this vector norm.

Furthermore, let Δ=𝛽₁(ℎ₁−𝜎₁)/[𝛽₁+(1−𝑆)(𝑆(𝑑₃−𝑏₃𝑓₁)−𝑓₁)]², then it can be computed that Γ₁𝐷₁₁=𝑆̇𝑆−𝑌̇𝑌+1−𝑢₁𝑏₃𝑆−𝑆−2𝑏₃𝑆²1+𝑏₃𝑆+𝑐₁+𝑎₃𝑆−𝑎₃+𝑒₁,𝐷₁₂=0,𝐷₁₃=𝑆1+𝑏₃𝑆,Γ₁𝐷₂₂=𝑆̇𝑆−𝑌̇𝑌𝑑+max³𝑆1+𝑏₃𝑆,𝛼₁Δ(1−𝑆)1+𝑏₃𝑆+𝛼₁𝑔₁+𝑆1−𝑢₁𝑏₃−1−2𝑏₃𝑆1+𝑏₃𝑆+max(1−𝑆)1+𝑏₃𝑆,𝛼₁Δ𝑑₃−𝑏₃𝑓₁𝑆−𝑓₁1+𝑏₃𝑆,𝐷₂₁=𝑑₃𝑆1+𝑏₃𝑆,𝐷₂₃𝑆𝑏=max⁴−1−𝑐₁−𝑏₄+𝑏₃+𝑏₃𝑐₁𝑆1+𝑏₃𝑆,𝑆1+𝑏₃𝑆,𝐷₃₁=𝐷₃₂=𝑑₃(1−𝑆)1+𝑏₃𝑆,Γ₁𝐷₃₃=𝑆̇𝑆−𝑌̇𝑌+𝛼₁Δ+1(1−𝑆)1+𝑏₃𝑆+𝑑₃𝑆1+𝑏₃𝑆+𝑐₁+𝑎₃𝑆−𝑎₃+𝑒₁+𝛼₁maxΔ𝑆𝑑₃−𝑏₃𝑓₁−𝑓₁,𝑔₁. (39)

According to Theorem 3 and definition of 𝑓_𝑖,𝑖=1,2,3, it follows that 𝑓₁≤𝑆̇𝑆−𝑌̇𝑌+2+𝑐₁+𝑒₁,𝑓₂≤𝑆̇𝑆−𝑌̇𝑌𝑏+2+max1,⁴−1−𝑐₁²4𝑏₄+𝑏₃1+𝑐₁+max1+𝑏₃²4𝑏₃,16𝛼₁𝛽₁ℎ₁−𝜎₁𝑑₃−𝑏₃𝑓₁²×𝑑₃+𝑓₁𝑏₃×+11+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹+max1,16𝛼₁𝛽₁ℎ₁−𝜎₁1+𝑏₃𝑑₃−𝑏₃𝑓₁²×1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹+𝛼₁𝑔₁,𝑓₃≤𝑆̇𝑆−𝑌̇𝑌+1+𝑐₁+𝑒₁+1+16𝛼₁𝛽₁ℎ₁−𝜎₁𝑑₃−𝑏₃𝑓₁²1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²×𝑏1+³+1²4𝑏₃+max𝛼₁𝑔₁,16𝛼₁𝛽₁ℎ₁−𝜎₁×𝑓₁𝑑₃+𝑏₃+1𝑑₃−𝑏₃𝑓₁²×1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹. (40)

Lemma 9. For element 𝑃(𝑇) of the disease-free equilibrium 𝑀^∗, some inequalities about liminf_𝑇→∞𝑃(𝑇) hold as follows: (a) if 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)≤0, then liminf_𝑇→∞𝑃(𝑇)≥𝑑₃𝑚₂+1−𝑢₁²4𝑚₂−𝑑²₃+𝑑₄𝑚₃𝑑₃; (41)(b) if 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)>0, then(i)liminf_𝑇→∞𝑃(𝑇)≥[𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)−𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄)]²/4𝑑₃𝑚₂𝑎²₃  −  (𝑑²₃+𝑑₄𝑚₃)/𝑑₃, provided that 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)≥𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄),(ii)liminf_𝑇→∞𝑃(𝑇)≥(𝑎₃𝑑₃(𝑚₂−𝑚₃)−𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄))/𝑎₃𝑚₂−(𝑑²₃+𝑑₄𝑚₃)/𝑑₃, provided that 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)<𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄), where 𝑚₂=max{𝑒₁,𝑓₁+𝑟(ℎ₁−𝜎₁)/𝑔₁} and 𝑚₃=𝑘(𝑒₁+1)²/4𝑒₁.

Proof. Based on Theorem 4 and nondimensionalized scaling of model system (20) mentioned in Section 4, it derives that 𝑌(𝑇)≤𝑟ℎ₁−𝜎₁𝑔₁(𝑇)≤𝑘,𝐼1+𝑒₁²4𝑒₁,𝑆(𝑇)≤1. (42)
According to model system (20), let 𝑊(𝑇)=𝑑₃𝑆(𝑇)+(𝑑₄/𝑎₃)𝐼(𝑇)+𝑃(𝑇). By calculating the time derivative of 𝑊(𝑇) along the solutions of model system (20), it can be obtained that d𝑊d𝑇=𝑑₃𝑆(1−𝑆−𝐼)−𝑐₁𝑆𝐼𝑑₃−𝑑₄𝑎₃−𝑑₄𝑒₁𝐼𝑎₃−𝑓₁𝑃−𝑃𝑌−𝑑₃𝑢₁𝑆=𝑑₃𝑆1−𝑢₁−𝑆−𝑑₃+𝑐₁𝑑₃−𝑑₄𝑎₃−𝑑𝑆𝐼₄𝑒₁𝐼𝑎₃−𝑓₁𝑃−𝑃𝑌≥𝑑₃1−𝑢₁+𝑚₂𝑆−𝑑₃𝑆²−𝑑₃+𝑐₁𝑑₃−𝑑₄𝑎₃𝑆𝐼−𝑚₂𝑊, (43) where 𝑚₂=max{𝑒₁,𝑓₁+𝑟(ℎ₁−𝜎₁)/𝑔₁}.
If 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)≤0, then d𝑊d𝑇≥𝑑₃1−𝑢₁+𝑚₂𝑆−𝑑₃𝑆²−𝑚₂𝑊. (44)
According to 𝑊=𝑑₃𝑆+(𝑑₄/𝑎₃)𝐼+𝑃, it follows from the standard comparison principle that liminf_𝑇→∞𝑃(𝑇)≥𝑑₃𝑚₂+1−𝑢₁²4𝑚₂−𝑑²₃+𝑑₄𝑚₃𝑑₃∶=𝑛₁. (45)
On the other hand, if 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)>0, then d𝑊≥𝑑d𝑇³1−𝑢₁+𝑚₂−𝑚₃−𝑐₁𝑚₃𝑑₃−𝑑₄𝑎₃𝑆−𝑑₃𝑆²−𝑚₂𝑊, (46) where 𝑚₃=𝑘(𝑒₁+1)²/4𝑒₁.
According to 𝑊=𝑑₃𝑆+(𝑑₄/𝑎₃)𝐼+𝑃, it follows from the standard comparison principle that liminf_𝑇→∞𝑃(𝑇)≥𝑎₃𝑑₃1−𝑢₁+𝑚₂−𝑚₃−𝑐₁𝑚₃𝑎₃𝑑₃−𝑑₄²4𝑑₃𝑚₂𝑎²₃−𝑑²₃+𝑑₄𝑚₃𝑑₃∶=𝑛₂, (47) provided that 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)≥𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄).
Furthermore, liminf_𝑇→∞𝑃(𝑇)≥𝑎₃𝑑₃𝑚₂−𝑚₃−𝑐₁𝑚₃𝑎₃𝑑₃−𝑑₄𝑎₃𝑚₂−𝑑²₃+𝑑₄𝑚₃𝑑₃∶=𝑛₃, (48) provided that 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)<𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄).

Let 𝑙₁=𝛼₁𝑔₁+1+𝑐₁+𝑒₁,𝑙₂𝑏=1+max1,⁴−1−𝑐₁²4𝑏₄+𝑏₃1+𝑐₁+max1+𝑏₃²4𝑏₃,16𝛼₁𝛽₁ℎ₁−𝜎₁𝑑₃−𝑏₃𝑓₁²×𝑑₃+𝑓₁𝑏₃×+11+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹+max1,16𝛼₁𝛽₁ℎ₁−𝜎₁1+𝑏₃𝑑₃−𝑏₃𝑓₁²×1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹+𝛼₁𝑔₁,𝑙₃=𝑐₁+𝑒₁+1+16𝛼₁𝛽₁ℎ₁−𝜎₁𝑑₃−𝑏₃𝑓₁²1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²×𝑏1+³+1²4𝑏₃+max𝛼₁𝑔₁,16𝛼₁𝛽₁ℎ₁−𝜎₁+𝑓₁𝑑₃+𝑏₃+1𝑑₃−𝑏₃𝑓₁²×1+𝑏₃𝑓₁−𝑑₃²+4𝛽₁𝑑₃−𝑏₃𝑓₁²⁻¹. (49)

Theorem 10. If the model system is locally stable around 𝑀^∗ and inequality (22) holds, then(a)if 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)≤0 and 𝛼₁𝑔₁𝑛₁(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₁(ℎ₁−𝜎₁))>𝑙₁+𝑙₂+𝑙₃, then model system (20) is globally stable around 𝑀^∗; (b)if 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)>0, then(i)model system (20) is globally stable around 𝑀^∗ provided that  𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)≥𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄),  𝛼₁𝑔₁𝑛₂(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₂(ℎ₁−𝜎₁))>𝑙₁  +  𝑙₂+𝑙₃,(ii)model system (20) is globally stable around 𝑀^∗ provided that  𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)<𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄),  𝛼₁𝑔₁𝑛₃(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₃(ℎ₁−𝜎₁))>𝑙₁  +  𝑙₂+𝑙₃, where 𝑛_𝑖,𝑖=1,2,3, and 𝑚_𝑖,𝑖=2,3, have been defined in Lemma 9. 𝑙_𝑖,𝑖=1,2,3, have been defined in (49)–(9).

Proof. According to model system (20), if inequality (22) holds, then there exists a unique disease-free equilibrium 𝑀^∗.
If 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)≤0 and 𝛼₁𝑔₁𝑛₁(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₁(ℎ₁−𝜎₁))>𝑙₁+𝑙₂+𝑙₃, then Γ(𝐷)≤sup𝑓₁,𝑓₂,𝑓₃≤𝑆̇𝑆−𝛼₁𝑔₁𝑛₁ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₁ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃, (50) which derives that 1𝑇^𝑇₀Γ(𝐷)1d𝑆≤𝑇𝑆(𝑇)log𝑆(0)−𝛼₁𝑔₁𝑛₁ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₁ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃,limsup_𝑇→+∞1𝑇^𝑇₀Γ(𝐷)𝛼d𝑆<−¹𝑔₁𝑛₁ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₁ℎ₁−𝜎₁−𝑙₁−𝑙₂−𝑙₃<0. (51)
Based on robust Bendixson criterion [31, 32], inequality (51) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ⁴, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around 𝑀^∗, then it follows from Lemma 8 that model system (20) is globally stable around 𝑀^∗.
If 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)>0, 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)≥𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄), and 𝛼₁𝑔₁𝑛₂(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₂(ℎ₁−𝜎₁))>𝑙₁+𝑙₂+𝑙₃, then Γ(𝐷)≤sup𝑓₁,𝑓₂,𝑓₃≤𝑆̇𝑆−𝛼₁𝑔₁𝑛₂ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₂ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃, (52) which derives that 1𝑇^𝑇₀Γ(𝐷)1d𝑆≤𝑇𝑆(𝑇)log𝑆(0)−𝛼₁𝑔₁𝑛₂ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₂ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃,limsup_𝑇→+∞1𝑇^𝑇₀Γ(𝐷)𝛼d𝑆<−¹𝑔₁𝑛₂ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₂ℎ₁−𝜎₁−𝑙₁−𝑙₂−𝑙₃<0. (53)
Based on robust Bendixson criterion [31, 32], inequality (53) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ⁴, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around 𝑀^∗, then it follows from Lemma 8 that model system (20) is globally stable around 𝑀^∗.
If 𝑎₃𝑑₃+𝑐₁(𝑎₃𝑑₃−𝑑₄)>0, 𝑎₃𝑑₃(1−𝑢₁+𝑚₂−𝑚₃)<𝑐₁𝑚₃(𝑎₃𝑑₃−𝑑₄), and 𝛼₁𝑔₁𝑛₃(ℎ₁−𝜎₁)/(𝛽₁𝑔₁+𝑟𝑛₃(ℎ₁−𝜎₁))>𝑙₁+𝑙₂+𝑙₃, then Γ(𝐷)≤sup𝑓₁,𝑓₂,𝑓₃≤𝑆̇𝑆−𝛼₁𝑔₁𝑛₃ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₃ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃. (54) It derives that 1𝑇^𝑇₀Γ(𝐷)1d𝑆≤𝑇𝑆(𝑇)log𝑆(0)−𝛼₁𝑔₁𝑛₃ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₃ℎ₁−𝜎₁+𝑙₁+𝑙₂+𝑙₃,limsup_𝑇→+∞1𝑇^𝑇₀Γ(𝐷)𝛼d𝑆<−¹𝑔₁𝑛₃ℎ₁−𝜎₁𝛽₁𝑔₁+𝑟𝑛₃ℎ₁−𝜎₁−𝑙₁−𝑙₂−𝑙₃<0. (55)
Based on robust Bendixson criterion [31, 32], inequality (55) ensures that there are no orbits (i.e., homoclinic orbits, heteroclinic cycles, and periodic orbits), which give rise to a simple closed rectifiable curve in ℝ⁴, invariant for the model system (20). Based on the above analysis, if model system (20) is locally stable around 𝑀^∗, then it follows from Lemma 8 that model system (20) is globally stable around 𝑀^∗.

5. Optimal Control of Model System

With the purpose of regulating harvesting and protecting population from infectious disease within harvested eco-epidemiological system, we design an optimal control to maximize total discounted net revenue to society by using taxation and vaccination as control instruments. The path traced out by (𝑆,𝐼,𝑃,𝑌) with optimal taxation 𝜎₁ and optimal vaccination 𝑢₁ is also investigated.

Total discounted net economic revenue to society 𝜋(𝑆,𝐼,𝑃,𝑌,𝜎₁,𝑢₁,𝑇)  = discounted economic revenue of harvesting + discounted economic revenue to regulatory agency − vaccination cost.

It follows from model system (20), that the following can be obtained: 𝜋𝑆,𝐼,𝑃,𝑌,𝜎₁,𝑢₁,𝑇=𝑒^−𝛿𝑇ℎ₁−𝜎₁𝑃𝑌𝛽₁+𝑃𝑌−𝑔₁𝑌+𝑒^−𝛿𝑇𝜎₁𝑃𝑌𝛽₁−1+𝑃𝑌2𝐵𝑢²₁=𝑒^−𝛿𝑇ℎ₁𝑃𝑌𝛽₁+𝑃𝑌−𝑔₁𝑌−12𝐵𝑢²₁, (56) where 𝛿 is the instantaneous annual rate of discount and the optimization problem is subject to the model system (20). 𝐵 is a positive weight parameter associated with the square of control variable to balance size of the terms; see [20, 21, 24, 25] and references therein. (1/2)𝐵𝑢²₁ represents vaccination cost of susceptible prey population.

Remark 11. Vaccination cost discussed in this section is a quadratic cost, which is the simplest and widest used nonlinear representation of vaccination cost; see [20, 21, 24, 25] and references therein. Actually, when nonlinear description is adopted, it may be not completely clear which nonlinear form must be appropriately chosen. In this case, choosing the simplest form compatible with mathematical requirements (existence, well-posedness of maximization problems, etc.) is a possible guideline; see [20, 21, 24] and references therein. Based on the above introduction, we consider a quadratic cost on the optimal control problem in this section. It should be noted that other nonlinear functions might provide a better description of the actual vaccination cost due to the increase of vaccination cost when most of population is already removed (vaccinated or immune); see [21, 24] and references therein.

5.1. Design of Optimal Control

Our objective is to maximize the following optimization problem: max^∞₀𝑒^−𝛿𝑇ℎ₁𝑃𝑌𝛽₁+𝑃𝑌−𝑔₁𝑌−12𝐵𝑢²₁𝑑𝑇. (57)

By using the Pontryagin’s maximum principle [12], the associated Hamiltonian function is constructed by 𝑅𝑆,𝐼,𝑃,𝑌,𝜎₁,𝑢₁,𝑇=𝑒^−𝛿𝑇ℎ₁𝑃𝑌𝛽₁+𝑃𝑌−𝑔₁𝑌−12𝐵𝑢²₁+𝜆₁(𝑇)×𝑆(1−𝑆−𝐼)−𝑆𝑃1+𝑏₃𝑆+𝑏₄𝐼−𝑐₁𝑆𝐼−𝑢₁𝑆+𝜆₂(𝑇)𝑐₁𝑎𝑆𝐼−³𝐼𝑃1+𝑏₃𝑆+𝑏₄𝐼−𝑒₁𝐼+𝜆₃(𝑇𝑑)³𝑆+𝑑₄𝐼𝑃1+𝑏₃𝑆+𝑏₄𝐼−𝑓₁𝑃−𝑃𝑌+𝜆₄(𝑇)𝛼₁ℎ𝑃𝑌¹−𝜎₁𝛽₁−𝑔+𝑃𝑌𝑃¹, (58) where 𝜆_𝑖(𝑇),𝑖=1,2,3,4, are adjoint variables. 𝜎₁ and 𝑢₁ are the control variables. According to [33], the conditions for a singular control to be optimal can be obtained by  𝜕𝑅/𝜕𝜎₁=0 and 𝜕𝑅/𝜕𝑢₁=0, from which we get 𝜆₄(𝑇)𝛼₁𝑃𝑌𝛽₁+𝑃𝑌=0,𝜆₁(𝑇)𝑆+𝐵𝑢₁=0. (59)

According to Theorem 3, it follows that 𝜆₄(𝑇)=0,𝑢₁𝜆=−¹(𝑇)𝑆𝐵. (60)

For adjoint variables 𝜆_𝑖(𝑇),𝑖=1,2,3,4, we have d𝜆₁(𝑇)d𝑇=−𝜕𝑅,𝜕𝑆d𝜆₂(𝑇)d𝑇=−𝜕𝑅,𝜕𝐼d𝜆₃(𝑇)d𝑇=−𝜕𝑅,𝜕𝑃d𝜆₄(𝑇)d𝑇=−𝜕𝑅.𝜕𝑌 (61)

By substituting 𝜆₄(𝑇)=0 into the Equation d𝜆₄(𝑇)/d𝑇=−𝜕𝑅/𝜕𝑌, it follows that 𝜆₃(𝑇)=𝑒^−𝛿𝑇ℎ₁𝛽₁𝛽₁+𝑃𝑌²−𝑔₁𝑃. (62)

By virtue of (62), it gives that d𝜆₁(𝑇)=𝑐d𝑇2𝑆+₁𝑃+1𝐼+1+𝑏₄𝐼1+𝑏₃𝑆+𝑏₄𝐼²−1+𝑢₁𝜆₁−𝐼𝑐₁+𝑎₃𝑏₃𝑃1+𝑏₃𝑆+𝑏₄𝐼²𝜆₂+𝑒^−𝛿𝑇𝑔₁−ℎ₁𝛽₁𝑃𝛽₁+𝑃𝑌²𝑏₄𝑑₃−𝑏₃𝑑₄𝐼+𝑑₃1+𝑏₃𝑆+𝑏₄𝐼². (63)

By associating with 𝑀^∗(𝑆^∗,0,𝑃^∗,𝑌^∗), the above equation can be deduced as follows: d𝜆₁(𝑇)d𝑇=𝑄₂𝜆₁−𝑑₃1−𝑆^∗𝑄₁1+𝑏₃𝑆^∗𝑒^−𝛿𝑇, (64) where 𝑄₁=ℎ₁𝛽₁/(𝛽₁+𝑃^∗𝑌^∗)²−𝑔₁/𝑃^∗ and 𝑄₂=[2𝑏₃𝑆^∗−𝑏₃+1+𝑏₃𝑢₁]𝑆^∗/(1+𝑏₃𝑆^∗).

By solving (64), it can be obtained that 𝜆₁(𝑇)=𝑑₃1−𝑆^∗𝑄₁𝛿+𝑄₂1+𝑏₃𝑆^∗𝑒^−𝛿𝑇. (65)

By virtue of (62), it gives that d𝜆₂(𝑇)d𝑇=𝑆1+𝑐₁−𝑏₄𝑃1+𝑏₃𝑆+𝑏₄𝐼²𝜆₁+𝑎₃𝑃1+𝑏₃𝑆1+𝑏₃𝑆+𝑏₄𝐼²+𝑒₁−𝑐₁𝑆𝜆₂+𝑒^−𝛿𝑇𝑔₁−ℎ₁𝛽₁𝑃𝛽₁+𝑃𝑌²𝑑₄+𝑏₃𝑑₄−𝑏₄𝑑₃𝑆1+𝑏₃𝑆+𝑏₄𝐼². (66)

By associating with 𝑀^∗(𝑆^∗,0,𝑃^∗,𝑌^∗), the above equation can be deduced as follows: d𝜆₂(𝑇)d𝑇=𝑄₃𝜆₁+𝑄₅𝜆₂+𝑆^∗+𝑢₁−1𝑑₄+𝑏₃𝑑₄−𝑏₄𝑑₃𝑆^∗𝑄₁1+𝑏₃𝑆^∗𝑒^−𝛿𝑇, (67) where 𝑄₃=[(𝑏₄+𝑏₃+𝑏₃𝑐₁)𝑆^∗−𝑏₄+𝑐₁+1+𝑏₄𝑢₁]𝑆^∗/(1+𝑏₃𝑆^∗) and 𝑄₅=[𝑎₃+𝑒₁−(𝑐₁+𝑎₃)𝑆^∗].

According to (65), (67) can be rewritten as follows: d𝜆₂(𝑇)d𝑇=𝑄₅𝜆₂−𝑆^∗−1𝑄₁𝑄₄1+𝑏₃𝑆^∗𝑒^−𝛿𝑇, (68) where 𝑄₄=𝑑₃𝑄₃/(𝛿+𝑄₂)(1+𝑏₃𝑆^∗)−[𝑑₄+(𝑏₃𝑑₄−𝑏₄𝑑₃)𝑆^∗].

By solving (68), it can be obtained that 𝜆₂(𝑇)=𝑆^∗+𝑢₁−1𝑄₁𝑄₄1+𝑏₃𝑆^∗𝛿+𝑄₅𝑒^−𝛿𝑇. (69)

Based on (65) and (69), it can be obtained that d𝜆₃(𝑇)=ℎd𝑇¹𝛽₁𝑌𝑒^−𝛿𝑇𝛽₁+𝑃𝑌²+𝑆𝜆₁1+𝑏₃𝑆+𝑏₄𝐼+𝑎₃𝐼𝜆₂1+𝑏₃𝑆+𝑏₄𝐼+𝑓₁+𝑌𝜆₃. (70)

By associating with 𝑀^∗(𝑆^∗,0,𝑃^∗,𝑌^∗), (70) can be deduced as follows: d𝜆₃(𝑇)=ℎd𝑇¹𝛽₁𝑌^∗𝑒^−𝛿𝑇𝛽₁+𝑃^∗𝑌^∗2+𝑆^∗𝜆₁1+𝑏₃𝑆^∗+𝑑₃𝑆^∗1+𝑏₃𝑆^∗𝜆₃=𝑑₃𝑆^∗1+𝑏₃𝑆^∗𝜆₃+ℎ₁𝛽₁𝑌^∗𝛽₁+𝑃^∗𝑌^∗2+𝑑₃𝑆^∗1−𝑆^∗𝑄₁𝛿+𝑄₂1+𝑏₃𝑆^∗2𝑒^−𝛿𝑇. (71)

By solving the above equation, it can be obtained that 𝜆₃(𝑇)=ℎ₁𝛽₁𝑑₃𝑆^∗𝑆^∗−1−ℎ₁𝛽₁𝑌^∗𝛿+𝑄₂1+𝑏₃𝑆^∗2𝛿+𝑄₂1+𝑏₃𝑆^∗𝛽₁+𝑃^∗𝑌^∗2𝛿+𝑏₃𝛿+𝑑₃𝑆^∗×𝑒^−𝛿𝑇. (72)

In order to obtain an optimal equilibrium solution, it can be obtained by considering the disease-free equilibrium 𝑀^∗ and substituting (72) into (62), and it derives that 𝑔1−¹𝛽+𝑃^∗𝑌^∗2𝑃^∗ℎ₁𝛽₁=𝑑₃𝑆^∗𝑆^∗−1−𝑌^∗𝛿+𝑄₂1+𝑏₃𝑆^∗2𝛿+𝑄₂1+𝑏₃𝑆^∗𝛿+𝑏₃𝛿+𝑑₃𝑆^∗. (73)

Equation (73) provides an equation to the singular path and gives the optimal equilibrium levels of population 𝑆^∗=𝑆_𝛿,  𝑃^∗=(1−𝑢₁−𝑆_𝛿)(1+𝑏₃𝑆_𝛿), 𝑌^∗=((𝑑₃−𝑏₃𝑓₁)𝑆_𝛿−𝑓₁)/(1+𝑏₃𝑆_𝛿). Then, the optimal equilibrium levels of taxation and vaccination control can be obtained as follows: 𝑢_1𝛿=𝑏₃−1−𝑏₃𝛿+2𝑏₃𝑆_𝛿+𝑏₃−1−𝑏₃𝛿+2𝑏₃𝑆_𝛿²+4𝑏₃𝑑₃𝑆_𝛿𝑆_𝛿−1𝑒^−𝛿𝑇𝐵^1/2×2𝑏₃⁻¹,𝜎_𝛿=𝑔₁𝑑₃−𝑏₃𝑓₁−ℎ₁𝑏₃𝑆²_𝛿+ℎ₁1−𝑢₁+𝑔₁𝑓₁1−𝑢₁−𝛽₁×1+𝑏₃𝑆_𝛿1−𝑢₁−𝑆_𝛿⁻¹+𝑏₃−𝑏₃𝑢₁−1ℎ₁−𝑔₁𝑓₁+1−𝑢₁×𝑑₃−𝑏₃𝑓₁𝑆_𝛿×1+𝑏₃𝑆_𝛿1−𝑢₁−𝑆_𝛿⁻¹. (74)

Remark 12 (according to [13]). 𝜆_𝑖(𝑇)𝑒^𝛿𝑇(𝑖=1,2,3,4) represent unusual shadow prices along the singular path. From (65), (69), and (72), it may be concluded that these shadow prices remain constant over time interval in an optimum equilibrium when they strictly satisfy the transversality condition at ∞ [34]. Furthermore, they remain bounded as 𝑇→∞.

Generally, there is a strong ecological relationship between predator population and susceptible prey population in the harvested eco-epidemiological system, and these two species are usually influenced by each other in several aspects. The imposition of taxation and vaccination defined in (74) act as a deterrent to regulate harvest effort on predator population and control the spread of infectious disease within eco-epidemiological prey predator system, which can not only guarantee the sustainable development of predator population, but also benefit the ecological balance in the ecosystem and sustainable development of other relevant species in the harvested eco-epidemiological system.

5.2. Numerical Simulation

With the help of MATLAB, simulation work with hypothetical set of parameters is performed to understand the theoretical results which have been established in Section 5.1 of this paper. The values of parameters are cited from [9, 12], which are given as follows: 𝑟=1,  𝑘=1,  𝑎₁=1,  𝑏₁=2.34,  𝑏₂=3,  𝑐=0.01,  𝑎₂=0.344,  𝑒=0.051,  𝑑₁=0.21,  𝑑₂=0.6774,  𝑓=0.001,  𝛼=0.4,  ℎ=1.67,  𝛽=1.1,  𝑔=0.41,  𝛿=0.05,  and   𝐵=1.

According to the parameters given above, it follows from (22) that the range of the taxation 𝜎₁∈[0.4496,1.6336] and vaccination 𝑢₁∈[−1.3663,1), which guarantees that there exists a unique disease-free equilibrium. Based on (74), it can be computed that the optimal taxation and optimal vaccination control is 𝜎_𝛿=1.0653 and 𝑢_1𝛿=0.483, respectively. Since 𝜎_𝛿∈[0.4496,1.6336] and 𝑢_1𝛿∈[−1.3663,1), it can be obtained that there exists a unique disease-free equilibrium 𝑀^∗ corresponding to the optimal taxation and vaccination control, that is, 𝑀^∗(0.5648,0,0.9406,0.0501). According to Theorem 10, model system (20) is globally stable around 𝑀^∗, whose dynamical responses can be shown in Figure 1, and the corresponding phase portraits are shown in Figure 2 (𝑆-𝐼-𝑃 space), Figure 3 (𝑆-𝑌-𝑃 space), and Figure 4 (𝑌-𝑃-𝐼 space) with different initial values.

Figure 1 

Dynamical responses of model system (20) with the optimal taxation and vaccination 𝜎𝛿=1.0653 and 𝑢1𝛿=0.483.

Figure 2 

Phase portrait (𝑆-𝐼-𝑃 space) of model system (20) with 𝜎𝛿=1.0653, and 𝑢1𝛿=0.483 and differential initial values, which indicates model system (20) is globally stable around 𝑀∗(0.5648,0,0.9406,0.0501).

Figure 3 

Phase portrait (𝑆-𝑌-𝑃 space) of model system (20) with 𝜎𝛿=1.0653 and 𝑢1𝛿=0.483 and differential initial values, which indicates model system (20) is globally stable around 𝑀∗(0.5648,0,0.9406,0.0501).

Figure 4 

Phase portrait (𝑌-𝑃-𝐼 space) of model system (20) with 𝜎𝛿=1.0653 and 𝑢1𝛿=0.483 and differential initial values, which indicates model system (20) is globally stable around 𝑀∗(0.5648,0,0.9406,0.0501).

6. Conclusion

It is well known that studies on ecology and epidemiology share some common features. Nowadays, taxation is extensively used as economic control instrument to protect the population of the prey predator system from overexploitation, and vaccination is widely adopted to control spread of infectious disease within harvested eco-epidemiological system. Based on superiority of taxation control in regulating harvesting and effectiveness of vaccination in protecting population from infectious disease, it motives the introduction of taxation and vaccination control into model proposed in this paper. In this paper, a harvested eco-epidemiological model is established to investigate the dynamics effects of vaccination and taxation control on the population dynamics. We extend the work in [9] by incorporating the harvest effort into eco-epidemiological prey predator system. Vaccination and taxation are chosen as controlling instruments to protect population from infectious disease and regulate harvesting on predator, respectively.

Since it is of more biological significance to consider the effect of interacting species when it comes to discussing dynamics of epidemiological model, conditions that influence nonnegativeness and boundedness of solutions are discussed, which can be found in Theorems 3 and 4. Local stability analysis around disease-free equilibrium is studied in Theorem 7. Based on robust Bendixson criterion, global stability of model system around the unique disease-free equilibrium is discussed in Theorem 10, which is theoretically beneficial to studying the coexistence and interaction mechanism of harvested eco-epidemiological prey predator system. By using Pontryagin’s maximum principle, optimal vaccination and taxation control are derived to ensure the sustainable development of biological resource and protect population from infectious disease in the fifth section of this paper. It reveals that the total user’s cost of harvest per unit effort must be equal to the discounted value of the future price at the steady-state level, and infectious disease among prey population can be removed with optimal vaccination strategy on susceptible prey population.

It should be noted that the optimal control is obtained based on an interesting reflection, which discusses the impact and efficiency of optimal control and deducts a tendency for the management of biological resource. Since infectious disease spreads within the eco-epidemiological prey predator system and the commercial harvesting exists extensively in the natural world, it is necessary to discuss optimal strategy associated with net harvesting revenue and vaccination cost. Hence, theoretical analysis obtained in the fifth section is constructive for regulatory agency to formulate constructive policies and maintain the sustainable harvesting and effective vaccination by means of optimal taxation and vaccination control; it makes the work studied in this paper have some new and positive features.

Acknowledgments

The authors are thankful to the learned reviewer and editors for his/her useful suggestions and comments. 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, supported by Hong Kong Admission Scheme for Mainland Talents and Professionals, Hong Kong Special Administrative Region.