- Email: [email protected]
Bifurcation analysis in a predator–prey system with stage-structure and harvesting
Journal of the Franklin Institute 347 (2010) 1097–1113 www.elsevier.com/locate/jfranklin
Bifurcation analysis in a predator–prey system with stage-structure and harvesting$ Ying Qua,, Junjie Weia,b a Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China Department of Mathematics, Harbin Institute of Technology (Weihai), Weihai, Shandong 264209, China
b
Received 18 November 2009; received in revised form 28 February 2010; accepted 29 March 2010
Abstract In this paper, we consider a predator–prey model with stage-structure and harvesting. This model is the same as the one developed by Kar and Pahari (2007) [9], but we make bifurcation analysis more general than their work. In particular, using the approach of Beretta and Kuang (2002) [4], we show that the positive steady state can be destabilized through a Hopf bifurcation. We also investigate the stability and direction of periodic solutions bifurcating from Hopf bifurcation by using the normal form theory and the center manifold theorem presented in Hassard et al. (1981) [8]. Numerical simulations are then carried out as supporting evidences of our analytical results. & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. Keywords: Predator–prey system; Harvesting; Time-delay; Stability; Hopf bifurcation
1. Introduction The dynamic of predator–prey model has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to the fact that predator– prey interaction is the fundamental structure in population dynamics. In most of the related studies, stage-structure is modeled by the introduction of a time delay, which leads to retarded functional differential equations. In general, equations of $
This research is supported by the NNSF of China (No. 10771045), the Program of Excellent Team in HIT and the Natural Science Foundation of Heilongjiang Province (A200806). Corresponding author. E-mail addresses: [email protected], [email protected] (Y. Qu). 0016-0032/$32.00 & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2010.03.017
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1098
this kind often exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the involved populations to fluctuate (see [2,3,9–14] and the references therein). Harvesting has a strong impact on the dynamic evolution of a population. To a certain extent, it can control the long-term stationary density of population efficiently. However, it can also lead to the incorporation of a positive extinction probability and, therefore, to potential extinction in finite time. The problems of predator–prey systems with harvesting have been discussed by many authors (see [5,13] for example). Recently, Kar and Pahari [9] proposed a predator–prey-type fishery model with Beddington–DeAngelis functional response and selective harvesting of predator species: x axy ; x_ ¼ rx 1 k a þ bx þ cy y_ ¼ b0 e
gt
bx yðttÞ d0 yqEy; a þ bx þ cy
ð1:1Þ
where x(t) and y(t) denote the population densities of prey and predator at time t, respectively. All coefficients are positive constants. Specially, r and k represent the intrinsic growth rate and carrying capacity of prey in the absence of predation, respectively. The predator consumes the prey with functional response axy=ða þ bx þ cyÞ, known as Beddington–DeAngelis response and contributes to its growth rate bxy=ða þ bx þ cyÞ. b0 is a constant birth rate and d0 is the maximum death rate in the absence of food. All the metabolic energy a predator obtains from its food goes into physiological maintenance, rather than into enhancing reproductive effort. In addition, they assumed that juveniles suffer a mortality rate of g (the through-stage death rate) and take t units of time to mature. Therefore, the term b0 egt yðttÞ represents the population of juveniles who were born at time tt and survive at time t (with the immature death rate g), and thus represents the transformation of juveniles to adults. Moreover, the catch rate function qEy is based on the catch-per-unit-effort (CPUE) hypothesis in [6]. Here, q is the catchability coefficient of the predator species and E is the harvesting effort. In [9], when t ¼ 0, the authors established the global asymptotical stability for the interior equilibrium provided the local asymptotical stability holds. When t40, the authors considered the special case where juveniles mortality rate is g ¼ 0. Under such a hypothesis, the dependence of coefficients on t in the characteristic equation associated with Eq. (1.1) vanishes, because the critical item egt leading to the dependence becomes the constant 1. For the case ga0, they stated that the analytical computations were difficult. Therefore, they only gave the numerical simulation to exhibit the phenomena that stability switches occurred when t is varied. As is well known, for some species (e.g., Antarctic whale and seal populations), delay (the length of time to maturity) is a function of the amount of food (mostly krill), therefore, taking t as parameter is biologically meaningful. From this point of view and as a significant extension of [9], we will give a detailed analysis for the case ga0 which Kar and Pahari failed to examine. Ultimately, the investigations will show the sensitivity of the model dynamics on maturation time delay. We remark that, besides the time delay parameter t, the other parameters in our model are also important biological factors. However, when discussing the dynamical behavior of our system, the key step is to analyze the characteristic equation (2.5) (see below), which is
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1099
a transcendental equation with coefficients depending on the parameters q, g, E,y, in a manner that is quite complicated. Therefore in the present study, we choose to add some simulation results to illustrate the influence of the harvesting effort E on the dynamical behavior of our system (see the last section) but leave out the corresponding analyses to future investigations. This paper is organized as follows. Section 2 gives the analysis of existence and uniqueness of positive equilibrium as well as its stability, by studying the characteristic equation associated with Eq. (1.1) which takes the form of a second degree exponential polynomial with delay-dependent coefficients. Section 3 is devoted to the investigation of the stability and direction of periodic solutions bifurcating from Hopf bifurcation by using the normal form theory and the center manifold theorem presented in Hassard et al. [8]. Finally, some numerical simulations will be given to illustrate our theoretical results in Section 4.
2. Interior equilibrium and its stability We rewrite Eq. (1.1) here for the sake of convenience: x axy ; x_ ¼ rx 1 k a þ bx þ cy y_ ¼ b0 egt yðttÞ d0
bx yqEy; a þ bx þ cy
ð2:1Þ
for tZ0 with the initial values yðyÞ ¼ fðyÞZ0; y 2 ½t;0 and xð0Þ40; yð0Þ40:
2.1. Existence and uniqueness of interior equilibrium We consider the biologically meaningful interior equilibrium E , which satisfies x ay r 1 ; ¼ k a þ bx þ cy TðtÞ ¼
bx ; a þ bx þ cy
where TðtÞ ¼ d0 þ qEb0 egt : Note from Eq. (2.2) that cbr 2 cbr x þ bTðtÞb xTðtÞa ¼ 0; ka a y¼
br ðkxx2 Þ: kaTðtÞ
ð2:2Þ
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1100
For convenience, define bk b cr ; þ1 Tm ¼ min kb þ a b a and denote x7 ¼
pffiffiffiffi k ½ðaTðtÞb þ cbrabÞ7a D; 2cbr
y7 ¼
br ½kx7 ðx7 Þ2 kaTðtÞ
ð2:3Þ
with
cbr D9 TðtÞb þ b a
2 þ 4TðtÞa
cbr : ka
Proposition 2.1. Assume (A1) 0oTðtÞoTm . Then system (2.1) has a unique positive equilibrium E ¼ ðxþ ; yþ Þ. Proof. TðtÞ40 implies that xþ 404x . Accordingly, in order to verify yþ 40, we only need to examine xþ ok. By the expression of xþ given by Eq. (2.3), it follows that pffiffiffiffi cbr TðtÞb þ b: ð2:4Þ Do a This shows that we need cbr=aTðtÞb þ b40; which is equivalent to b cr TðtÞo þ1 : b a Hence, the inequality (2.4) is equivalent to TðtÞo
bk : kb þ a
This completes the proof.
&
Note that TðtÞ is a monotone increasing function of t such that d0 þ qEb0 ¼ Tð0ÞrTðtÞoTð1Þ ¼ d0 þ qE: Therefore, (A1) holds if 0
(A1 ) b0 od0 þ qErTm is satisfied. In what follows, we always assume (A1) holds and we also denote ðxþ ; yþ Þ by ðx ; y Þ.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1101
2.2. Stability analysis Linearizing Eq. (2.1) around ðx ; y Þ gives
2x ay ða þ cy Þ ax ða þ bx Þ _ ¼ r 1 x y; xðtÞ k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 by ða þ cy Þ bx ða þ bx Þ _yðtÞ ¼ b0 egt yðttÞ þ xþ d0 qE y; ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 and its characteristic equation is l2 þ AðtÞl þ BðtÞ þ ½CðtÞl þ DðtÞelt ¼ 0;
ð2:5Þ
where
2x ay ða þ cy Þbx ða þ bx Þ AðtÞ ¼ d0 þ qEr 1 þ ; k ða þ bx þ cy Þ2
ay ða þ cy Þ 2x bx ða þ bx Þ d ; r 1 þ qE BðtÞ ¼ ðd0 þ qEÞ 0 k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 CðtÞ ¼ b0 egt o0;
ay ða þ cy Þ 2x DðtÞ ¼ b0 egt : þ r 1 k ða þ bx þ cy Þ2
Note that these coefficients are all dependent on t. In addition, AðtÞ; BðtÞ; CðtÞ and DðtÞ are all continuous and differentiable functions from ½0; 1Þ to R. We now investigate the local asymptotic stability of ðx ; y Þ for system (2.1). It is well known that the trivial solution of system (2.1) is locally asymptotically stable if all roots l of the characteristic equation (2.5) satisfy ReðlÞo0: Eq. (2.5) takes the general form Pðl;tÞ þ Qðl;tÞelt ¼ 0; where
ð2:6Þ
Pðl;tÞ9l2 þ AðtÞl þ BðtÞ; Qðl;tÞ9CðtÞl þ DðtÞ:
ð2:7Þ
Let us first examine the case t ¼ 0. This case is of importance, because if the nontrivial positive steady-state of Eq. (2.1) is stable when t ¼ 0, then we are able to obtain the local stability for all nonnegative values of the delay, or to find a critical value which could destabilize the steady-state. When t ¼ 0, the characteristic equation (2.5) reduces to l2 þ ½Að0Þ þ Cð0Þl þ ½Bð0Þ þ Dð0Þ ¼ 0: Clearly, all roots of Eq. (2.6) with t ¼ 0 have negative real parts provided (H1) Að0Þ þ Cð0Þ40 and Bð0Þ þ Dð0Þ40: Hence, we get the following conclusion.
1102
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
Proposition 2.2 (Kar and Pahari [9, Theorem 2.1]). Assume (A1) and (H1) hold. Then the positive equilibrium E ¼ ðx ; y Þ of Eq. (2.1) is globally asymptotically stable if t ¼ 0. Now we assume t40 and regard it as a bifurcation parameter to obtain finer results on the stability of ðx ; y Þ. Note that Eq. (2.5) takes the form of a second degree exponential polynomial in l, with all the coefficients of P and Q depending on t. Thus, we use the method introduced by Beretta and Kuang [4], which gives the existence of purely imaginary roots of a characteristic equation with delay-dependent coefficients (see also [1,11]). In order to apply the criterion in [4], we need to verify the following properties for o40 and t 2 I with I defined in Eq. (2.8): (i) (ii) (iii) (iv) (v)
Pð0; tÞ þ Qð0; tÞa0; Pðio; tÞ þ Qðio; tÞa0; lim supjlj-1;RelZ0 jQðl; tÞ=Pðl; tÞjo1; F ðoÞ9jPðio; tÞj2 jQðio; tÞj2 has a finite number of zeros; each positive root oðtÞ of F ðoðtÞÞ ¼ 0 is continuous and differentiable in t whenever it exists.
Here Pðl; tÞ and Qðl; tÞ are defined by Eq. (2.7) and I9ftZ0 : there exists oðtÞ40 such that F ðoðtÞÞ ¼ 0g:
ð2:8Þ
Let F be defined as in (iv). By Eq. (2.7), we obtain F ðoÞ ¼ o4 þ a1 ðtÞo2 þ a2 ðtÞ
ð2:9Þ
with a1 ðtÞ ¼ A2 ðtÞ2BðtÞC 2 ðtÞ; a2 ðtÞ ¼ B2 ðtÞD2 ðtÞ: Before proceeding further, let us analyze the structure of I described by Eq. (2.8). Set qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 ðtÞ7 a21 ðtÞ4a2 ðtÞ : z7 ðtÞ ¼ 2 We need the following hypotheses: (H2) AðtÞ þ CðtÞa0 and BðtÞ þ DðtÞa0: (H3) 0oa2 ðtÞoa21 ðtÞ=4 and a1 ðtÞo0: (H4) a2 ðtÞr0: Proposition 2.3. Assume (A1) holds. pffiffiffiffiffiffiffiffiffiffiffi (i) If (H3) is satisfied, then F ðoÞ ¼ 0 has two different positive roots z7 ðtÞ denoted by o7 , respectively. pffiffiffiffiffiffiffiffiffiffiffi (ii) If (H4) is met, then F ðoÞ ¼ 0 has a unique positive root oþ ¼ zþ ðtÞ.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1103
Define I1 ¼ ftZ0 : ðA1Þ; ðH2Þ and ðH3Þ holdg; I2 ¼ ftZ0 : ðA1Þ; ðH2Þ and ðH4Þ holdg: Then I ¼ I1 [ I2 :
ð2:10Þ
Now, we verify the properties (i)–(v) for t 2 I. Firstly, (i) and (ii) are satisfied. Indeed, Pð0;tÞ þ Qð0;tÞ ¼ BðtÞ þ DðtÞa0 and Pðio;tÞ þ Qðio;tÞ ¼ ½BðtÞ þ DðtÞo2 þ io½AðtÞ þ CðtÞa0: Then, from Eq. (2.7) we know Qðl;tÞ ¼ 0: lim jlj-1 Pðl;tÞ Therefore (iii) follows. Finally, Eq. (2.9) implies (iv) and (v). Let l ¼ io (o40) be a root of Eq. (2.5). Substituting it into Eq. (2.5) and separating the real and imaginary parts yields o2 BðtÞ ¼ DðtÞcosot þ CðtÞosinot; AðtÞo ¼ DðtÞsinotCðtÞocosot; then we have sinot ¼
ðo2 BÞCo þ oAD o2 C 2 þ D2
and
cosot ¼
ðBo2 ÞD þ o2 AC : o2 C 2 þ D2
ð2:11Þ
Here the dependence of the coefficients on t is implicitly assumed. By the definitions of Pðl; tÞ and Qðl; tÞ in Eq. (2.7) and applying the property (ii), Eq. (2.11) can be written as
Pðio;tÞ Pðio;tÞ and cosot ¼ Re ; sinot ¼ Im Qðio;tÞ Qðio;tÞ which lead us to F ðoÞ ¼ 0:
ð2:12Þ
Noting that Proposition 2.3 has given the explicit expressions of oðtÞ that satisfies F ðoðtÞÞ ¼ 0 on t 2 I, we can define yðtÞ 2 ½0; 2p for t 2 I: sinyðtÞ ¼
ðo2 BÞCo þ oAD o2 C 2 þ D2
and
cosyðtÞ ¼
ðBo2 ÞD þ o2 AC : o2 C 2 þ D2
Such yðtÞ defined above is well and uniquely defined for all t 2 I (see [4]).
ð2:13Þ
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1104
Therefore, io with o ¼ oðt Þ40 is a purely imaginary root of Eq. (2.5) if and only if t is a zero of the Sn, where
Sn ðtÞ9t
yðtÞ þ 2np ; t2I oðtÞ
with n 2 N0 9f0; 1; 2; . . .g. The following result is due to Beretta and Kuang [4]. Lemma 2.4. The characteristic equation (2.5) has a pair of simple pure imaginary roots l ¼ 7ioðt Þ at t 2 I, provided Sn ðt Þ ¼ 0 for n 2 N0 . Moreover, if oðt Þ ¼ oþ ðt Þ, then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if dþ ðt Þ40, and crosses the imaginary axis from right to left if dþ ðt Þo0, where ) dReðlÞ dSn ðtÞ ; dþ ðt Þ ¼ Sign ¼ Sign dt l¼ioþ ðt Þ dt t¼t and if oðt Þ ¼ o ðt Þ, then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if d ðt Þ40, and crosses the imaginary axis from right to left if d ðt Þo0, where ) dReðlÞ dSn ðtÞ : d ðt Þ ¼ Sign ¼ Sign dt l¼io ðt Þ dt t¼t If I ¼ I2 , then only oþ is feasible. In this case, we can easily observe that Sn ð0Þr0 and Sn ðtÞ4Snþ1 ðtÞ for t 2 I and n 2 N0 . Therefore, Sn has no zeros in I for any n 2 N0 if S0 has no zeros in I. As a result, Eq. (2.1) has no scenario at ðx ; y Þ. If Sn ðtÞ has positive zeros, denoted by tjn , for some n 2 N0 , then without loss of generality, we may suppose dSn ðtjn Þ a0 dt
with Sn ðtjn Þ ¼ 0:
ð2:14Þ
Then for one thing, stability switches occur at the zeros of S0 ðtÞ, denoted by tj0 , if (H1) holds. For another, applying Proposition 2.2 and Hopf bifurcation theorem for functional differential equations (see Hale’s book [7]), we can conclude the existence of a Hopf bifurcation. Before stating the main theorem, denote tm ¼ minft 2 I : S0 ðtÞ ¼ 0g
and
tM ¼ maxft 2 I : S0 ðtÞ ¼ 0g:
Theorem 2.5. Assume (A1), (H1), (H2) and (H4) hold. (i) If S0 ðtÞ has no positive zeros in I, then ðx ; y Þ is locally asymptotically stable for all tZ0. (ii) If Sn ðtÞ has at least one positive zero in I and Eq. (2.14) is met, then ðx ; y Þ is locally asymptotically stable for t 2 ½0; tm Þ [ ðtM ; 1Þ and unstable for t 2 ðtm ; tM Þ. That is, stability switches of stability–instability–stability occur. Hopf bifurcation takes place when t ¼ tjn , n 2 N 0 . If I ¼ I1 , then both oþ and o are feasible for t 2 I. We have the following two sequences of functions on I: Snþ ðtÞ ¼ tðyþ ðtÞ þ 2npÞ=oþ ðtÞ
and
Sn ðtÞ ¼ tðy ðtÞ þ 2npÞ=o ðtÞ;
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1105
where y7 ðtÞ is the solution of Eq. (2.13) when o ¼ o7 . Similarly, it is obtained for t 2 I 7 ðtÞ with n 2 N0 . Furthermore, if S0þ ðtÞ4S0 ðtÞ, then that Sn7 ð0Þr0 and Sn7 ðtÞ4Snþ1 þ Sn ðtÞ4Sn ðtÞ for n 2 N0 . Other than the case above, stability switch may depend on all real roots of Snþ ðtÞ ¼ 0 and Sn ðtÞ ¼ 0. In addition, Hopf bifurcations can also appear at each root of Sn7 ðtÞ ¼ 0. Naturally, we can also obtain results similar to Theorem 2.5. Here, we omit the corresponding statements. If I ¼ I1 [ I2 with I1 a| and I2 a|, then we can discuss different cases on their own sets, respectively. The remaining discussion is the same as above.
3. The properties of Hopf bifurcation In the previous section, we obtained some sufficient conditions ensuring system (2.1) undergoes a Hopf bifurcation at E ¼ ðx ; y Þ. In this section, under these conditions such as Theorem 2.5(iii), we shall establish the explicit formula for determining the direction and stability of periodic solutions bifurcating from E at each zero t of Sn , by using normal form theory and center manifold argument presented in Hassard et al. [8] (see also [7,15,16]). By the change of variables xðtÞx /xðtÞ and yðtÞy /yðtÞ and by rescaling the time t/ðt=tÞ, system (2.1) can be written as ! ! ! _ xðtÞ xðtÞ xðt1Þ ¼ tMðtÞ þ tNðtÞ þ tF ðx;yÞ; ð3:1Þ _ yðtÞ yðtÞ yðt1Þ where
0 2x ay ða þ cy Þ r 1 B k ða þ bx þ cy Þ2 B MðtÞ ¼ B @ by ða þ cy Þ ða þ bx þ cy Þ2
NðtÞ ¼
F ðx; yÞ9
!
0
0
0
b0 egt
f ðx; yÞ lðx; yÞ
1 ax ða þ bx Þ C ða þ bx þ cy Þ2 C C; bx ða þ bx Þ A d qE 0 ða þ bx þ cy Þ2
!
;
0
1 rx2 axy B k a þ bx þ cy C B C ¼B C bxy @ A a þ bx þ cy 0 1 1 1 1 2 2 3 x þ f xy þ y þ x þ f f f 11 02 30 B 2 20 C 2 6 C ¼B @1 A 1 1 2 3 2 l20 x þ l11 xy þ l02 y þ l30 x þ 2 2 6
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1106
with fij ¼
@iþj f ðx ;y Þ @i x@j y
and
lij ¼
@iþj lðx ;y Þ ; i;j ¼ 1;2; . . . : @i x@j y
Let t ¼ t þ m. Thus m ¼ 0 is a Hopf bifurcation value of Eq. (3.1). For f ¼ colðf1 ; f2 Þ 2 Cð½1; 0; R2 Þ, define Lm ðfÞ ¼ ðt þ mÞ½Mðt þ mÞfð0Þ þ Nðt þ mÞfð1Þ: By the Riesz representation theorem, there exists a 2 2 matrix Zðy; mÞ; y 2 ½1; 0, whose elements are functions of bounded variation such that Z 0 Lm ðfÞ ¼ ½dZðy;mÞfðyÞ; f 2 Cð½1;0;R2 Þ: ð3:2Þ 1
In fact, we can choose Zðy; mÞ ¼ ðt þ mÞðMðt þ mÞdðyÞCðt þ mÞdðy þ 1ÞÞ; that is
8 > < ðt þ mÞBðt þ mÞ; Zðy; mÞ ¼ 0; > : ðt þ mÞCðt þ mÞ;
y ¼ 0; y 2 ð1; 0Þ; y ¼ 1:
Then Eq. (3.2) is satisfied. For f 2 C 1 ð½1; 0; R2 Þ, define the operator AðmÞ as 8 > < dfðyÞ; y 2 ½1; 0Þ; AðmÞfðyÞ ¼ R dy > : 0 ½dZðx; mÞfðxÞ; y ¼ 0 1 and
( RðmÞfðyÞ ¼
with
0; hðm;
y 2 ½1;0Þ; fÞ;
y¼0
0
1 1 1 1 2 2 3 B 2 f20 f1 ð0Þ þ f11 f1 ð0Þf2 ð0Þ þ 2 f02 f2 ð0Þ þ 6 f30 f1 ð0Þ þ C C: hðm; fÞ ¼ ðt þ mÞB @1 A 1 1 2 2 3 l20 f1 ð0Þ þ l11 f1 ð0Þf2 ð0Þ þ l02 f2 ð0Þ þ l30 f1 ð0Þ þ 2 2 6
Then system (2.1) is equivalent to the following operator equation: u_t ¼ AðmÞut þ RðmÞut ; where uðtÞ ¼ colðxðtÞ; yðtÞÞ and ut ¼ uðt þ yÞ for y 2 ½1; 0: For c 2 C 1 ð½0; 1; ðR2 Þ Þ, define 8 > < dcðsÞ; s 2 ð0; 1; A cðsÞ ¼ R ds > : 0 cðxÞdZðx;0Þ; s ¼ 0 1
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1107
and a bilinear form /cðsÞ;fðyÞS ¼ cð0Þfð0Þ
Z
0
Z
y
cðxyÞ dZðyÞfðxÞ dx; 1
x¼0
where ZðyÞ ¼ Zðy; 0Þ. Then Að0Þ and A are adjoint operators. From the discussion in Section 2, we know that 7io t are eigenvalues of Að0Þ and therefore they are also eigenvalues of A . Moreover, it is not difficult to verify that vectors qðyÞ ¼ colðq1 ; q2 Þeio t y ðy 2 ½1; 0Þ and q ðsÞ ¼ ð1=K Þðq1 ; q2 Þeio t s ðs 2 ½0; 1Þ are the eigenvectors of Að0Þ and A corresponding to the eigenvalues io t and io t , respectively, where q1 ¼ q1 ¼ 1; q2 ¼
rða þ bx þ cy Þ2 ðk2x Þaky ða þ cy Þ o ða þ bx þ cy Þ2 i; ax ða þ bx Þ kax ða þ bx Þ
q2 ¼
aky ða þ cy Þ þ rða þ bx þ cy Þ2 ð2x kÞ o ða þ bx þ cy Þ2 i: by ða þ cy Þ kby ða þ cy Þ
Let K ¼ 1 þ q2 q2 ð1 þ b0 t et
ðio þgÞ
Þ:
Then /q ðsÞ; qðyÞS ¼ 1 and /q ðsÞ; qðyÞS ¼ 0. Following the same algorithms and using a similar computation process as in Hassard et al. [8], we obtain the coefficients determining the important quantities: g20 ¼ g11 ¼
t ½q ðf20 q21 þ 2f11 q1 q2 þ f02 q22 Þ þ q 2 ðl20 q21 þ 2l11 q1 q2 þ l02 q22 Þ; K 1 t fq ðf20 jq1 j2 þ f11 ðq1 q 2 þ q 1 q2 Þ þ f02 jq2 j2 Þ K 1 þ q 2 ðl20 jq1 j2 þ l11 ðq1 q 2 þ q 1 q2 Þ þ l02 jq2 j2 Þg;
g02 ¼
t ½q ðf20 q 21 þ 2f11 q 1 q 2 þ f02 q 22 Þ þ q 2 ðl20 q 21 þ 2l11 q 1 q 2 þ l02 q 22 Þ; K 1
g21 ¼
t ð1Þ ð1Þ fq ½f20 ðW20 ð0Þq 1 þ 2W11 ð0Þq1 Þ K 1 ð1Þ ð1Þ ð2Þ ð2Þ þ f11 ðW20 ð0Þq 2 þ 2W11 ð0Þq2 þ 2W11 ð0Þq1 þ W20 ð0Þq 1 Þ ð2Þ ð2Þ ð0Þq 2 þ 2W11 ð0Þq2 Þ þ f30 q21 q 1 þ f21 ðq21 q 2 þ 2jq1 j2 q2 Þ þ f02 ðW20
þ f12 ðq22 q 1 þ 2jq2 j2 q1 Þ þ f03 q22 q 2 ð1Þ ð1Þ ð0Þq 1 þ 2W11 ð0Þq1 Þ þ q 2 ½l20 ðW20 ð1Þ ð1Þ ð2Þ ð2Þ ð0Þq 2 þ 2W11 ð0Þq2 þ 2W11 ð0Þq1 þ W20 ð0Þq 1 Þ þ l11 ðW20 ð2Þ ð2Þ þ l02 ðW20 ð0Þq 2 þ 2W11 ð0Þq2 Þ þ l30 q21 q 1 þ l21 ðq21 q 2 þ 2jq1 j2 q2 Þ
þ l12 ðq22 q 1 þ 2jq2 j2 q1 Þ þ l03 q22 q 2 g;
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1108
where W20 ðyÞ ¼
ig20 ig qðyÞ þ 02 qðyÞ þ E1 e2io t y ; o t 3o t
W11 ðyÞ ¼
ig11 ig qðyÞ þ 11 qðyÞ þ E2 o t o t
and
0 2x ay ða þ cy Þ r 1 2io B k ða þ bx þ cy Þ2 B E1 ¼ 2B by ða þ cy Þ @ ða þ bx þ cy Þ2 0 1 1 1 2 2 f f q þ f q q þ q 11 1 2 02 2 C B 2 20 1 2 C; B @1 A 1 l20 q21 þ l11 q1 q2 þ l02 q22 2 2
11 ax ða þ bx Þ ða þ bx þ cy Þ2 C C C A M
0 11 2x ay ða þ cy Þ ax ða þ bx Þ r 1 B C k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 B C E 2 ¼ B C by ða þ cy Þ bx ða þ bx Þ @ A d0 qE þ b0 egt 2 2 ða þ bx þ cy Þ ða þ bx þ cy Þ 0 1 1 1 2 2 f f jq j þ f ðq q þ q q Þ þ jq j Þ 20 1 11 1 2 02 2 2 1 B2 C 2 C B @1 A 1 2 2 l20 jq1 j þ l11 ðq1 q 2 þ q 1 q2 Þ þ l02 jq2 j Þ 2 2 with M¼
bx ða þ bx Þ d0 qE þ b0 egt2io t 2io : 2 ða þ bx þ cy Þ
Consequently, we can compute the following values: 2
i g21 2 jg02 j c1 ð0Þ ¼ þ ; g11 g20 2jg11 j 3 2 2o t m2 ¼
Reðc1 ð0ÞÞ ; 0 Reðl ðt ÞÞ
b2 ¼ 2 Reðc1 ð0ÞÞ; 0
T2 ¼
Imðc1 ð0ÞÞ þ m2 Imðl ðt ÞÞ ; o t
which determine the properties of bifurcating periodic solutions at the critical value t . To be concrete, m2 determines the direction of the Hopf bifurcation: if m2 40 (resp. o0), then Hopf bifurcations are forward (resp. backward); b2 determines the stability of bifurcating
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1109
periodic solutions: the bifurcating periodic solutions on the center manifold are stable (resp. unstable) if b2 o0 (resp. 40); T2 determines the periods of the bifurcating periodic solutions: the periods increase (resp. decrease) if T2 40 (resp. o0). Combining these properties with discussions in Section 2, we have: Theorem 3.1. Under the conditions given in Theorem 2.5(iii), Hopf bifurcation at each t is 0 forward (resp: backward) if S ðt ÞReðc1 ð0ÞÞ40 ðresp:o0Þ, and the bifurcating periodic solutions on the center manifold are stable (resp: unstable) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. In particular, the periodic solution bifurcated from the first critical value tm is stable (resp: unstable) in phase space and the direction is forward (resp: backward) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. Moreover, when t ¼ tm , ðx ; y Þ is stable (resp: unstable) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. 4. Simulations In this section, we shall carry out some numerical simulations to support our theoretical analysis. With the data given by Kar and Pahari [9], we investigate the impact of harvesting effort on the dynamics, that is ðIÞ r ¼ 2; k ¼ 120; a ¼ 5; b ¼ a ¼ 4; b ¼ c ¼ 1; b0 ¼ 1:5; d0 ¼ 2:5; q ¼ 0:4; g ¼ 0:009: 0 Firstly, let E=3. We have d0 þ qE ¼ 3:7 and Tm ¼ 3:84, hence (A1 ) is satisfied. This implies there exists a unique inner equilibrium E . When t ¼ 0, (H1) is met. Therefore, Proposition 2.3 tells E is globally asymptotically stable. When t40, we restrict t on [0,1000] due to biological meaning. On this interval, (H2) holds and Fig. 1 shows: (i) when t 2 ½0; 33:155Þ, a1 ðtÞ40 and a2 ðtÞr0; (ii) when t 2 ½33:155; 1000, a1 ðtÞ40 and a2 ðtÞZ0. (i) implies (H3) holds and (ii) indicates (H4) is not satisfied. Hence, we obtain I ¼ ½0; 33:155 and for any t 2 I, F ðoðtÞ; tÞ ¼ 0 has only one root oþ ðtÞ40. Fig. 2 describes the curves of S0 ; S1 and S2 versus t on I ¼ ½0; 33:155Þ. One can see that only S0 and S1 have zeros on taxis and there are four critical values of t, denoted by tm 0:1665, t1 9:793; t2 31:875 and tM 33:048, respectively. Theorem 2.5 shows E is asymptotically stable for t 2 ½0; tm Þ [ ðtM ; 1000 and unstable for t 2 ðtm ; tM Þ. Furthermore, system (2.2) undergoes a Hopf bifurcation at the positive equilibrium 1
a2(τ)
a1(τ)
3
1.5
0
0
500 τ
1000
0
0
500 τ
Fig. 1. Curves of a1 ðtÞ and a2 ðtÞ on t 2 ½0; 1000 with E=3.
1000
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1110
Fig. 2. Graphs of functions S0, S1 and S2 for t 2 ½0; t1 Þ with parameters given in (I).
y
18
12
6 10
22 x
34
Fig. 3. Periodic solutions bifurcated from ðx ; y Þ with E=3 and t ¼ 0:1674tm 0:1665.
when t reaches each critical value. In addition, we obtain Reðc1 ð0ÞÞjt¼tm 0:0037o0 and Reðc1 ð0ÞÞjt¼tM 0:0182o0, according to the formula given in Section 3. Therefore, Hopf bifurcation of the system (2.1) at tm (resp. tM ) is forward (resp. backward) and the bifurcating periodic solution is orbitally asymptotically stable (see Figs. 3 and 4). Secondly, we change the value of E from 3 to 5 and the remainder are given by (I). In 0 this case, (A1 ) does not hold. Noting that Tm ¼ 3:84, (A1) implies to91:2201. Therefore, we know from Proposition 2.1 that system (2.2) has a unique positive equilibrium if t 2 ½0; 91:2201Þ. In addition, (H1) and (H2) hold on this interval (see Fig. 5). But Fig. 6 shows a1 ðtÞ and a2 ðtÞ are all positive, which implies both (H3) and (H4) are not fulfilled. Therefore, F ðoÞ ¼ 0 has no positive roots in [0,91.2201) and the unique positive equilibrium E of Eq. (2.1) is asymptotically stable for all t 2 ½0; 91:2201Þ. 0 Lastly, choose E=1. Under such data, (A1 ) holds and F ðoÞ has a positive zero oþ ðtÞ on [0,109.719) and another positive zero o ðtÞ on (102.119,109.719), therefore, I=[0,109.719). Fig. 7 describes the curves of both Snþ ðtÞ and Sn ðtÞ, which displays an interesting phenomenon. In particular, each Si7 ðtÞ ði ¼ 0; . . . ; 7Þ has one root on I, denoted
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
20 x(t)
100 50 0 1.3
y
1111
10
1.45 t
1.6 x 104
1.45 t
x 104
y(t)
20
0
50 x
10 0 1.3
100
1.6
Fig. 4. Periodic solutions bifurcated from ðx ; y Þ with E=3 and t ¼ 32:8otM 33:0479. Left: phase diagram; Right: oscillogram.
0.8 B(τ)+D(τ)
A(τ)+C(τ)
2
1.4
0.8
0
τ
0.4
0
91.2201
0
τ
91.2201
Fig. 5. Curves of AðtÞ þ CðtÞ and BðtÞ þ DðtÞ on t 2 ½0; 91:2201 with E=5.
1.8
a2(τ)
a1(τ)
4
2.5
1
0
τ
91.2201
0.9
0
0
τ
91.2201
Fig. 6. Curves of a1 ðtÞ and a2 ðtÞ on t 2 ½0; 91:2201 with E=5. þ þ þ þ þ by t7 i , respectively, both S8 ðtÞ and S9 ðtÞ have two zeros on I, denoted by t8;1 , t8;2 , t9;1 and þ t9;2 , respectively. Careful observation gives that þ þ þ þ þ tm ¼ t þ 0 o t7 ot8;1 ot9;1 ot0 ot1 ot9;2 ot2 o ot6 ot8;2 ot7 ¼ tM
with tm 6:935 and tM 108:928. However, Að0Þ þ Cð0Þo0 and Bð0Þ þ Dð0Þ40 hold. Therefore, there exists only one positive root of the characteristic equation (2.5) in ½0; tm Þ [ ðtM ; 1Þ. The number of positive roots increases by two whenever t passes the
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1112
120
s−0(τ)
80
s+0(τ)
Sn(τ)
40 0 −40 −80
s+9(τ) 0
35.0323
τ
s−7(τ)
102.12 109.719
Fig. 7. Curves of Sn on t 2 ½0; 109:719 with E=1.
critical value from left to right on ½0; t 0 Þ, and then decreases by two in ½t0 ; t7 . Accordingly, ðx ; y Þ is stable for any tZ0. To sum up, we have the following results.
(a) When E ¼ 1, ðx ; y Þ is always unstable for t 2 ½0; 1Þ. (b) When E ¼ 3, there exists stability switches of stability–instability–stability and Hopf bifurcation, and (i) E is locally asymptotically stable when t 2 ½0; 0:1665 [ ½33:0479; 1000 and unstable when t 2 ð0:1665; 33:0479Þ; (ii) if t ¼ 0:1665 or 33.0479, then the cyclic behavior happens near the positive equilibrium. (c) When E ¼ 5, if t 2 ½0; 91:2201Þ, then E is locally asymptotically stable. Kar and Pahari [9] have carried out simulations and found that if the harvesting effort E increases from 3 to 5, then the oscillatory dynamics (E=3) becomes stable (E=5). Obviously, this phenomenon can be explained by our analytical results. Combining the above simulations and theoretical analysis, it appears that when E increases, the stable interval of E where t varies could extend to a larger one and the phenomenon of existence of periodic solution may disappear. In other words, by adjusting the harvesting effort, cyclic behavior can be prevented and the system may approach a stable steady state. Therefore, harvesting effort does play an important role in the dynamical properties of Eq. (2.1).
Acknowledgments The authors are grateful to Professor Sui Sun Cheng and the anonymous reviewers for their helpful comments and valuable suggestions.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1113
References [1] M. Adimy, F. Crauste, S. Ruan, Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases, Bull. Math. Biol. 68 (2006) 2321–2351. [2] W.G. Aiello, H.I. Freedman, A time-delay model of single-species growth with stage structure, Math. Biosci. 101 (1990) 139–153. [3] W.G. Aiello, H.I. Freedman, J. Wu, Analysis of a model representing stage-structured population growth with state-dependent time delay, SIAM J. Appl. Math. 52 (1992) 855–869. [4] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay-dependent parameters, SIAM J. Math. Anal. 33 (2002) 1144–1165. [5] F. Brauer, A. Soudack, Stability regions in predator–prey systems with constant rate prey harvesting, J. Math. Biol. 8 (1979) 55–71. [6] C. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, second ed., Wiley, New York, 1990. [7] J.K. Hale, M.V. Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993. [8] B. Hassard, N. Kazarinoff, Y. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [9] T. Kar, U. Pahari, Modelling and analysis of a prey–predator system with stage-structure and harvesting, Nonlinear Anal. Real World Appl. 8 (2007) 601–609. [10] S. Liu, L. Chen, G. Luo, Y. Jiang, Asymptotic behaviors of competitive Lotka–Volterra system with stage structure, J. Math. Anal. Appl. 271 (2002) 124–138. [11] Y. Qu, J. Wei, Bifurcation analysis in a time-delay model for prey–predator growth with stage-structure, Nonlinear Dyn. 49 (2007) 285–294. [12] S. Ruan, On nonlinear dynamics of predator–prey models with discrete delay, Math. Model. Nat. Phenom. 4 (2009) 140–187. [13] X. Song, L. Chen, Optimal harvesting and stability for a two species competitive system with stage structure, Math. Biosci. 170 (2001) 173–186. [14] Y. Saito, Y. Takeuchi, A time-delay model for prey–predator growth with stage structure, Can. Appl. Math. Q. 11 (2003) 293–302. [15] J. Wei, M.Y. Li, Hopf bifurcation analysis in a delayed Nicholson blowflies equation, Nonlinear Anal. 60 (2005) 1351–1367. [16] J. Wei, S. Ruan, Stability and bifurcation in a neural network model with two delays, Physica D 130 (1999) 255–272.
Bifurcation analysis in a predator–prey system with stage-structure and harvesting$ Ying Qua,, Junjie Weia,b a Department of Mathematics, Harbin Institute of Technology, Harbin 150001, China Department of Mathematics, Harbin Institute of Technology (Weihai), Weihai, Shandong 264209, China
b
Received 18 November 2009; received in revised form 28 February 2010; accepted 29 March 2010
Abstract In this paper, we consider a predator–prey model with stage-structure and harvesting. This model is the same as the one developed by Kar and Pahari (2007) [9], but we make bifurcation analysis more general than their work. In particular, using the approach of Beretta and Kuang (2002) [4], we show that the positive steady state can be destabilized through a Hopf bifurcation. We also investigate the stability and direction of periodic solutions bifurcating from Hopf bifurcation by using the normal form theory and the center manifold theorem presented in Hassard et al. (1981) [8]. Numerical simulations are then carried out as supporting evidences of our analytical results. & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. Keywords: Predator–prey system; Harvesting; Time-delay; Stability; Hopf bifurcation
1. Introduction The dynamic of predator–prey model has long been and will continue to be one of the dominant themes in both ecology and mathematical ecology due to the fact that predator– prey interaction is the fundamental structure in population dynamics. In most of the related studies, stage-structure is modeled by the introduction of a time delay, which leads to retarded functional differential equations. In general, equations of $
This research is supported by the NNSF of China (No. 10771045), the Program of Excellent Team in HIT and the Natural Science Foundation of Heilongjiang Province (A200806). Corresponding author. E-mail addresses: [email protected], [email protected] (Y. Qu). 0016-0032/$32.00 & 2010 The Franklin Institute. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.jfranklin.2010.03.017
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1098
this kind often exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the involved populations to fluctuate (see [2,3,9–14] and the references therein). Harvesting has a strong impact on the dynamic evolution of a population. To a certain extent, it can control the long-term stationary density of population efficiently. However, it can also lead to the incorporation of a positive extinction probability and, therefore, to potential extinction in finite time. The problems of predator–prey systems with harvesting have been discussed by many authors (see [5,13] for example). Recently, Kar and Pahari [9] proposed a predator–prey-type fishery model with Beddington–DeAngelis functional response and selective harvesting of predator species: x axy ; x_ ¼ rx 1 k a þ bx þ cy y_ ¼ b0 e
gt
bx yðttÞ d0 yqEy; a þ bx þ cy
ð1:1Þ
where x(t) and y(t) denote the population densities of prey and predator at time t, respectively. All coefficients are positive constants. Specially, r and k represent the intrinsic growth rate and carrying capacity of prey in the absence of predation, respectively. The predator consumes the prey with functional response axy=ða þ bx þ cyÞ, known as Beddington–DeAngelis response and contributes to its growth rate bxy=ða þ bx þ cyÞ. b0 is a constant birth rate and d0 is the maximum death rate in the absence of food. All the metabolic energy a predator obtains from its food goes into physiological maintenance, rather than into enhancing reproductive effort. In addition, they assumed that juveniles suffer a mortality rate of g (the through-stage death rate) and take t units of time to mature. Therefore, the term b0 egt yðttÞ represents the population of juveniles who were born at time tt and survive at time t (with the immature death rate g), and thus represents the transformation of juveniles to adults. Moreover, the catch rate function qEy is based on the catch-per-unit-effort (CPUE) hypothesis in [6]. Here, q is the catchability coefficient of the predator species and E is the harvesting effort. In [9], when t ¼ 0, the authors established the global asymptotical stability for the interior equilibrium provided the local asymptotical stability holds. When t40, the authors considered the special case where juveniles mortality rate is g ¼ 0. Under such a hypothesis, the dependence of coefficients on t in the characteristic equation associated with Eq. (1.1) vanishes, because the critical item egt leading to the dependence becomes the constant 1. For the case ga0, they stated that the analytical computations were difficult. Therefore, they only gave the numerical simulation to exhibit the phenomena that stability switches occurred when t is varied. As is well known, for some species (e.g., Antarctic whale and seal populations), delay (the length of time to maturity) is a function of the amount of food (mostly krill), therefore, taking t as parameter is biologically meaningful. From this point of view and as a significant extension of [9], we will give a detailed analysis for the case ga0 which Kar and Pahari failed to examine. Ultimately, the investigations will show the sensitivity of the model dynamics on maturation time delay. We remark that, besides the time delay parameter t, the other parameters in our model are also important biological factors. However, when discussing the dynamical behavior of our system, the key step is to analyze the characteristic equation (2.5) (see below), which is
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1099
a transcendental equation with coefficients depending on the parameters q, g, E,y, in a manner that is quite complicated. Therefore in the present study, we choose to add some simulation results to illustrate the influence of the harvesting effort E on the dynamical behavior of our system (see the last section) but leave out the corresponding analyses to future investigations. This paper is organized as follows. Section 2 gives the analysis of existence and uniqueness of positive equilibrium as well as its stability, by studying the characteristic equation associated with Eq. (1.1) which takes the form of a second degree exponential polynomial with delay-dependent coefficients. Section 3 is devoted to the investigation of the stability and direction of periodic solutions bifurcating from Hopf bifurcation by using the normal form theory and the center manifold theorem presented in Hassard et al. [8]. Finally, some numerical simulations will be given to illustrate our theoretical results in Section 4.
2. Interior equilibrium and its stability We rewrite Eq. (1.1) here for the sake of convenience: x axy ; x_ ¼ rx 1 k a þ bx þ cy y_ ¼ b0 egt yðttÞ d0
bx yqEy; a þ bx þ cy
ð2:1Þ
for tZ0 with the initial values yðyÞ ¼ fðyÞZ0; y 2 ½t;0 and xð0Þ40; yð0Þ40:
2.1. Existence and uniqueness of interior equilibrium We consider the biologically meaningful interior equilibrium E , which satisfies x ay r 1 ; ¼ k a þ bx þ cy TðtÞ ¼
bx ; a þ bx þ cy
where TðtÞ ¼ d0 þ qEb0 egt : Note from Eq. (2.2) that cbr 2 cbr x þ bTðtÞb xTðtÞa ¼ 0; ka a y¼
br ðkxx2 Þ: kaTðtÞ
ð2:2Þ
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1100
For convenience, define bk b cr ; þ1 Tm ¼ min kb þ a b a and denote x7 ¼
pffiffiffiffi k ½ðaTðtÞb þ cbrabÞ7a D; 2cbr
y7 ¼
br ½kx7 ðx7 Þ2 kaTðtÞ
ð2:3Þ
with
cbr D9 TðtÞb þ b a
2 þ 4TðtÞa
cbr : ka
Proposition 2.1. Assume (A1) 0oTðtÞoTm . Then system (2.1) has a unique positive equilibrium E ¼ ðxþ ; yþ Þ. Proof. TðtÞ40 implies that xþ 404x . Accordingly, in order to verify yþ 40, we only need to examine xþ ok. By the expression of xþ given by Eq. (2.3), it follows that pffiffiffiffi cbr TðtÞb þ b: ð2:4Þ Do a This shows that we need cbr=aTðtÞb þ b40; which is equivalent to b cr TðtÞo þ1 : b a Hence, the inequality (2.4) is equivalent to TðtÞo
bk : kb þ a
This completes the proof.
&
Note that TðtÞ is a monotone increasing function of t such that d0 þ qEb0 ¼ Tð0ÞrTðtÞoTð1Þ ¼ d0 þ qE: Therefore, (A1) holds if 0
(A1 ) b0 od0 þ qErTm is satisfied. In what follows, we always assume (A1) holds and we also denote ðxþ ; yþ Þ by ðx ; y Þ.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1101
2.2. Stability analysis Linearizing Eq. (2.1) around ðx ; y Þ gives
2x ay ða þ cy Þ ax ða þ bx Þ _ ¼ r 1 x y; xðtÞ k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 by ða þ cy Þ bx ða þ bx Þ _yðtÞ ¼ b0 egt yðttÞ þ xþ d0 qE y; ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 and its characteristic equation is l2 þ AðtÞl þ BðtÞ þ ½CðtÞl þ DðtÞelt ¼ 0;
ð2:5Þ
where
2x ay ða þ cy Þbx ða þ bx Þ AðtÞ ¼ d0 þ qEr 1 þ ; k ða þ bx þ cy Þ2
ay ða þ cy Þ 2x bx ða þ bx Þ d ; r 1 þ qE BðtÞ ¼ ðd0 þ qEÞ 0 k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 CðtÞ ¼ b0 egt o0;
ay ða þ cy Þ 2x DðtÞ ¼ b0 egt : þ r 1 k ða þ bx þ cy Þ2
Note that these coefficients are all dependent on t. In addition, AðtÞ; BðtÞ; CðtÞ and DðtÞ are all continuous and differentiable functions from ½0; 1Þ to R. We now investigate the local asymptotic stability of ðx ; y Þ for system (2.1). It is well known that the trivial solution of system (2.1) is locally asymptotically stable if all roots l of the characteristic equation (2.5) satisfy ReðlÞo0: Eq. (2.5) takes the general form Pðl;tÞ þ Qðl;tÞelt ¼ 0; where
ð2:6Þ
Pðl;tÞ9l2 þ AðtÞl þ BðtÞ; Qðl;tÞ9CðtÞl þ DðtÞ:
ð2:7Þ
Let us first examine the case t ¼ 0. This case is of importance, because if the nontrivial positive steady-state of Eq. (2.1) is stable when t ¼ 0, then we are able to obtain the local stability for all nonnegative values of the delay, or to find a critical value which could destabilize the steady-state. When t ¼ 0, the characteristic equation (2.5) reduces to l2 þ ½Að0Þ þ Cð0Þl þ ½Bð0Þ þ Dð0Þ ¼ 0: Clearly, all roots of Eq. (2.6) with t ¼ 0 have negative real parts provided (H1) Að0Þ þ Cð0Þ40 and Bð0Þ þ Dð0Þ40: Hence, we get the following conclusion.
1102
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
Proposition 2.2 (Kar and Pahari [9, Theorem 2.1]). Assume (A1) and (H1) hold. Then the positive equilibrium E ¼ ðx ; y Þ of Eq. (2.1) is globally asymptotically stable if t ¼ 0. Now we assume t40 and regard it as a bifurcation parameter to obtain finer results on the stability of ðx ; y Þ. Note that Eq. (2.5) takes the form of a second degree exponential polynomial in l, with all the coefficients of P and Q depending on t. Thus, we use the method introduced by Beretta and Kuang [4], which gives the existence of purely imaginary roots of a characteristic equation with delay-dependent coefficients (see also [1,11]). In order to apply the criterion in [4], we need to verify the following properties for o40 and t 2 I with I defined in Eq. (2.8): (i) (ii) (iii) (iv) (v)
Pð0; tÞ þ Qð0; tÞa0; Pðio; tÞ þ Qðio; tÞa0; lim supjlj-1;RelZ0 jQðl; tÞ=Pðl; tÞjo1; F ðoÞ9jPðio; tÞj2 jQðio; tÞj2 has a finite number of zeros; each positive root oðtÞ of F ðoðtÞÞ ¼ 0 is continuous and differentiable in t whenever it exists.
Here Pðl; tÞ and Qðl; tÞ are defined by Eq. (2.7) and I9ftZ0 : there exists oðtÞ40 such that F ðoðtÞÞ ¼ 0g:
ð2:8Þ
Let F be defined as in (iv). By Eq. (2.7), we obtain F ðoÞ ¼ o4 þ a1 ðtÞo2 þ a2 ðtÞ
ð2:9Þ
with a1 ðtÞ ¼ A2 ðtÞ2BðtÞC 2 ðtÞ; a2 ðtÞ ¼ B2 ðtÞD2 ðtÞ: Before proceeding further, let us analyze the structure of I described by Eq. (2.8). Set qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a1 ðtÞ7 a21 ðtÞ4a2 ðtÞ : z7 ðtÞ ¼ 2 We need the following hypotheses: (H2) AðtÞ þ CðtÞa0 and BðtÞ þ DðtÞa0: (H3) 0oa2 ðtÞoa21 ðtÞ=4 and a1 ðtÞo0: (H4) a2 ðtÞr0: Proposition 2.3. Assume (A1) holds. pffiffiffiffiffiffiffiffiffiffiffi (i) If (H3) is satisfied, then F ðoÞ ¼ 0 has two different positive roots z7 ðtÞ denoted by o7 , respectively. pffiffiffiffiffiffiffiffiffiffiffi (ii) If (H4) is met, then F ðoÞ ¼ 0 has a unique positive root oþ ¼ zþ ðtÞ.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1103
Define I1 ¼ ftZ0 : ðA1Þ; ðH2Þ and ðH3Þ holdg; I2 ¼ ftZ0 : ðA1Þ; ðH2Þ and ðH4Þ holdg: Then I ¼ I1 [ I2 :
ð2:10Þ
Now, we verify the properties (i)–(v) for t 2 I. Firstly, (i) and (ii) are satisfied. Indeed, Pð0;tÞ þ Qð0;tÞ ¼ BðtÞ þ DðtÞa0 and Pðio;tÞ þ Qðio;tÞ ¼ ½BðtÞ þ DðtÞo2 þ io½AðtÞ þ CðtÞa0: Then, from Eq. (2.7) we know Qðl;tÞ ¼ 0: lim jlj-1 Pðl;tÞ Therefore (iii) follows. Finally, Eq. (2.9) implies (iv) and (v). Let l ¼ io (o40) be a root of Eq. (2.5). Substituting it into Eq. (2.5) and separating the real and imaginary parts yields o2 BðtÞ ¼ DðtÞcosot þ CðtÞosinot; AðtÞo ¼ DðtÞsinotCðtÞocosot; then we have sinot ¼
ðo2 BÞCo þ oAD o2 C 2 þ D2
and
cosot ¼
ðBo2 ÞD þ o2 AC : o2 C 2 þ D2
ð2:11Þ
Here the dependence of the coefficients on t is implicitly assumed. By the definitions of Pðl; tÞ and Qðl; tÞ in Eq. (2.7) and applying the property (ii), Eq. (2.11) can be written as
Pðio;tÞ Pðio;tÞ and cosot ¼ Re ; sinot ¼ Im Qðio;tÞ Qðio;tÞ which lead us to F ðoÞ ¼ 0:
ð2:12Þ
Noting that Proposition 2.3 has given the explicit expressions of oðtÞ that satisfies F ðoðtÞÞ ¼ 0 on t 2 I, we can define yðtÞ 2 ½0; 2p for t 2 I: sinyðtÞ ¼
ðo2 BÞCo þ oAD o2 C 2 þ D2
and
cosyðtÞ ¼
ðBo2 ÞD þ o2 AC : o2 C 2 þ D2
Such yðtÞ defined above is well and uniquely defined for all t 2 I (see [4]).
ð2:13Þ
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1104
Therefore, io with o ¼ oðt Þ40 is a purely imaginary root of Eq. (2.5) if and only if t is a zero of the Sn, where
Sn ðtÞ9t
yðtÞ þ 2np ; t2I oðtÞ
with n 2 N0 9f0; 1; 2; . . .g. The following result is due to Beretta and Kuang [4]. Lemma 2.4. The characteristic equation (2.5) has a pair of simple pure imaginary roots l ¼ 7ioðt Þ at t 2 I, provided Sn ðt Þ ¼ 0 for n 2 N0 . Moreover, if oðt Þ ¼ oþ ðt Þ, then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if dþ ðt Þ40, and crosses the imaginary axis from right to left if dþ ðt Þo0, where ) dReðlÞ dSn ðtÞ ; dþ ðt Þ ¼ Sign ¼ Sign dt l¼ioþ ðt Þ dt t¼t and if oðt Þ ¼ o ðt Þ, then this pair of simple conjugate purely imaginary roots crosses the imaginary axis from left to right if d ðt Þ40, and crosses the imaginary axis from right to left if d ðt Þo0, where ) dReðlÞ dSn ðtÞ : d ðt Þ ¼ Sign ¼ Sign dt l¼io ðt Þ dt t¼t If I ¼ I2 , then only oþ is feasible. In this case, we can easily observe that Sn ð0Þr0 and Sn ðtÞ4Snþ1 ðtÞ for t 2 I and n 2 N0 . Therefore, Sn has no zeros in I for any n 2 N0 if S0 has no zeros in I. As a result, Eq. (2.1) has no scenario at ðx ; y Þ. If Sn ðtÞ has positive zeros, denoted by tjn , for some n 2 N0 , then without loss of generality, we may suppose dSn ðtjn Þ a0 dt
with Sn ðtjn Þ ¼ 0:
ð2:14Þ
Then for one thing, stability switches occur at the zeros of S0 ðtÞ, denoted by tj0 , if (H1) holds. For another, applying Proposition 2.2 and Hopf bifurcation theorem for functional differential equations (see Hale’s book [7]), we can conclude the existence of a Hopf bifurcation. Before stating the main theorem, denote tm ¼ minft 2 I : S0 ðtÞ ¼ 0g
and
tM ¼ maxft 2 I : S0 ðtÞ ¼ 0g:
Theorem 2.5. Assume (A1), (H1), (H2) and (H4) hold. (i) If S0 ðtÞ has no positive zeros in I, then ðx ; y Þ is locally asymptotically stable for all tZ0. (ii) If Sn ðtÞ has at least one positive zero in I and Eq. (2.14) is met, then ðx ; y Þ is locally asymptotically stable for t 2 ½0; tm Þ [ ðtM ; 1Þ and unstable for t 2 ðtm ; tM Þ. That is, stability switches of stability–instability–stability occur. Hopf bifurcation takes place when t ¼ tjn , n 2 N 0 . If I ¼ I1 , then both oþ and o are feasible for t 2 I. We have the following two sequences of functions on I: Snþ ðtÞ ¼ tðyþ ðtÞ þ 2npÞ=oþ ðtÞ
and
Sn ðtÞ ¼ tðy ðtÞ þ 2npÞ=o ðtÞ;
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1105
where y7 ðtÞ is the solution of Eq. (2.13) when o ¼ o7 . Similarly, it is obtained for t 2 I 7 ðtÞ with n 2 N0 . Furthermore, if S0þ ðtÞ4S0 ðtÞ, then that Sn7 ð0Þr0 and Sn7 ðtÞ4Snþ1 þ Sn ðtÞ4Sn ðtÞ for n 2 N0 . Other than the case above, stability switch may depend on all real roots of Snþ ðtÞ ¼ 0 and Sn ðtÞ ¼ 0. In addition, Hopf bifurcations can also appear at each root of Sn7 ðtÞ ¼ 0. Naturally, we can also obtain results similar to Theorem 2.5. Here, we omit the corresponding statements. If I ¼ I1 [ I2 with I1 a| and I2 a|, then we can discuss different cases on their own sets, respectively. The remaining discussion is the same as above.
3. The properties of Hopf bifurcation In the previous section, we obtained some sufficient conditions ensuring system (2.1) undergoes a Hopf bifurcation at E ¼ ðx ; y Þ. In this section, under these conditions such as Theorem 2.5(iii), we shall establish the explicit formula for determining the direction and stability of periodic solutions bifurcating from E at each zero t of Sn , by using normal form theory and center manifold argument presented in Hassard et al. [8] (see also [7,15,16]). By the change of variables xðtÞx /xðtÞ and yðtÞy /yðtÞ and by rescaling the time t/ðt=tÞ, system (2.1) can be written as ! ! ! _ xðtÞ xðtÞ xðt1Þ ¼ tMðtÞ þ tNðtÞ þ tF ðx;yÞ; ð3:1Þ _ yðtÞ yðtÞ yðt1Þ where
0 2x ay ða þ cy Þ r 1 B k ða þ bx þ cy Þ2 B MðtÞ ¼ B @ by ða þ cy Þ ða þ bx þ cy Þ2
NðtÞ ¼
F ðx; yÞ9
!
0
0
0
b0 egt
f ðx; yÞ lðx; yÞ
1 ax ða þ bx Þ C ða þ bx þ cy Þ2 C C; bx ða þ bx Þ A d qE 0 ða þ bx þ cy Þ2
!
;
0
1 rx2 axy B k a þ bx þ cy C B C ¼B C bxy @ A a þ bx þ cy 0 1 1 1 1 2 2 3 x þ f xy þ y þ x þ f f f 11 02 30 B 2 20 C 2 6 C ¼B @1 A 1 1 2 3 2 l20 x þ l11 xy þ l02 y þ l30 x þ 2 2 6
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1106
with fij ¼
@iþj f ðx ;y Þ @i x@j y
and
lij ¼
@iþj lðx ;y Þ ; i;j ¼ 1;2; . . . : @i x@j y
Let t ¼ t þ m. Thus m ¼ 0 is a Hopf bifurcation value of Eq. (3.1). For f ¼ colðf1 ; f2 Þ 2 Cð½1; 0; R2 Þ, define Lm ðfÞ ¼ ðt þ mÞ½Mðt þ mÞfð0Þ þ Nðt þ mÞfð1Þ: By the Riesz representation theorem, there exists a 2 2 matrix Zðy; mÞ; y 2 ½1; 0, whose elements are functions of bounded variation such that Z 0 Lm ðfÞ ¼ ½dZðy;mÞfðyÞ; f 2 Cð½1;0;R2 Þ: ð3:2Þ 1
In fact, we can choose Zðy; mÞ ¼ ðt þ mÞðMðt þ mÞdðyÞCðt þ mÞdðy þ 1ÞÞ; that is
8 > < ðt þ mÞBðt þ mÞ; Zðy; mÞ ¼ 0; > : ðt þ mÞCðt þ mÞ;
y ¼ 0; y 2 ð1; 0Þ; y ¼ 1:
Then Eq. (3.2) is satisfied. For f 2 C 1 ð½1; 0; R2 Þ, define the operator AðmÞ as 8 > < dfðyÞ; y 2 ½1; 0Þ; AðmÞfðyÞ ¼ R dy > : 0 ½dZðx; mÞfðxÞ; y ¼ 0 1 and
( RðmÞfðyÞ ¼
with
0; hðm;
y 2 ½1;0Þ; fÞ;
y¼0
0
1 1 1 1 2 2 3 B 2 f20 f1 ð0Þ þ f11 f1 ð0Þf2 ð0Þ þ 2 f02 f2 ð0Þ þ 6 f30 f1 ð0Þ þ C C: hðm; fÞ ¼ ðt þ mÞB @1 A 1 1 2 2 3 l20 f1 ð0Þ þ l11 f1 ð0Þf2 ð0Þ þ l02 f2 ð0Þ þ l30 f1 ð0Þ þ 2 2 6
Then system (2.1) is equivalent to the following operator equation: u_t ¼ AðmÞut þ RðmÞut ; where uðtÞ ¼ colðxðtÞ; yðtÞÞ and ut ¼ uðt þ yÞ for y 2 ½1; 0: For c 2 C 1 ð½0; 1; ðR2 Þ Þ, define 8 > < dcðsÞ; s 2 ð0; 1; A cðsÞ ¼ R ds > : 0 cðxÞdZðx;0Þ; s ¼ 0 1
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1107
and a bilinear form /cðsÞ;fðyÞS ¼ cð0Þfð0Þ
Z
0
Z
y
cðxyÞ dZðyÞfðxÞ dx; 1
x¼0
where ZðyÞ ¼ Zðy; 0Þ. Then Að0Þ and A are adjoint operators. From the discussion in Section 2, we know that 7io t are eigenvalues of Að0Þ and therefore they are also eigenvalues of A . Moreover, it is not difficult to verify that vectors qðyÞ ¼ colðq1 ; q2 Þeio t y ðy 2 ½1; 0Þ and q ðsÞ ¼ ð1=K Þðq1 ; q2 Þeio t s ðs 2 ½0; 1Þ are the eigenvectors of Að0Þ and A corresponding to the eigenvalues io t and io t , respectively, where q1 ¼ q1 ¼ 1; q2 ¼
rða þ bx þ cy Þ2 ðk2x Þaky ða þ cy Þ o ða þ bx þ cy Þ2 i; ax ða þ bx Þ kax ða þ bx Þ
q2 ¼
aky ða þ cy Þ þ rða þ bx þ cy Þ2 ð2x kÞ o ða þ bx þ cy Þ2 i: by ða þ cy Þ kby ða þ cy Þ
Let K ¼ 1 þ q2 q2 ð1 þ b0 t et
ðio þgÞ
Þ:
Then /q ðsÞ; qðyÞS ¼ 1 and /q ðsÞ; qðyÞS ¼ 0. Following the same algorithms and using a similar computation process as in Hassard et al. [8], we obtain the coefficients determining the important quantities: g20 ¼ g11 ¼
t ½q ðf20 q21 þ 2f11 q1 q2 þ f02 q22 Þ þ q 2 ðl20 q21 þ 2l11 q1 q2 þ l02 q22 Þ; K 1 t fq ðf20 jq1 j2 þ f11 ðq1 q 2 þ q 1 q2 Þ þ f02 jq2 j2 Þ K 1 þ q 2 ðl20 jq1 j2 þ l11 ðq1 q 2 þ q 1 q2 Þ þ l02 jq2 j2 Þg;
g02 ¼
t ½q ðf20 q 21 þ 2f11 q 1 q 2 þ f02 q 22 Þ þ q 2 ðl20 q 21 þ 2l11 q 1 q 2 þ l02 q 22 Þ; K 1
g21 ¼
t ð1Þ ð1Þ fq ½f20 ðW20 ð0Þq 1 þ 2W11 ð0Þq1 Þ K 1 ð1Þ ð1Þ ð2Þ ð2Þ þ f11 ðW20 ð0Þq 2 þ 2W11 ð0Þq2 þ 2W11 ð0Þq1 þ W20 ð0Þq 1 Þ ð2Þ ð2Þ ð0Þq 2 þ 2W11 ð0Þq2 Þ þ f30 q21 q 1 þ f21 ðq21 q 2 þ 2jq1 j2 q2 Þ þ f02 ðW20
þ f12 ðq22 q 1 þ 2jq2 j2 q1 Þ þ f03 q22 q 2 ð1Þ ð1Þ ð0Þq 1 þ 2W11 ð0Þq1 Þ þ q 2 ½l20 ðW20 ð1Þ ð1Þ ð2Þ ð2Þ ð0Þq 2 þ 2W11 ð0Þq2 þ 2W11 ð0Þq1 þ W20 ð0Þq 1 Þ þ l11 ðW20 ð2Þ ð2Þ þ l02 ðW20 ð0Þq 2 þ 2W11 ð0Þq2 Þ þ l30 q21 q 1 þ l21 ðq21 q 2 þ 2jq1 j2 q2 Þ
þ l12 ðq22 q 1 þ 2jq2 j2 q1 Þ þ l03 q22 q 2 g;
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1108
where W20 ðyÞ ¼
ig20 ig qðyÞ þ 02 qðyÞ þ E1 e2io t y ; o t 3o t
W11 ðyÞ ¼
ig11 ig qðyÞ þ 11 qðyÞ þ E2 o t o t
and
0 2x ay ða þ cy Þ r 1 2io B k ða þ bx þ cy Þ2 B E1 ¼ 2B by ða þ cy Þ @ ða þ bx þ cy Þ2 0 1 1 1 2 2 f f q þ f q q þ q 11 1 2 02 2 C B 2 20 1 2 C; B @1 A 1 l20 q21 þ l11 q1 q2 þ l02 q22 2 2
11 ax ða þ bx Þ ða þ bx þ cy Þ2 C C C A M
0 11 2x ay ða þ cy Þ ax ða þ bx Þ r 1 B C k ða þ bx þ cy Þ2 ða þ bx þ cy Þ2 B C E 2 ¼ B C by ða þ cy Þ bx ða þ bx Þ @ A d0 qE þ b0 egt 2 2 ða þ bx þ cy Þ ða þ bx þ cy Þ 0 1 1 1 2 2 f f jq j þ f ðq q þ q q Þ þ jq j Þ 20 1 11 1 2 02 2 2 1 B2 C 2 C B @1 A 1 2 2 l20 jq1 j þ l11 ðq1 q 2 þ q 1 q2 Þ þ l02 jq2 j Þ 2 2 with M¼
bx ða þ bx Þ d0 qE þ b0 egt2io t 2io : 2 ða þ bx þ cy Þ
Consequently, we can compute the following values: 2
i g21 2 jg02 j c1 ð0Þ ¼ þ ; g11 g20 2jg11 j 3 2 2o t m2 ¼
Reðc1 ð0ÞÞ ; 0 Reðl ðt ÞÞ
b2 ¼ 2 Reðc1 ð0ÞÞ; 0
T2 ¼
Imðc1 ð0ÞÞ þ m2 Imðl ðt ÞÞ ; o t
which determine the properties of bifurcating periodic solutions at the critical value t . To be concrete, m2 determines the direction of the Hopf bifurcation: if m2 40 (resp. o0), then Hopf bifurcations are forward (resp. backward); b2 determines the stability of bifurcating
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1109
periodic solutions: the bifurcating periodic solutions on the center manifold are stable (resp. unstable) if b2 o0 (resp. 40); T2 determines the periods of the bifurcating periodic solutions: the periods increase (resp. decrease) if T2 40 (resp. o0). Combining these properties with discussions in Section 2, we have: Theorem 3.1. Under the conditions given in Theorem 2.5(iii), Hopf bifurcation at each t is 0 forward (resp: backward) if S ðt ÞReðc1 ð0ÞÞ40 ðresp:o0Þ, and the bifurcating periodic solutions on the center manifold are stable (resp: unstable) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. In particular, the periodic solution bifurcated from the first critical value tm is stable (resp: unstable) in phase space and the direction is forward (resp: backward) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. Moreover, when t ¼ tm , ðx ; y Þ is stable (resp: unstable) if Reðc1 ð0ÞÞo0 ðresp: 40Þ. 4. Simulations In this section, we shall carry out some numerical simulations to support our theoretical analysis. With the data given by Kar and Pahari [9], we investigate the impact of harvesting effort on the dynamics, that is ðIÞ r ¼ 2; k ¼ 120; a ¼ 5; b ¼ a ¼ 4; b ¼ c ¼ 1; b0 ¼ 1:5; d0 ¼ 2:5; q ¼ 0:4; g ¼ 0:009: 0 Firstly, let E=3. We have d0 þ qE ¼ 3:7 and Tm ¼ 3:84, hence (A1 ) is satisfied. This implies there exists a unique inner equilibrium E . When t ¼ 0, (H1) is met. Therefore, Proposition 2.3 tells E is globally asymptotically stable. When t40, we restrict t on [0,1000] due to biological meaning. On this interval, (H2) holds and Fig. 1 shows: (i) when t 2 ½0; 33:155Þ, a1 ðtÞ40 and a2 ðtÞr0; (ii) when t 2 ½33:155; 1000, a1 ðtÞ40 and a2 ðtÞZ0. (i) implies (H3) holds and (ii) indicates (H4) is not satisfied. Hence, we obtain I ¼ ½0; 33:155 and for any t 2 I, F ðoðtÞ; tÞ ¼ 0 has only one root oþ ðtÞ40. Fig. 2 describes the curves of S0 ; S1 and S2 versus t on I ¼ ½0; 33:155Þ. One can see that only S0 and S1 have zeros on taxis and there are four critical values of t, denoted by tm 0:1665, t1 9:793; t2 31:875 and tM 33:048, respectively. Theorem 2.5 shows E is asymptotically stable for t 2 ½0; tm Þ [ ðtM ; 1000 and unstable for t 2 ðtm ; tM Þ. Furthermore, system (2.2) undergoes a Hopf bifurcation at the positive equilibrium 1
a2(τ)
a1(τ)
3
1.5
0
0
500 τ
1000
0
0
500 τ
Fig. 1. Curves of a1 ðtÞ and a2 ðtÞ on t 2 ½0; 1000 with E=3.
1000
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1110
Fig. 2. Graphs of functions S0, S1 and S2 for t 2 ½0; t1 Þ with parameters given in (I).
y
18
12
6 10
22 x
34
Fig. 3. Periodic solutions bifurcated from ðx ; y Þ with E=3 and t ¼ 0:1674tm 0:1665.
when t reaches each critical value. In addition, we obtain Reðc1 ð0ÞÞjt¼tm 0:0037o0 and Reðc1 ð0ÞÞjt¼tM 0:0182o0, according to the formula given in Section 3. Therefore, Hopf bifurcation of the system (2.1) at tm (resp. tM ) is forward (resp. backward) and the bifurcating periodic solution is orbitally asymptotically stable (see Figs. 3 and 4). Secondly, we change the value of E from 3 to 5 and the remainder are given by (I). In 0 this case, (A1 ) does not hold. Noting that Tm ¼ 3:84, (A1) implies to91:2201. Therefore, we know from Proposition 2.1 that system (2.2) has a unique positive equilibrium if t 2 ½0; 91:2201Þ. In addition, (H1) and (H2) hold on this interval (see Fig. 5). But Fig. 6 shows a1 ðtÞ and a2 ðtÞ are all positive, which implies both (H3) and (H4) are not fulfilled. Therefore, F ðoÞ ¼ 0 has no positive roots in [0,91.2201) and the unique positive equilibrium E of Eq. (2.1) is asymptotically stable for all t 2 ½0; 91:2201Þ. 0 Lastly, choose E=1. Under such data, (A1 ) holds and F ðoÞ has a positive zero oþ ðtÞ on [0,109.719) and another positive zero o ðtÞ on (102.119,109.719), therefore, I=[0,109.719). Fig. 7 describes the curves of both Snþ ðtÞ and Sn ðtÞ, which displays an interesting phenomenon. In particular, each Si7 ðtÞ ði ¼ 0; . . . ; 7Þ has one root on I, denoted
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
20 x(t)
100 50 0 1.3
y
1111
10
1.45 t
1.6 x 104
1.45 t
x 104
y(t)
20
0
50 x
10 0 1.3
100
1.6
Fig. 4. Periodic solutions bifurcated from ðx ; y Þ with E=3 and t ¼ 32:8otM 33:0479. Left: phase diagram; Right: oscillogram.
0.8 B(τ)+D(τ)
A(τ)+C(τ)
2
1.4
0.8
0
τ
0.4
0
91.2201
0
τ
91.2201
Fig. 5. Curves of AðtÞ þ CðtÞ and BðtÞ þ DðtÞ on t 2 ½0; 91:2201 with E=5.
1.8
a2(τ)
a1(τ)
4
2.5
1
0
τ
91.2201
0.9
0
0
τ
91.2201
Fig. 6. Curves of a1 ðtÞ and a2 ðtÞ on t 2 ½0; 91:2201 with E=5. þ þ þ þ þ by t7 i , respectively, both S8 ðtÞ and S9 ðtÞ have two zeros on I, denoted by t8;1 , t8;2 , t9;1 and þ t9;2 , respectively. Careful observation gives that þ þ þ þ þ tm ¼ t þ 0 o t7 ot8;1 ot9;1 ot0 ot1 ot9;2 ot2 o ot6 ot8;2 ot7 ¼ tM
with tm 6:935 and tM 108:928. However, Að0Þ þ Cð0Þo0 and Bð0Þ þ Dð0Þ40 hold. Therefore, there exists only one positive root of the characteristic equation (2.5) in ½0; tm Þ [ ðtM ; 1Þ. The number of positive roots increases by two whenever t passes the
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1112
120
s−0(τ)
80
s+0(τ)
Sn(τ)
40 0 −40 −80
s+9(τ) 0
35.0323
τ
s−7(τ)
102.12 109.719
Fig. 7. Curves of Sn on t 2 ½0; 109:719 with E=1.
critical value from left to right on ½0; t 0 Þ, and then decreases by two in ½t0 ; t7 . Accordingly, ðx ; y Þ is stable for any tZ0. To sum up, we have the following results.
(a) When E ¼ 1, ðx ; y Þ is always unstable for t 2 ½0; 1Þ. (b) When E ¼ 3, there exists stability switches of stability–instability–stability and Hopf bifurcation, and (i) E is locally asymptotically stable when t 2 ½0; 0:1665 [ ½33:0479; 1000 and unstable when t 2 ð0:1665; 33:0479Þ; (ii) if t ¼ 0:1665 or 33.0479, then the cyclic behavior happens near the positive equilibrium. (c) When E ¼ 5, if t 2 ½0; 91:2201Þ, then E is locally asymptotically stable. Kar and Pahari [9] have carried out simulations and found that if the harvesting effort E increases from 3 to 5, then the oscillatory dynamics (E=3) becomes stable (E=5). Obviously, this phenomenon can be explained by our analytical results. Combining the above simulations and theoretical analysis, it appears that when E increases, the stable interval of E where t varies could extend to a larger one and the phenomenon of existence of periodic solution may disappear. In other words, by adjusting the harvesting effort, cyclic behavior can be prevented and the system may approach a stable steady state. Therefore, harvesting effort does play an important role in the dynamical properties of Eq. (2.1).
Acknowledgments The authors are grateful to Professor Sui Sun Cheng and the anonymous reviewers for their helpful comments and valuable suggestions.
Y. Qu, J. Wei / Journal of the Franklin Institute 347 (2010) 1097–1113
1113
References [1] M. Adimy, F. Crauste, S. Ruan, Modelling hematopoiesis mediated by growth factors with applications to periodic hematological diseases, Bull. Math. Biol. 68 (2006) 2321–2351. [2] W.G. Aiello, H.I. Freedman, A time-delay model of single-species growth with stage structure, Math. Biosci. 101 (1990) 139–153. [3] W.G. Aiello, H.I. Freedman, J. Wu, Analysis of a model representing stage-structured population growth with state-dependent time delay, SIAM J. Appl. Math. 52 (1992) 855–869. [4] E. Beretta, Y. Kuang, Geometric stability switch criteria in delay differential systems with delay-dependent parameters, SIAM J. Math. Anal. 33 (2002) 1144–1165. [5] F. Brauer, A. Soudack, Stability regions in predator–prey systems with constant rate prey harvesting, J. Math. Biol. 8 (1979) 55–71. [6] C. Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources, second ed., Wiley, New York, 1990. [7] J.K. Hale, M.V. Lunel, Introduction to Functional Differential Equations, Springer, New York, 1993. [8] B. Hassard, N. Kazarinoff, Y. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, Cambridge, 1981. [9] T. Kar, U. Pahari, Modelling and analysis of a prey–predator system with stage-structure and harvesting, Nonlinear Anal. Real World Appl. 8 (2007) 601–609. [10] S. Liu, L. Chen, G. Luo, Y. Jiang, Asymptotic behaviors of competitive Lotka–Volterra system with stage structure, J. Math. Anal. Appl. 271 (2002) 124–138. [11] Y. Qu, J. Wei, Bifurcation analysis in a time-delay model for prey–predator growth with stage-structure, Nonlinear Dyn. 49 (2007) 285–294. [12] S. Ruan, On nonlinear dynamics of predator–prey models with discrete delay, Math. Model. Nat. Phenom. 4 (2009) 140–187. [13] X. Song, L. Chen, Optimal harvesting and stability for a two species competitive system with stage structure, Math. Biosci. 170 (2001) 173–186. [14] Y. Saito, Y. Takeuchi, A time-delay model for prey–predator growth with stage structure, Can. Appl. Math. Q. 11 (2003) 293–302. [15] J. Wei, M.Y. Li, Hopf bifurcation analysis in a delayed Nicholson blowflies equation, Nonlinear Anal. 60 (2005) 1351–1367. [16] J. Wei, S. Ruan, Stability and bifurcation in a neural network model with two delays, Physica D 130 (1999) 255–272.