Dynamics analysis of a predator–prey system with harvesting prey and disease in prey species

ABSTRACT In this paper, a predator–prey system with harvesting prey and disease in prey species is given. In the absence of time delay, the existence and stability of all equilibria are investigated. In the presence of time delay, some sufficient conditions of the local stability of the positive equilibrium and the existence of Hopf bifurcation are obtained by analysing the corresponding characteristic equation, and the properties of Hopf bifurcation are given by using the normal form theory and centre manifold theorem. Furthermore, an optimal harvesting policy is investigated by applying the Pontryagin's Maximum Principle. Numerical simulations are performed to support our analytic results.


Introduction
In the twentieth century, the pioneer work of Lotka [27] and Volterra [36] opened new door to biological species. Many researchers have made a lot of achievements in this field (see [13,29,31,32] and references therein). In 1927, Kermack and MacKendrick [22] proposed the classical Susceptible, Infectious, Recovered model which has attracted more scholars attention (see [3,37,45] and references therein). Ecology and epidemiology are two different and independent fields for years. However, as time goes by, these two fields are more and more closer, and a new cross field which is called eco-epidemiology emerges. In 1986, Anderson and May [20] firstly investigated the dynamics of eco-epidemiology model where the predator interacts with infected prey population.
In recent years, the eco-epidemiology problem has been pained attention to many scholars and experts. Generally speaking, prey-predator models could be divided into three cases by the infective disease in population. The first one is that there is only the infected prey in models (see [3][4][5][6][7]14,19,23,34,45]). Johri et al. [19] studied a Lotka-Volterra type predator-prey model with disease in prey and without harvesting, and assumed that the conversion rate of the susceptible prey is the same to that of the infected prey. Sharma and Samanta [34] studied an eco-epidemiology model with two prey population where one prey species is infected by an infectious disease. The second one is that there is only the infected predator in models [18,40,42,44]. Zhang and Sun [42] proposed a predator-prey model with disease in the predator and general functional response. The last one is disease in both populations for prey-predator models [9,21]. Kant and Kumar [21] considered a predator-prey system with migrating prey and disease infection in both populations.
It is well known that time delays of one type or another have been incorporated into mathematical models of population dynamics by many scholars. In general, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause unstable of the equilibrium. Time delay due to incubation period of infectious disease is a common example. Hu and Li [14] investigated a delayed predator-prey model. Xu and Zhang [40] considered a delayed predator-prey model with disease in the predator. Zhou et al. [44] proposed a three-dimensional eco-epidemiological model with delay.
In addition, many of the species in the natural ecosystem are being mined and harvested. For instance, in an Irish sea, at least four populations of skate and shark have been exploited to extinction [2]. Naturally, optimal harvesting policy is very much needed for such system (see [1,[4][5][6]8,9,18,23,25,28] and references therein). Biswas et al. [4] studied an eco-epidemiological system with weak Allee effect and harvesting in prey population, and found that optimal harvesting policy is dependent on both gestation delay and Allee effect. Chattopadhy et al. [6] investigated a harvesting predator-prey model with infection in the prey population and concluded that harvesting on infected prey prevents the limit cycle oscillations. Khan et al. [23] investigated a predator-prey model with density constraints for susceptible prey population, and considered harvesting to both susceptible and infected prey species. Huang et al. [15] presented a predator-prey system of Holling and Leslie type with constant-yield prey harvesting, and obtained various bifurcations, such as saddle-node bifurcation, Hopf-bifurcation, repelling and attracting Bogdanov-Takens bifurcation. For the other models with optimal harvesting, we refer to [15][16][17]38,39] and the references cited therein.
In our model, the predator not only captures the infective prey but also catches the susceptible prey based on the work of Biswas et al. [4]. Due to incubation period of infectious disease [14,44], we will consider the affection of time delay in our model. In addition, the harvesting of ecological resources and population is a common behaviour, such as harvest of wild animals and fish. From the perspective of reality, capture cannot be ignored in the ecosystem. Based on the above analysis and work of Johri et al. [19], we will propose and analyse an eco-epidemiological model with time delay, harvesting prey population and constant input of susceptible prey population. More detailed assumptions and descriptions can be seen in the next section. This paper is organized as follows. In Section 2, we will formulate an ecoepidemiological model with harvesting prey, time delay and disease in prey species, and discuss the basic properties of system (1). In Section 3, by analysing the corresponding characteristic equations, we will investigate the local stability of the non-negative equilibria of system (1), and the global stability of system (1) at the disease-free equilibrium and the positive equilibrium, respectively. In Section 4, we will show that Hopf bifurcation occurs in delayed model when the delay passes though its critical value and show the properties of the Hopf bifurcation. We will constitute the optimal control problem by choosing effort as the control variable, and give the optimal harvesting policy of system by using Pontryagin's Maximum Principle in Section 5. To support our theoretical analysis, some numerical simulations are given in Section 6. We end with a brief conclusion in Section 7.

Model description
Let N(t) be the population density of the total prey species and P(t) be the population density of the predator species. Now, in order to formulate our eco-epidemiological model, some following assumptions are given: (A 1 ) The total prey is composed of two classes: susceptible prey (S) and infective prey (I).
Therefore, the total biomass of the prey is S(t) + I(t). (A 2 ) The disease spreads among the prey population only by contact, and can not be transmitted vertically. The infection prey species do not recover or become immune. The incidence of the disease in the prey is bilinear incidence form βSI. (A 3 ) The susceptible prey and the infective prey are both captured. Generally speaking, susceptible prey are more stronger than infected prey. Therefore, the probability of predation of the infected prey is more bigger than that of susceptible prey. That is, The death rate of the population is considered, where the death rate of the infective prey contains natural mortality and mortality due to infection of infected prey. (A 5 ) The prey species is harvested by a linear rate.
The model structure is shown in Figure 1. The transfer diagram leads to the following system of ordinary differential equations and the initial conditions are given as We assume all parameters of system (1) are positive. The detailed biological meanings of parameters are given in Table 1. In addition, due to incubation period of infectious disease, we consider system (1) with feedback delay, which makes the model more appropriate to biological environment. Therefore, system (1) can be rewritten as the following:  The constant recruitment rate of susceptible prey β The disease transmission coefficient α 1 Capture rate of the susceptible prey α 2 Capture rate of the infected prey d 1 The death rate of the susceptible prey d 2 The death rate of the infected prey d 3 The death rate of the predator E The harvesting effort e Conversion coefficient from the prey to the predator q 1 The catchability coefficient of susceptible prey q 2 The catchability coefficient of infected prey where τ ≥ 0 is time interval of infected prey species. The rest of the parameters are the same as system (1). The initial conditions of system (2) are given by

Basic properties
In this subsection, we will show the basic properties which are necessary for the understanding of subsequent results. Therefore, we have the following lemmas. Proof: Let S(t), I(t) and P(t) be any solution of system (1). Assuming that one solution of the system (1) is at least not positive, then we have the following three cases: (i) there exists a t 1 , such that (ii) there exists a t 2 , such that (iii) there exists a t 3 , such that If case (i) holds on, then we obtain This contradicts with S (t 1 ) < 0. If case (ii) holds on, then we obtain This contradicts with I (t 2 ) < 0. If case (iii) holds on, then we obtain This contradicts with P (t 3 ) < 0. From the arbitrariness of S(t), I(t) and P(t), all solutions of system (1) remain positive for all t > 0.

Proof: Defining a function
we have that, Further, we can get If we suppose that μ = min{d 1 + q 1 E, d 2 + q 2 E, d 3 } and 0 < e < 1, then we get The solution of the ordinary differential equation dZ/dt + μZ = is where C is an arbitrary constant. Therefore, we have It implies that the region = (S(t), I(t), P(t)) ∈ R 3 + : is a positively invariant set of system (1). Therefore, we will consider dynamics of system (1) on the set in this paper. The proof is completed.

Boundary equilibria and their stability
In this subsection, we will study the existence of the boundary equilibria and their stability of the system (1) by analysing the corresponding characteristic equations, respectively. In order to make convenience of computation, some notations are given as follows , Firstly, we will analyse the existence of the boundary equilibria of the system (1). System (1) always has a trivial equilibrium E 0 (0, 0, 0) and an axial equilibrium E 1 ( d 1 +q 1 E , 0, 0).
If R 00 > 1, then system (1) has a predator-extinction equilibrium E 2 (S,Ī, 0), wherē If R 01 > 1, then system (1) has a disease-free equilibrium E 3 (S, 0,P), wherẽ Nextly, we will investigate the stability of the boundary equilibria of the system (1). The Jacobian matrix of system (1) is given by With help of this matrix we will analyse the local stability of equilibria in our model. It is easy to prove that the trivial equilibrium E 0 (0, 0, 0) is always stable. The eigenvalues of the Jacobian matrix at E 1 are It is clear that the axial equilibrium E 1 is locally asymptotically stable when R 00 < 1 and R 01 < 1.
One of the eigenvalues of the Jacobian matrix at E 2 is λ 1 = eα 1S + eα 2Ī − d 3 , and the others satisfy the following equation where Thus the two roots of the Equation (3) are negative. Therefore, the predator-extinction equilibrium E 2 is locally asymptotically stable when R 02 < 1. From the above discussion, we have the following results.
Theorem 3.1: (i) The trivial equilibrium E 0 is locally asymptotically stable for all parameters and the axial equilibrium E 1 is a locally asymptotically stable if R 00 < 1 and R 01 < 1; (ii) The predator-extinction equilibrium E 2 is locally asymptotically stable if R 02 < 1.

Remark 3.1:
The number of susceptible prey tends to be constant, and the predator tends to extinction, when R 01 < 1; the susceptible prey and the predator populations are coexistence, when R 01 > 1. Consequently, the parameter R 01 controls the coexistence of the prey and the predator. Proof: Firstly, we prove that the disease-free equilibrium E 3 is locally asymptotically stable. One of the eigenvalues of the Jacobian matrix at E 3 is λ 1 = βS − α 2P − d 2 − q 2 E, and the others satisfy the following equation: Nextly, we prove that the disease-free equilibrium E 3 is globally asymptotically stable. By the second equation of system (1), we obtain Structural comparison equation By the comparison theorem, then we can get Hence, we consider the limit system of the system (1) as follows It is easy to prove that is locally asymptotically stable.
If we want to prove that E 3 is globally asymptotically stable in 1 = {(S, P) ∈ R 2 + : S + P ≤ M}, we only need to prove that limit system has no periodic orbit in 1 . Now, we chose a Dulac function B(S, P) = P −1 , and denote two functions By calculation, we can obtain The Bendixon-Dulac criterion [24] shows that the system (4) has no periodic orbit in 1 , and the equilibrium E 3 is globally asymptotically stable in 1 . Therefore, the disease-free equilibrium E 3 is globally asymptotically stable.

Positive equilibrium and its stability
In this subsection, we will discuss the stability of the positive equilibrium. It is easy to show that if the following assumption holds, then system (1) has a positive equilibrium E 4 (S * , I * , P * ), where .

Remark 3.2:
If R 03 > 1, then the disease-free equilibrium point E 3 of system (1) is unstable and there exists the positive equilibrium E 4 , that is, disease will be endemic. R 03 is equal to the basic reproductive number in epidemic model.
From the above analysis, we can get the following conclusions.
Proof: Firstly, we show that the positive equilibrium E 4 is locally asymptotically stable. The characteristic equation of the Jacobian matrix at E 4 is given by It is clear that C 1 > 0, C 2 > 0 and C 3 > 0. Consequently, according to the Routh-Hurwitz criterion, E 4 will be locally asymptotically stable if the assumption (H2) is true. Nextly, we will prove that the positive equilibrium E 4 is globally asymptotically stable. We define the following Lyapunov function Differentiating of V along the solutions of system (1) with respect to t, we obtaiṅ Substituting Equation (1) into Equation (5), we can geṫ The maximum invariant set in {(S, I, P) :V = 0} is the singleton. By LaSalle's invariance principle in [26], E 4 is globally asymptotically stable.

Local stability and existence of Hopf bifurcation
From the previous analysis, we know the system (2) has a positive equilibrium E 4 (S * , I * , P * ).
Linearizing system (2) at the positive equilibrium E 4 , we have the following system where Then, we can get that the characteristic equation of system (6) at E 4 (S * , I * , P * ) is given by where Firstly, if τ = 0, then Equation (7) can be written as By the Routh-Hurwitz criteria, all roots of Equation (8) must have negative real parts if the assumption (H3) m 0 + n 0 > 0, m 1 + n 1 > 0, m 2 + n 2 > 0 and (m 2 + n 2 )(m 1 + n 1 ) > (m 0 + n 0 ) is satisfied. Then the positive equilibrium E 4 of system (2) is locally asymptotically stable when τ = 0 if the condition (H3) holds on.
Secondly, if τ = 0, let λ = iω(ω > 0) be the solution of Equation (7). Therefore, we can obtain the following form Separating the real and imaginary parts, we have It follows that ω satisfies Defining u = ω 2 , Equation (11) can be written as where Discussion about the roots of Equation (12) is similar to that in [35], then we have the following lemma. Furthermore, we assume that the coefficients in f (u) satisfy the condition (H4) p 0 < 0 or p 0 ≥ 0, > 0 and f (u * ) ≤ 0. If condition (H 4 ) holds on, then Equation (12) has at least one positive root. Without loss of generality, we may assume that it has three positive roots, which are denoted as u 1 , u 2 and u 3 , respectively. Then Equation (11) has three positive roots ω k = √ u k , k = 1,2,3.
Through the above analysis and the Hopf bifurcation theorem in [12], we have the following theorem. (H3) and (H4) hold on, then system (2) at the positive equilibrium E 4 is locally asymptotically stable when τ ∈ [0, τ 0 ), and unstable when τ > τ 0 . That is, system (2) undergoes a Hopf bifurcation at the positive equilibrium E 4 when τ = τ 0 .

Properties of the Hopf bifurcation
According to the above analysis, we know that system (2) undergoes a Hopf bifurcation at the positive equilibrium E 4 when τ = τ 0 . In this subsection, we will investigate the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions of system (2) by using the normal form approach theory and centre manifold theorem due to Hassard et al. [12].
We assume that μ = τ − τ 0 is a new bifurcation parameter. Then μ = 0 is the Hopf bifurcation value. Let Then system (2) is transformed into a functional differential equation where L μ : C → R 3 and F : R × C → R 3 are defined, respectively, by By the Reisz representation theorem, there exists a 3 × 3 matrix function In fact, we choose where δ is the Dirac delta function. and Then system (14) can be rewritten aṡ where η(θ) = η(θ, 0). From the above analysis, we have the following lemma.  (16), (19) and (20), we can get Namely, A = A(0) and A * are adjoint operator.
From the above discussion, we know that ±iω 0 are eigenvalues of A(0) and ∓iω 0 are eigenvalues of A * (0). Thus, we have the following lemma.
Proof: Let q(θ ) be the eigenvector of A(0) corresponding to iω 0 and q * (s) be the eigenvector of A * (0) corresponding to −iω 0 . It means A(0)q(θ ) = iω 0 q(θ ) and A * (0)q * T (θ ) = −iω 0 q * T (θ ). From Equations (16) and (19) we can get Then By a simple calculation, we can get v 1 = (a 12 a 23 − a 13 a 22 ) + a 13 Since By a simple calculation, we can get Since From Equation (20), we can obtain Next, by using the same notation as in Hassard et al. [12] and using a computation process similar to that in [30], we calculate the coordinates of the centre manifold C 0 at μ = 0. Let u t be the solution of Equation (14) when μ = 0. We define On the centre manifold, we have Since z andz are local coordinates for the centre manifold in the direction of q * andq * . Note that W is real if u t is real. So, we consider only real solutions. Thuṡ where g(z,z) = g 20 z 2 2 + g 11 zz + g 02z According to Equations (18) and (21), we havė On the other hand, on C 0Ẇ Substituting Equations (23) and (24) Since u t = x(t + θ) = W(z,z, θ) + zq +zq, then ⎛ ⎜ ⎝ that is, where Thus, we have that Sinceq * T =D(v * 1 , 1,v * 2 ), we can obtain Therefore, we have In order to determine g 21 , we need to calculate W 20 (θ ) and W 11 (θ ). For θ ∈ [−τ , 0), we can get Comparing the coefficients of the above equation with (25), we have From Equations (18), (27) and (29), we havė Therefore, we can get Similarly, since Equations (18), (27) and (29), we obtain Where 1 , E Now, we will seek appropriate E 1 and E 2 in Equations (30) and (31). In H(z,z, θ), if we assume that θ = 0, we have Comparing the coefficients of the above equation with (26), we can get From Equations (15), (16) and (27), we can obtain Notice that Substituting Equations (30) and (33) into the first equation of (32), we have It follows that Similarly, substituting Equations (31) and (34) into the second equation of (32), we can get Based on the analysis above, we can find that g 20 , g 11 , g 02 , g 21 are determined by the parameters and delay in system (2). Therefore, we can compute the following quantities: Here, the sign of μ 2 determines the stability of the bifurcating periodic solutions. The sign of β 2 determines the direction of the Hopf bifurcation. The sign of T 2 determines the period of the bifurcating periodic solutions. Furthermore, according to the above analysis and using the method of Hassard et al. [12], we can get the following theorem. (35),

The existence of optimal control
In this subsection, we will study the optimal harvesting policy for the delayed model. In business development of renewable resources, the fundamental problem from the economic point of view is to determine the optimal trade-off between the present and future harvests. Therefore, we will investigate the optimal harvesting policy by considering the profit earned by harvesting, and focusing on quadratic cost [33] and conservation of prey species. The main reason for using quadratic costs is that it allows us to derive an analytical expression for the optimal harvesting, whose result is different from the Bang-Bang control obtained in the case of a linear cost function. It is assumed that the price is a function which is inversely proportional to the available prey population, that is, the price function decreases when prey population increases. Therefore, we give the optimal harvesting as follows where c, p 1 and p 2 are the constant harvesting cost per unit effort, the constant price per unit biomass of the susceptible and infected prey species, respectively. v 1 and v 2 are economic constants and δ is the instantaneous annual discount rate.
The problem (36), subject to system (2) and control constraints 0 ≤ E ≤ E max , where E max is the maximum effort capacity of the harvesting industry. Thus, our object is to find an optimal control E * such that where U is the control set defined by In order to verify these conditions, we use the result in Fleming and Rishel [10]. We note that the solutions are bounded. The set of all the control variables E(t) ∈ U is also convex and closed by definition. The optimal system is bounded which determines the compactness needed for the existence of the optimal control. In addition, the integrand in Equation (36) is given by e −δt [(p 1 − v 1 q 1 ES)q 1 ES + (p 2 − v 2 q 2 EI)q 2 EI − cE] and is convex on the control set U. Moreover, the integrand in Equation (36) can be written as we can easily see that there exist a constant α > 0 and positive numbers d 1 and d 2 such that which completes the existence of an optimal control.

The characterization of optimal control
To find the maximum value of J(E), we assume that x(t) = (S(t), I(t), P(t)) and x τ (t) = (S τ (t), I τ (t), P τ (t)), where S τ (t) = S(t − τ ), I τ (t) = I(t − τ ), P τ (t) = P(t − τ ). By using the Pontryagin's Maximum Principle to control the optimal harvesting [1], we formulate the Hamiltonian H(t) = H(x, x τ , E, λ) given by where λ 1 , λ 2 and λ 3 are the adjoint variables. λ 1 (t f ) = 0, λ 2 (t f ) = 0, and λ 3 (t f ) = 0 are transversality conditions. In order to obtain the necessary condition for the optimal control, Pontryagin's Maximum Principle with delay is given in [41]. The state equation are given by the optimality condition is and the adjoint equation is where H λ , H E , H x and H x τ denote the derivative with respect to λ, E, x and x τ , respectively. Now, we apply the necessary conditions to the Hamiltonian H(t) in Equation (37).

Theorem 5.2:
Let S * (t), I * (t), P * (t) be optimal state solutions associated with the optimal control variable E * (t) for the corresponding state system (2), Then, there exist adjoint variables λ 1 , λ 2 and λ 3 satisfying with transversality condition (or boundary conditions) Furthermore, optimal control pairs are given as follows Proof: In order to determine the adjoint equations and the transversality conditions, we use the Hamiltonian (37). By using the adjoint Equation (39) and differentiating the Hamiltonian (37) with respect to x(t) and x τ (t), letting x(t) = x * (t) and x τ (t) = x * τ (t). Then, the adjoint system can be written as Through the calculation, it is easy to verify the adjoint Equation (40). Using the optimality conditions (38), we obtain This implies that .
Using the property of control set in U, we have Therefore, we have the the optimal control E * (t) is given by , E max .

Numerical simulations
To observe the rich dynamical behaviour, we will take different parameter values in Table 2.
Finally, we will investigate numerically illustrate the optimal solutions to system (1) by numerical method from [11]. We use fourth-order Runge-Kutta forward iterative method to solve the state variables of the system (1) and then we solve the system (37) for the adjoint variables by backward fourth-order Runge-Kutta iterative method [11]. Now, we choose a set of the following parameter values: = 0.49; β = 0.75; α 1 =  Figure 5 shows the solution curves of the state variables in the presence of the optimal harvesting. Figure 6 represents the graphs of adjoint variables λ 1 , λ 2 and λ 3 , respectively. It is interesting to observe that all the adjoint variables are monotone with time increasing(λ 1 is monotonically decreasing, on the contrary, λ 2 and λ 3 are monotonically increasing) and their final values are zero. From Figure 7(a), we observe that the optimal harvesting effort is minimum and tends to be stable when α 1 = 0.8. From the viewpoint of biology, if α 1 increases, then the value of α 1 SP will increase, which leads to cut down the maximum optimal harvesting effort. To see it more clearly, we cut out a paragraph in the time axis arbitrarily as Figure 7(b) shown. In addition, we choose β = 0.35, e = 0.5, τ = 1 and assume that 0 < E < 1. We investigate the optimal harvesting policy for system (2). Figure 8 shows the optimal trajectories for system (2) with population densities (S(0), I(0), P(0)). Figure 9 describes the variations of adjoint variables. From Figure 10(a), we can find that the optimal harvesting effort is strongly dependent on incubation period of infectious disease. Maximum value of optimal  harvesting decreases as time delay τ increases. The optimal harvesting effort decreases monotonically. To see it more clearly, we cut out a paragraph in the time axis arbitrarily as Figure 10(b) depicted. From Figure 11(a) and 11(b), it is clear that the effect of α 2 on the predator is greater than that of α 1 on predator. It means that the infected prey are more to be caught.

Conclusion
In this paper, under the assumption that the predator populations can still survive in the disease spread in the prey populations, a delayed eco-epidemiology model with harvesting prey is proposed. For the model without time delay, by the analysing the corresponding characteristic equations, the local stability of each of feasible equilibria of system (1) has been established, respectively. Theorem 3.2 showed that the disease would be eliminated, when R 03 < 1. Theorem 3.3 showed that the disease would epidemic when R 03 > 1. For the model with time delay, the time delay plays an important role. Time delay could switch the stability of the equilibrium form stable to unstable. That is, the positive equilibrium E 4 is stable when τ < τ 0 (Theorem 4.1), but it can lose its stability and a Hopf bifurcation occurs at the positive equilibrium when the delay crosses the critical value τ 0 . Furthermore, we have derived an explicit algorithm and sufficient conditions for the stability of the bifurcating periodic solutions by using the normal form method and centre manifold theorem. From the biological point of view, the time delay has influence on the transmission of the infectious diseases, that is, the time delay is harmful. In addition, we have investigated an optimal harvesting policy by using the Pontryagin's Maximum Principle. We could find that the optimal harvesting policy is dependent on both the incubation period of infectious disease and the capture rate of the susceptible prey. To verify the analytical results obtained in this paper, we also have given suitable numerical simulations and numerical solutions for a special case of system under consideration.
In our model, we have assumed that the incidence of the disease in the prey is bilinear incidence form βSI. However, as the number of susceptible prey is large, owing to the number of susceptible prey with which every infective contacts within a certain time being limited, it is unreasonable to consider the bilinear incidence rate [43]. It is an interesting thing to consider the following model with saturation incidence and double delay.
where α, τ 1 and τ 2 are the semi-saturated coefficient, the incubation period of infectious disease and the gestation period of the predator, respectively. We leave these work in the future.