Hopf Bifurcation and Optimal Control in a Diffusive Predator-prey System with Time Delay and Prey Harvesting *

In this paper, we investigated the dynamics of a diffusive delayed predator-prey system with Holling type II functional response and nozero constant prey harvesting on no-flux boundary condition. At first, we obtain the existence and the stability of the equilibria by analyzing the distribution of the roots of associated characteristic equation. Using the time delay as the bifurcation parameter and the harvesting term as the control parameter, we get the existence and the stability of Hopf bifurcation at the positive constant steady state. Applying the normal form theory and the center manifold argument for partial functional differential equations, we derive an explicit formula for determining the direction and the stability of Hopf bifurcation. Finally, an optimal control problem has been considered.

Recently, as the continuous exploitation of biological resources, the harvest of population in fishery, forestry, and wildlife management become more crucial.Concerning the conservation for the long-term benefits of humanity, there is a wide-range of interest in the use of bioeconomic modeling to gain insight in the scientific management of renewable resources like fisheries and forestry [13].In biological, both harvesting in prey and in predator could lead some dangers in real-life harvesting such as no equilibrium exists and either prey or predator goes to extinction for some values of harvesting rate.It would be interesting that the nonzero constant harvesting in both models can prevent mutual extinction of prey and predator, and remove the singularity of the origin, which was regarded as "pathological behavior" by some researchers (see [14,15]).And in mathematics, the study of predator-prey models with harvesting have attracted the attention of many researchers which have rich bifurcating phenomenons (see e.g.[13,14,[16][17][18] and references therein).
Xiao and Jennings [13] considered the Michaelis-Menten type (ratio-dependent) predator-prey model with a constant rate prey harvesting.It has been showed that system has four equilibria and undergoes two saddle-node bifurcations, the separatrix connecting a saddle and a saddle-node bifurcation and heteroclinic bifurcation.When the spatial dispersal is considered in ratio-dependent predator-prey model, system has been investigated by Zhang [19].They derived the conditions for Hopf and Turing bifurcation on the spatial domain.Their results have shown that modeling by reaction-diffusion equations is an appropriate tool for investigating fundamental mechanisms of complex spatiotemporal dynamics.
Many researchers focus on constant-yield harvesting of a predator-prey system (1) with Holling type II functional response and its variants (see [20][21][22][23][24]), where h 1 > 0 and h 2 > 0 are harvesting or removal rate for the prey and the predator, respectively.They observed very rich and interesting dynamical behaviors such as the existence of multiple equilibria, homoclinic loop, Hopf bifurcation etc. Brauer and Soudack [25] developed the theory of global behavior of system (1) and shown by examples which of the theoretically possible transitions can actually occur for a class of biologically motivated models.Brauer and Soudack [26] carried out the phase portrait analysis and considered the general form of system (1).They have shown how to approximate the region of asymptotic stability in biological terms the initial states which lead to coexistence of the two species by efficient computer simulations and how to identify values of the harvest rates for which the region of asymptotic stability disappears corresponding to collapse of the biological system for every initial state.If a time delay in the predator response term is considered in system (1) with h 1 = H and h 2 = 0, then it becomes ( The time delay represents a gestation time of the predators.The reproduction of predating the prey is not instantaneous, but will be mediated by some discrete time lag required for www.mii.lt/NAgestation of the predators (see also [27]).Xia et al. in [27] have investigated system (2), applying the normal form theory of retarded functional differential equations developed by Faria and Magalhães [28].Also, Martin and Ruan in [16] considered dynamics of the general model of system above.
If system (2) is further considered a spatial diffusion with no-flux boundary conditions, then it can be described by the following system: Here the habitat Ω is a bounded domain in R n (n ≥ 1), ν is the outer normal direction for x ∈ ∂Ω and the imposed no-flux boundary condition means the system is a closed one.Also, r is the intrinsic growth rate of the prey; K is the carry capacity of the prey; A is the half saturation constant for the predators which is the prey density at which the functional response is half maximal; B is the death rate of predators.System (3) with H = 0 has been investigated by many researchers (see e.g.[1,[29][30][31][32][33] and references therein).System (3) with H = 0 have no results as far as we know.We assume that the predator in system (3) is not of commercial importance.The prey is continuously being harvested at a constant rate by a harvesting agency.The harvesting activity does not affect the predator population directly.It is obvious that the harvesting activity does reduce the predator population indirectly by reducing the availability of the prey to the predator.The rest of the paper is organized as follows.In Section 2, we study the distribution of the roots of the characteristic equation and the stability of a positive constant steady state and show the conditions under which Hopf bifurcation occurs.In Section 3, we derived an algorithm for determining the direction and stability of the Hopf bifurcation by applying the normal form theory and the center manifold theory of partial functional differential equations ( [34] and [35]).Some numerical analysis is given which support the analysis of Section 2 and Section 3 in Section 4. In Section 5, an optimal harvesting policy is considered by using Pontryagin's maximal principle.

The boundary equilibria
Now, we consider the stability of two bounded equilibria in 0 < H < H 1 .

Consider
Linearizating system (3) around E 1 , we get the eigenvalues are Obviously, the eigenvalues are independent of the time delay, and the characteristic equation is coincide with system (3) with τ = 0. u + > K/2 implies λ 1n < 0. From and H = H 0 , then u + = u * , λ 1n < 0, λ 20 = 0, and λ 2n < 0, whose stability is determined by center manifold theory and normal form theory, and which is postponed to Appendix.Hence, we get the following lemma.
Theorem 2. For the boundary equilibrium E 1 of system (3), we have:

Consider
Using the same method above, we linearize the system around E 2 , and the eigenvalues are This implies that E 2 is unstable.

The interior equilibrium
Now, we focus on considering the dynamics behaviors of the interior equilibrium , also let u and v denote û and v, respectively.System (3) is transformed into the following equivalent system: www.mii.lt/NADenote u •), and U = (u 1 , u 2 ) T .Then system (4) can be rewritten as an abstract differential equation in the phase space C τ , where The linearized equation of system (4) at the origin has the form: and its characteristic equation from [35] is where • represents a time delay variable.It is well known that the operator u → ∆u with ∂ ν u = 0 at 0 and lπ has eigenvalues −n 2 /l 2 , n = 0, 1, 2, . . ., with corresponding eigenfunctions cos(nπx/l).Let be an eigenfunction for ∆ + L with eigenvalue λ, see also [36].Hence, Eq. ( 7) can be written as where Assume that hypotheses (A1) and (A3) hold.Then λ = 0 is not a root of Eq. (5 n ) for any n ∈ N 0 .
From the result of Ruan and Wei [37], we know that the sum of the orders of the zeros of Eq. (5 n ) in the open right half plane can change only if a zero appears on or crosses the imaginary axis, as τ varies.Now, we focus on the conditions under which the cases above occur.Denote Lemma 3. Eq. (5 n ) has a pair of purely imaginary roots ±iω n (0 where Separating the real and imaginary parts, we have which is equivalent to Setting z = ω 2 , Eq. ( 10) is transformed into where From hypothesis (A3), B n + C > 0, For any u * ∈ (0, K], Eq. ( 11) has no positive roots for n > N 1 , and 0 ≤ n ≤ N 1 is the necessary condition of Eq. ( 11) having positive roots.For 0 ≤ n ≤ N 1 , a unique positive root z n of Eq. ( 11) is and is the imaginary part of the purely imaginary root.At Eq. (5 n ) has a pair of purely imaginary roots ±iω n (0 ≤ n ≤ N 1 ).Completing the proof.
Remark 2. Notice that Eq. (5 0 ) is the characteristic equation of the linearization of the system (2) at the positive equilibrium, which has been considered in [27].From the results in [27] listed above, we know when τ = τ k , system (2) has simple imaginary roots ±iσ + , which are just coincide with τ k 0 and ±iω 0 .That is, when n = 0, the conclusions of Lemma 2.3 are also the results in [27].
From Eq. ( 8) and Eq. ( 9), we have τ j n < τ j+1 n for any 0 ≤ n ≤ N 1 , j ∈ N 0 and the following lemma point out the increasing proposition of τ j n in n.Lemma 4. Let τ j n be defined as Eq. ( 8) and Eq.(9).Then for any 0 ≤ n ≤ N 1 , j ∈ N 0 , τ j n < τ j n+1 .Proof.From Eq. ( 8) and Eq. ( 9), we have , where Simple computation shows that ) is strictly increasing in n.Hence, τ j n is strictly increasing in n.
From the above proposition, we have Notice that what possibility will occur is τ j m = τ k n for some m > n, j < k.In this paper, we do not consider this case.In other words, we consider Then we have the following transversality condition.
Proof.From Eq. (5 n ), we have Hence, Substituting τ j n into the above equation, we obtain where cos(ω n τ From the above analysis, we have the following theorem. Theorem 3. Suppose that hypotheses (A1) and (A3) hold.For system (3), the following statements are true: are Hopf bifurcation values of system (3) and the bifurcating periodic solutions are all spatially homogeneous, which coincides with the periodic solution of the corresponding ODE system.And when τ ∈ D \ {τ k 0 : k ∈ N 0 }, system (3) also undergoes a Hopf bifurcation of spatially inhomogeneous periodic orbits.

Direction and stability of spatial Hopf bifurcation
In this section, we shall study the directions, stability and the period of bifurcating periodic solutions by applying the normal formal theory and the center manifold theorem of partial functional differential equations presented in [34,35].For fixed j ∈ N 0 , 0 ≤ n ≤ N 1 , we denote τ = τ j n .Let ũ(t, x) and ṽ(t, x) denote u(τ t, x) and v(τ t, x), respectively.And drop the tilde as a matter of convenience.Then system (4) can be transformed into Then Eq. ( 12) can be rewritten in an abstract form in the phase space where Consider the linear equation According to results in Section 2, we know the origin (0, 0) is an equilibrium of system (12), and Λ n := {iω n τ , −iω n τ } are characteristic values of the system (12) of ordinary functional differential equation www.mii.lt/NABy the Riesz representation theorem, there exists a 2 × 2 matrix function η n (θ, τ ), −1 ≤ θ ≤ 0, whose entries are of bounded variation such that In fact, we can choose Let A(τ ) denote the infinitesimal generators of the semigroup induced by the solutions of Eq. ( 15) and A * be the formal adjoint of A(τ ) under the bilinear pairing [36,39]).A(τ ) has a pair of simple purely imaginary eigenvalues ±iω n τ , and they are also eigenvalues of A * .Let P and P * be the center subspace, that is, the generalized eigenspace of A(τ ) and A * associated with Λ n , respectively.Then P * is the adjoint space of P and dim P = dim P * = 2.
From (20), we have By [35], W (z, z) satisfies where Hence, That is, 22), we know that for θ ∈ [−1, 0), Therefore, by (23), for −1 ≤ θ < 0, By the definition of A τ and ( 24), we have That is, where Using the definition of A τ and ( 24), we get that for −1 ≤ θ < 0, we get for n = 0, 1, 2, . . ., where • represents a time delay variable.That is where Similarly, from (24), we get Similar to the computing method of W 20 , we have where Thus we can compute the following quantities which determine the direction and stability of bifurcating periodic orbits: and we have Theorem 4. For any critical value τ j n , we have: (i) µ 2 determines the directions of the Hopf bifurcation: if µ 2 > 0 (< 0), then the direction of the Hopf bifurcation is forward (backward), that is, the bifurcating periodic solutions exists for τ > τ j n (τ < τ j n ); (ii) β 2 determines the stability of the bifurcating periodic solutions on the center manifold: if β 2 < 0 (> 0), then the bifurcating periodic solutions are orbitally asymptotically stable (unstable); (iii) T 2 determines the period of the bifurcating periodic solutions: if T 2 > 0 (< 0), then the period increases (decreases).

Optimal control strategy
In this section, applying the method and notations of Fister [38], we consider optimal control strategy for the system (3) with τ = 0.For the cause of prey conservation and economic benefit, we add predator harvesting term to the optimal control system and obtain the following system u(x, 0) = u 0 (x), v(x, 0) = v 0 (x), x ∈ Ω, where Ω is smooth bounded domain in R n , the functions E 1 and E 2 are controls that represent harvesting a proportion of the population.Our class of admissible control is Our payoff function is where K 1 E 1 u, K 2 E 2 v represent the revenue of harvesting and M 1 E 2 1 , M 2 E 2 2 denote the cost of the controls.We want to maximize the functional over the admissible class of controls, that is, there exists E * 1 and E * 2 such that J(E * 1 , E * 2 ) = max (E1,E2)∈E J(E 1 , E 2 ).For any (E 1 , E 2 ) ∈ E , we assume that there exists a pair of solution (u, v) of the system (25).Applying the method of Fister [38] and the estimates: www.mii.lt/NASimilarly, we rewrite the φ 1 , φ 2 PDE systems as We define the adjoint PDE system as , where L * p q = L * 1 p L * 2 q + M T p q , and L * 1 p = −p t − ∆p, L * 2 q = −q t − ∆q.For the adjoint system, we have the appropriate bounded conditions, namely, zero Neumann conditions and the transversality conditions p(x, T ) = 0, and q(x, T ) = 0 for x ∈ Ω.
then system (3) has two boundary equilibrium points E 1 and E 2 ; z),