Abstract and Applied Analysis

Volume 2013 (2013), Article ID 950396, 16 pages

http://dx.doi.org/10.1155/2013/950396

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

^{1}Institute of Systems Science, Northeastern University, Shenyang 110004, China^{2}State Key Laboratory of Integrated Automation of Process Industry, Northeastern University, Shenyang 110004, China^{3}Shenyang Institute of Automation, Chinese Academy of Sciences, Shenyang 110004, China

Received 28 April 2013; Accepted 2 July 2013

Academic Editor: Luca Guerrini

Copyright © 2013 Chao Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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: 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 and intrinsic growth rate constant 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 , where 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 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 respectively, where represents the harvest effort on predator population at any time , , , and are all positive constants. To conserve the population in the harvested eco-epidemiological prey predator system, the regulatory agency imposes a tax per unit biomass of the predator, and 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 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: The initial conditions of model system (4) are given as follows: 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 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:
where and is given as follows:Let be the nonnegative octant in , then is locally Lipschitz and satisfy the condition .

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 for some . 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 and , then all solutions of model system (4) are bounded within a region for all ,
**On the other hand, if and , then all solutions of model system (4) are bounded within a region for all ,
**
where , , and .*

*Proof. *Firstly, we consider the first equation of model system (4). We have
which implies that

Now, defining the function by calculating the time derivative of along the solutions of model system (4), we get
By utilizing the standard comparison principle, it follows that

Furthermore, by using nonnegativity of solution of model system (4), if , then it derives that
which implies that

Let , and by calculating the time derivative of along the solutions of model system (4), it can be obtained that

Hence, it derives that
where .

Consequently, If and , then all solutions of model system (4) are bounded within a region for all ,

On the other hand, if and , then all solutions of model system (4) are bounded within a region for all ,

*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: where , , , , , , , , , , , , and .

By simple computation, disease-free equilibrium can be obtained as follows.

, , and satisfies the following equation:

Conditions for existence of unique positive real root of (21) are as follows:

*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
**
then model system (20) is locally stable around disease-free equilibrium .*

*Proof. *For disease-free equilibrium , the variational matrix of model system (20) around takes the following form:
where

It follows from that the characteristic equation around is as follows:
where

It follows from Routh-Hurwitz criteria [13] that there are four eigenvalues with negative real part provided that , and , which derives that

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 be the second additive compound matrix [31, 32]. In the case of , can be written based on variational matrix of model system (20) around , which is as follows:where , , can be found in Theorem 7. Since , , , , , and , can be deduced as follows:

Furthermore, we have where and are defined based on model system (20), and the model system (20) can be interpreted as the matrix form where and is given as follows:

By simple computation, it can be obtained that

Let where

Consider the Lozinskii measure of with respect to a vector norm in (see, [32]) and (identity matrix), then it can be obtained that

Hence, it follows from simple computation that where for , and . , , , , , and are matrix norms with respect to the vector norm, and is the Lozinskii measure with respect to this vector norm.

Furthermore, let , then it can be computed that

According to Theorem 3 and definition of , it follows that

Lemma 9. *For element of the disease-free equilibrium , some inequalities about hold as follows: *(a)* if , then
*(b)* if , then(i) − , provided that ,(ii), provided that ,*

*where and .*

*Proof. *Based on Theorem 4 and nondimensionalized scaling of model system (20) mentioned in Section 4, it derives that

According to model system (20), let . By calculating the time derivative of along the solutions of model system (20), it can be obtained that
where .

If , then

According to , it follows from the standard comparison principle that

On the other hand, if , then
where .

According to , it follows from the standard comparison principle that
provided that .

Furthermore,
provided that .

Let

Theorem 10. *If the model system is locally stable around and inequality (22) holds, then*(a)*if and , then model system (20) is globally stable around ; *(b)*if , then(i) model system (20) is globally stable around provided that , + ,(ii)model system (20) is globally stable around provided that , + ,*

*where , and , have been defined in Lemma 9. , 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 and , then
which derives that

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 , , and , then
which derives that

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 , , and , then
It derives that

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

By using the Pontryagin’s maximum principle [12], the associated Hamiltonian function is constructed by where , are adjoint variables. and are the control variables. According to [33], the conditions for a singular control to be optimal can be obtained by and , from which we get

According to Theorem 3, it follows that

For adjoint variables , we have

By substituting into the Equation , it follows that

By virtue of (62), it gives that

By associating with , the above equation can be deduced as follows: where and .

By solving (64), it can be obtained that

By virtue of (62), it gives that

By associating with , the above equation can be deduced as follows: where and .

According to (65), (67) can be rewritten as follows: where .

By solving (68), it can be obtained that

Based on (65) and (69), it can be obtained that

By associating with , (70) can be deduced as follows: