Stability and Hopf Bifurcation of Delayed Predator-Prey System Incorporating Harvesting

Abstract

A kind of delayed predator-prey system with harvesting is considered in this paper. The influence of harvesting and delay is investigated. Our results show that Hopf bifurcations occur as the delay τ passes through critical values. By using of normal form theory and center manifold theorem, the direction of Hopf bifurcation and the stability of the bifurcating periodic solutions are obtained. Finally, numerical simulations are given to support our theoretical predictions.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Chun-Lei Tang

College of Mathematics and Computer Science, Fuzhou University, Fuzhou, Fujian 350116, China

Received 28 September 2013; Accepted 6 January 2014; 6 March 2014

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.

1. Introduction

The classical predator-prey systems have been extensively investigated in recent years, and they will continue to be one of the dominant themes in the future due to their universal existence and importance. Many biological phenomena are always described by differential equations, difference equations, and other type equations. In general, delay differential equations exhibit more complicated dynamical behaviors than ordinary ones; for example, the delay can induce the loss of stability, various oscillations, and periodic solutions. The dynamical behaviors of delay differential equations, stability, bifurcation and chaos, and so forth have been paid much attention by many researchers. Especially, the direction and stability of Hopf bifurcation to delay differential equations have been investigated extensively in recent work (see [1-7] and references therein).

After the classical predator-prey model was first proposed and discussed by May in [8], there were some similar topics, regarding persistence, local and global stabilities of equilibria, and other dynamical behaviors (see [5, 9, 10] and references therein). Recently, Song and Wei in [7] had considered a delayed predator-prey system as follows: [figure omitted; refer to PDF] where x ( t ) and y ( t ) were the densities of prey species and predator species at time t , respectively. The local Hopf bifurcation and the existence of the periodic solution bifurcating of system (1) was investigated in [7]. When selective harvesting was put into the predator-prey model similar to (1), Kar [11] studied two predator-prey models with selective harvesting; that is, in the first model, selective harvesting of predator species: [figure omitted; refer to PDF] and, in the second model, selective harvesting of prey species: [figure omitted; refer to PDF] had been considered by incorporating time delay on the harvesting term. They found that the delay for selective harvesting could induce the switching of stability and Hopf bifurcation occurred at τ = _{τ 0} .

Recently, Kar and Ghorai [9] had investigated a predator-prey model with harvesting: [figure omitted; refer to PDF] They obtained the local stability, global stability, influence of the harvesting, direction of Hopf bifurcation and the stability to system (4). Motivated by models (1)-(4), we will consider a predator-prey system with delay incorporating harvests to predator and prey: [figure omitted; refer to PDF] where x ( t ) and y ( t ) represent the population densities of prey species and predator species, respectively, at time t ; a , b , _{h 1} , _{h 2} , _{k 1} , and _{k 2} are model parameters assuming only positive values; _{k 1} measures the scale whose environment provides protection to prey x ; _{k 2} denotes the scale whose environment provides protection to predator y ; τ means the period of pregnancy; x ( t - τ ) represents the number of prey species which was born at time t - τ and still survived at time t ; _{h 1} and _{h 2} represent the coefficients of prey species and predator species, respectively. We always assume that 0 ...4; _{h 1} ...4; _{h 2} < 1 in this paper.

The organization of the paper is as follows. The stability of the positive equilibrium and the existence of the Hopf bifurcation are discussed in Section 2. The effect of harvesting to prey species and predator species is investigated in Section 3. The direction of Hopf bifurcation and stability of the corresponding periodic solution are obtained in Section 4. Numerical simulations are carried out to illustrate our results in Section 5.

2. Stability of Positive Equilibrium and Hopf Bifurcation

By simple computation, if _{k 1} + _{k 2} a ( _{h 2} - 1 ) > 0 holds, system (5) admits a unique positive equilibrium ^{E *} ( ^{x *} , ^{y *} ) : [figure omitted; refer to PDF] Let _{x 1} = x - ^{x *} , _{x 2} = y - ^{y *} , and then we get the linear system of (5): [figure omitted; refer to PDF] where _{a 11} = ^{x *} / ( _{k 1} - a ^{y *} ) , _{a 12} = a ^{x * 2} / ( _{k 1} - a ^{y * ) 2} , _{a 21} = b ^{y * 2} / ( _{k 2} + b ^{x * ) 2} , _{a 22} = ^{y *} / ( _{k 2} + b ^{x *} ) . From linear system (5) the characteristic equation is as follows: [figure omitted; refer to PDF] Roots of system (8) imply the stability of the equilibrium ^{E *} and Hopf bifurcation of system (5). Obviously, λ = 0 is not a root of system (8). For τ = 0 , system (8) becomes [figure omitted; refer to PDF] It is obvious that the root of system (9) has negative real part. Now, for τ > 0 , if λ = i ω ( ω > 0 ) is a root of (8), then we have [figure omitted; refer to PDF] Furthermore, [figure omitted; refer to PDF] which lead to polynomial equation [figure omitted; refer to PDF] It is easy to see that (12) has one positive root [figure omitted; refer to PDF] where Δ = ( ^{a 11 2} - ^{a 22 2 ) 2} + 4 ( _{a 11 a 22} + _{a 12 a 21} ) . By (11), one gets that [figure omitted; refer to PDF] Let [figure omitted; refer to PDF] be a pair of purely imaginary roots of (8), such that [figure omitted; refer to PDF] Next, we will prove λ ( _{τ j} ) meets the transversality conditions; taking the derivative of system (8) with respect to τ , one derives that [figure omitted; refer to PDF] which, together with (11), leads to [figure omitted; refer to PDF] So, we have [figure omitted; refer to PDF] Thus, we can obtain the following lemma.

Lemma 1.

If _{k 1} + _{k 2} a ( _{h 2} - 1 ) > 0 holds, then the following results are true:

(i) when τ = 0 , the positive equilibrium of ^{E *} of system (5) is locally asymptotically stable;

(ii) when 0 < τ < _{τ 0} , the positive equilibrium of ^{E *} of system (5) is locally asymptotically stable, and ^{E *} is unstable when τ > _{τ 0} , where ω , _{τ j} ( j = 0,1 , ... ) can be defined in (13), (14).

3. The Influence of Harvesting

Next, we will discuss the influence of the harvesting on system (5).

Case 1 (only predator species is harvested).

For _{h 1} = 0 , and the positive equilibrium of system (5) changes to ^{E 1 *} ( ^{x 1 *} , ^{y 1 *} ) , where [figure omitted; refer to PDF] it is obvious that ^{y 1 *} > 0 and ^{x 1 *} > 0 if and only if _{k 1} + _{k 2} a ( _{h 2} - 1 ) > 0 . Obviously, ^{x 1 *} and ^{y 1 *} are the continuous differentiable functions with respect to _{h 2} ; then, we have [figure omitted; refer to PDF]

Theorem 2.

If _{k 1} + _{k 2} a ( _{h 2} - 1 ) > 0 holds, then ^{x 1 *} is the monotonic increasing function of _{h 2} , ^{y 1 *} is the monotonic decreasing function of _{h 2} ; that is, when _{h 2} increases, the density of prey species will increase, the density of predator species will decrease.

Case 2 (only prey species is harvested).

For _{h 2} = 0 , and the positive equilibrium of system (5) changes to ^{E 2 *} ( ^{x 2 *} , ^{y 2 *} ) , where [figure omitted; refer to PDF] it is obvious that ^{y 2 *} > 0 and ^{x 2 *} > 0 if and only if _{k 1} - _{k 2} a > 0 . Obviously, ^{x 2 *} and ^{y 2 *} are the continuous differentiable functions with respect to _{h 1} ; then, one get that [figure omitted; refer to PDF]

Theorem 3.

If _{k 1} - _{k 2} a > 0 holds, then ^{x 2 *} and ^{y 2 *} are the monotonic decreasing functions of _{h 1} ; that is, if _{h 1} increases, then the density of prey species and predator species will decrease; on the contrary, if _{h 1} decreases, the density of prey species and predator species will increase.

Case 3 (predator species and prey species are harvested simultaneously).

For _{h 1 h 2} ...0; 0 , the mixed derivative of ^{x *} and ^{y *} are given by [figure omitted; refer to PDF]

Theorem 4.

If _{k 1} + _{k 2} a ( _{h 2} - 1 ) > 0 is valid, then the densities of prey species and predator species will both decrease when harvesting rate _{h 1} increases; on the contrary, the density of prey species will increase and predator species will decrease when harvesting rate _{h 2} increases.

4. Direction and Stability of Hopf Bifurcation

Motivated by the ideas of Hassard et al. [12], by applying the normal form theory and the center manifold theorem, the properties of the Hopf bifurcation at the critical value τ = _{τ j} are derived in this section.

Let t = s τ , _{x i} ( s τ ) = _{x ^ i} ( s ) , i = 1,2 , τ = _{τ 0} + μ , μ ∈ R ; _{τ 0} is defined by (14), we still denote _{x ^ i} ( s ) = _{u i} ( s ) and s = t , then system (5) is transformed into functional differential equations in C ( [ - 1,0 ] , ^{R 2} ) as [figure omitted; refer to PDF] where u ( t ) = ( _{u 1} ( t ) , _{u 2} ( t ) ^{) T} ∈ ^{R 2} , _{u t} ( θ ) = u ( t + θ ) , θ ∈ [ - 1,0 ] , and _{L μ} : C ( [ - 1,0 ] ; ^{R 2} ) [arrow right] R , f : R × C ( [ - 1,0 ] ; ^{R 2} ) [arrow right] R are given by [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

By Riesz representation theorem, there exists a function η ( θ , μ ) of bounded variation for θ ∈ [ - 1,0 ] , such that [figure omitted; refer to PDF] We choose [figure omitted; refer to PDF] where δ is the Dirac delta function. For [varphi] ∈ ^{C 1} ( [ - 1,0 ] , ^{R 2} ) , we define [figure omitted; refer to PDF] Then, system (25) can be transformed into an operator differential equation of the form [figure omitted; refer to PDF] where _{u t} ( θ ) = u ( t + θ ) , for θ ∈ [ - 1,0 ] . For ψ ∈ ^{C 1} ( [ 0,1 ] , ( ^{R 2 ) *} ) , we define [figure omitted; refer to PDF] and a bilinear inner product [figure omitted; refer to PDF] where η ( θ ) = η ( θ , 0 ) ; then, A ( 0 ) and ^{A *} are adjoint operators. Noting that ± i ω _{τ 0} are eigenvalues of A ( 0 ) , thus, they are also eigenvalues of ^{A *} . In order to calculate the eigenvector q ( θ ) of A ( 0 ) corresponding to the eigenvalue i ω _{τ 0} and p ( s ) of ^{A *} corresponding to the eigenvalue - i ω _{τ 0} , let q ( θ ) = ( 1 , α ^{) T e i ω _{τ 0} θ} be the eigenvector of A ( 0 ) corresponding to i ω _{τ 0} ; then, A ( 0 ) q ( θ ) = i ω _{τ 0} q ( θ ) .

By the definition of A ( 0 ) and (26), (30), then, [figure omitted; refer to PDF] Thus, we can get [figure omitted; refer to PDF] Similarly, let p ( s ) = D ( 1 , β ^{) T e i ω _{τ 0} s} be the eigenvector of ^{A *} corresponding to - i ω _{τ 0} ; by similar discussion, we get β = ( _{a 11} - i ω ) / _{a 21} ^{e i ω _{τ 0}} .

In view of standardization of p ( s ) and q ( θ ) ; that is, Y9; p ( s ) , q ( θ ) YA; = 1 , we have [figure omitted; refer to PDF] Thus, choose D = [ 1 + β α ¯ + _{τ 0} β ^{e i ω _{τ 0}} ( _{a 21} - α ¯ _{a 22} ) ^{] - 1} . Next, we will quote the same notation (see [13]), we first compute the coordinates to describe the center manifold _{C 0} at μ = 0 . Define [figure omitted; refer to PDF] On the center manifold _{C 0} , we have [figure omitted; refer to PDF] z and z ¯ are local coordinates for center manifold _{C 0} in the direction p and p ¯ ; noting that W is real if _{u t} is real, we only consider real solution _{u t} ∈ _{C 0} of (25). Since μ = 0 , then we have [figure omitted; refer to PDF] We rewrite this equation as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Noting _{u t} ( θ ) = ( _{[varphi] 1} ( θ ) , _{[varphi] 2} ( θ ) ^{) T} = W ( t , θ ) + z q ( θ ) + z ¯ q ¯ ( θ ) and q ( θ ) = ( 1 , α ^{) T e i ω _{τ 0} θ} , we have [figure omitted; refer to PDF] From (27), (42), we obtain that [figure omitted; refer to PDF] Because _{g 21} contains _{W 20} and _{W 11} , from (32) and (38), we have [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Substituting the corresponding series into (45) and comparing the coefficients, we have [figure omitted; refer to PDF] From (45), we know that for θ ∈ [ - 1,0 ) , we have [figure omitted; refer to PDF] Comparing the coefficient with (46) yields that for θ ∈ [ - 1,0 ) [figure omitted; refer to PDF] From (47), (49) and the definition of A , it follows that [figure omitted; refer to PDF] taking notice of q ( θ ) = ( 1 , α ^{) T e i ω _{τ 0} θ} ; hence, [figure omitted; refer to PDF] where _{E 1} = ( ^{E 1 ( 1 )} , ^{E 1 ( 2 )} ) ∈ ^{R 2} is a constant vector. By the similar way, we have [figure omitted; refer to PDF] where _{E 2} = ( ^{E 2 ( 1 )} , ^{E 2 ( 2 )} ) ∈ ^{R 2} is a constant vector.

Next, computing _{E 1} and _{E 2} , from the definition of A and (47), one then obtains [figure omitted; refer to PDF] where η ( θ ) = η ( 0 , θ ) . Furthermore, we have [figure omitted; refer to PDF] Substituting (52) and (56) into (54) and noting that [figure omitted; refer to PDF] it implies that [figure omitted; refer to PDF] Namely, [figure omitted; refer to PDF] Then it yields that [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Similarly, we get [figure omitted; refer to PDF] where [figure omitted; refer to PDF] Through simple computation, we determine _{W 20} , _{W 11} from (52) and (53); further, we can determine _{g 21} . Therefore, _{g i j} in (44) can be expressed by the parameter and delay; hence, [figure omitted; refer to PDF] which determine the qualities of bifurcation periodic solution of the critical value _{τ 0} .

Theorem 5.

(i) _{μ 2} determines the direction of Hopf bifurcation: if _{μ 2} > 0 ( < 0 ), then Hopf bifurcation is supercritical (subcritical), and the bifurcating periodic solutions exist for τ > _{τ 0} ( τ < _{τ 0} ).

(ii) ζ determines the stability of the bifurcating periodic solutions: the bifurcating periodic solutions are stable (unstable) if ζ < 0 ( ζ > 0 ). T determines the period of the bifurcating periodic solution: the period increases (decrease) if T > 0 ( T < 0 ).

5. Numerical Simulations

In this section, we consider a delayed predator-prey system with harvesting as follows: [figure omitted; refer to PDF] Because ( _{H 1} ) holds, from (14), we obtain that [figure omitted; refer to PDF] The unique positive equilibrium is ^{E *} = ( 2.907,3.436 ) .

If _{h 1} = 0.4 , when _{h 2} decreases, then prey species decreases and predator species increases (see Figure 1); when _{h 2} increases, prey species increases and predator species decreases (see Figure 2); If _{h 2} = 0.5 , when the values of harvesting _{h 1} decreases, then both predator species and prey species will increase (see Figure 3); on the other hand, when _{h 1} increases, then both predator species and prey species will decrease (see Figure 4).

Figure 1: When _{h 1} = 0.4 and _{h 2} decreases, prey species x decreases and predator species y increases.

[figure omitted; refer to PDF]

Figure 2: When _{h 1} = 0.4 and _{h 2} increases, prey species x increases and predator species y decreases.

[figure omitted; refer to PDF]

Figure 3: When _{h 2} = 0.5 and _{h 1} decreases, prey species x and predator species y increase.

[figure omitted; refer to PDF]

Figure 4: When _{h 2} = 0.5 and _{h 1} increases, prey species x and predator species y decrease.

[figure omitted; refer to PDF]

When parameter τ is little bigger than the critical value _{τ 0} , system (5) will become unstable and predator species and prey species can coexist; when τ increases much more, prey species will go to extinct (see Figure 5). Moreover, from Figure 6, we can see that system (5) is unstable when τ passes through the critical value _{τ 0} . By controlling the harvesting rates _{h 1} and _{h 2} , respectively, the stability of positive equilibrium to system (5) can been changed. Similarly, when τ < _{τ 0} , system (5) is stable; if we decrease the harvesting rate _{h 2} , then the stable system becomes unstable one (see Figure 7).

Figure 5: When τ = 2.9 > _{τ 0} [=, single dot above] 2.8015 , prey species x and predator species y coexist; when τ = 5 > _{τ 0} [=, single dot above] 2.8015 , prey species x goes to extinct.

[figure omitted; refer to PDF]

Figure 6: When τ = 2.9 > _{τ 0} [approximate] 2.8015 and _{h 1} increases, prey species x and predator species y become stable; when _{h 2} increases, prey species x and predator species y also become stable.

[figure omitted; refer to PDF]

Figure 7: When τ = 2.7 < _{τ 0} [approximate] 2.8015 and _{h 2} decreases, prey species x and predator species y become unstable.

[figure omitted; refer to PDF]

Since _{μ 2} < 0 , ζ < 0 , Hopf bifurcation is subcritical and the positive equilibrium ^{E *} is asymptotically stable for 0 < τ < _{τ 0} (see Figure 8); when τ > _{τ 0} , ^{E *} loses its stability and Hopf bifurcation occurs; that is, a family of periodic solutions bifurcate from ^{E *} (see Figure 9).

Figure 8: When τ = 2.7 < _{τ 0} [=, single dot above] 2.8015 , the positive equilibrium ^{E *} of system (5) is asymptotically stable.

[figure omitted; refer to PDF]

Figure 9: When τ = 2.9 > _{τ 0} [=, single dot above] 2.8015 , the positive equilibrium ^{E *} of system (5) loses its stability and a Hopf bifurcations occurs.

[figure omitted; refer to PDF]

As discussed, our results show that the delay τ affects the stability of system (5) and harvesting rates _{h 1} and _{h 2} also affect the stability of system (5).

6. Conclusion

In our model, the harvesting term has been introduced into the model (5); by applying the normal form theorem and the center manifold theorem, we investigate the dynamical behaviors of the delayed predator-prey model with harvesting term and obtain the influence of harvesting term on the prey species and predator species. Further, we prove that the influence of the harvesting rates _{h 1} and _{h 2} to the stability of system (5), by controlling harvesting rates _{h 1} and _{h 2} of prey species and predator species, which makes the unstable (stable) system become stable (unstable).

Our results show that Hopf bifurcations occur as the delay τ passes through critical values _{τ 0} [approximate] 2.8015 . When τ < _{τ 0} , the positive equilibrium ^{E *} of system (5) is asymptotically stable; when τ > _{τ 0} , the positive equilibrium ^{E *} of system (5) loses its stability and Hopf bifurcations occur.

Acknowledgments

This work is supported by National Natural Science Foundation of China (no. 11201075) and Natural Science Foundation of Fujian Province of China (no. 2010J01005).

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] K. Gopalsamy Stability and Oscillations in Delay Differential Equations of Population Dynamics , vol. 74, of Mathematics and its Applications, Kluwer Academic Publishers Group, Dordrecht, The Netherlands, 1992.

[2] D. Xiao, S. Ruan, "Multiple bifurcations in a delayed predator-prey system with nonmonotonic functional response," Journal of Differential Equations , vol. 176, no. 2, pp. 494-510, 2001.

[3] J.-F. Zhang, "Bifurcation analysis of a modified Holling-Tanner predator-prey model with time delay," Applied Mathematical Modelling , vol. 36, no. 3, pp. 1219-1231, 2012.

[4] S. Yuan, Y. Song, "Stability and Hopf bifurcations in a delayed Leslie-Gower predator-prey system," Journal of Mathematical Analysis and Applications , vol. 355, no. 1, pp. 82-100, 2009.

[5] N. Bairagi, D. Jana, "On the stability and Hopf bifurcation of a delay-induced predator-prey system with habitat complexity," Applied Mathematical Modelling , vol. 35, no. 7, pp. 3255-3267, 2011.

[6] C. Çelik, "The stability and Hopf bifurcation for a predator-prey system with time delay," Chaos, Solitons and Fractals , vol. 37, no. 1, pp. 87-99, 2008.

[7] Y. Song, J. Wei, "Local Hopf bifurcation and global periodic solutions in a delayed predator-prey system," Journal of Mathematical Analysis and Applications , vol. 301, no. 1, pp. 1-21, 2005.

[8] R. M. May, "Time delay versus stability in population model with two or three trophic levels," Ecology , vol. 54, pp. 315-325, 1973.

[9] T. K. Kar, A. Ghorai, "Dynamic behaviour of a delayed predator-prey model with harvesting," Applied Mathematics and Computation , vol. 217, no. 22, pp. 9085-9104, 2011.

[10] S. B. Hsu, T. W. Huang, "Global stability for a class of predator-prey systems," SIAM Journal on Applied Mathematics , vol. 55, no. 3, pp. 763-783, 1995.

[11] T. K. Kar, "Selective harvesting in a prey-predator fishery with time delay," Mathematical and Computer Modelling , vol. 38, no. 3-4, pp. 449-458, 2003.

[12] B. D. Hassard, N. D. Kazarinoff, Y. H. Wan Theory and Applications of Hopf Bifurcation , vol. 41, Cambridge University Press, Cambridge, UK, 1981.

[13] K. Gopalsamy, "Harmless delays in model systems," Bulletin of Mathematical Biology , vol. 45, no. 3, pp. 295-309, 1983.