Stability, Nonlinear Oscillations and Bifurcation in a ... Nonlinear Oscillations and Bifurcation in a Delay-Induced Predator-Prey System with ... mechanism which takes τunits of time to respond to changes

  • Published on

  • View

  • Download

Embed Size (px)


<ul><li><p>Stability, Nonlinear Oscillations and Bifurcation ina Delay-Induced Predator-Prey System with</p><p>HarvestingDebaldev Jana, Swapan Chakraborty and Nandadulal Bairagi</p><p>AbstractA harvested predator-prey system that incorpo-rates feedback delay in prey growth rate is studied. Consideringdelay as a parameter, we investigate the effect of delay on thestability of the coexisting equilibrium. It is observed that thereexists a critical length of the delay parameter below whichthe coexistence equilibrium is stable and above which it isunstable. A Hopf bifurcation exists when delay crosses thecritical value. By applying the normal form theory and thecenter manifold theorem, we determine the explicit formulaethat demonstrate the stability and direction of the bifurcatingperiodic solutions. Computer simulations have been carriedout to illustrate different analytical findings. Our simulationresults indicate that the Hopf bifurcation is supercritical and thebifurcating periodic solutions are unstable for the consideredparameter values. The system, however, can be stabilized fromits unstable oscillatory state only by lowering the intrinsicgrowth rate of prey population.</p><p>Index Termspredator-prey model, harvesting, delay, stabil-ity, oscillations, direction of Hopf bifurcation.</p><p>I. INTRODUCTION</p><p>THE most studied topics in theoretical ecology is thepredator-prey interaction and its possible outcomes.This is because of its all over existence in the naturalworld. Various mathematical models have been proposedfrom different biological point of views to study the in-teraction between prey and predator [1][3]. Harvesting isan old technique used for biological resource exploitation.However, overexploitation is the most responsible factor tothe persistence of harvested stocks and causes not onlyeconomic loss but also changes the ecosystem structure andfunctioning [4][9]. Recently, Bairagi et al. [10] consideredthe following predator-prey model with harvesting:</p><p>dxdt</p><p>= rx(1 xk) xy</p><p>ay+x ,dydt</p><p>= b0xyay+x d0y </p><p>qEybE+ly ,</p><p>(1)</p><p>where x(t) and y(t) be, respectively, the prey and predatordensities at time t. It is assumed that the prey populationgrows logistically to its carrying capacity k with intrinsicgrowth rate r. b0 and d0 be, respectively, the conversion effi-ciency and natural death rate of predator. Functional responseof predator is assumed to follow ratio-dependent Type II formwith maximum prey consumption rate and half-saturationconstant a. Predator is harvested following the Michaelis-Menten type catch rate with catchability coefficient q and</p><p>Manuscript received August 21, 2012. This work was supported by DRSProgramme (UGC No. F510/2/DRS/2012 (SAP-1)).</p><p>D. Jana and N. Bairagi are with the Centre for Mathematical Biology andEcology, Department of Mathematics, Jadavpur University, Kolkata-700032,INDIA, e-mail:</p><p>S. Chakraborty is with Baruipara High School, Halisahar, 24 Parganas(N), West Bengal, India.</p><p>harvesting effort E. All parameters including b and l areassumed to be positive. Using Melnikovs method, they ob-served that the system (1) exhibits heteroclinic bifurcations.They also showed that the system may exhibit monostability,bistability and tristability depending on the initial values ofthe system populations and the harvesting effort.</p><p>Delay is frequently used in a predator-prey model torepresent the biological process more accurately. Predator-prey models with discrete delay exhibit different interestingdynamics, such as the existence of Hopf bifurcation, andBogdanov-Takens bifurcation [11], [12]. An excellent reviewwork of predator-prey models with single discrete delay canbe seen in [13]. Simultaneous effects of harvesting and delayon predator-prey system were studied by several researchers[14][18]. The general observation is that the time delaymakes a stable equilibrium unstable and harvesting, on theother hand, helps to regain the stability. In this study, weconsider a discrete delay in the specific growth rate of preyto incorporate the effect of density dependence feedbackmechanism which takes units of time to respond to changesin the prey population [19]. Thus the model system (1)reduces to the following delay-induced predator-prey modelwith harvesting:</p><p>dxdt</p><p>= rx(1 x(t)k</p><p>) xyay+x ,</p><p>dydt</p><p>= y[ b0xay+x d0 </p><p>qEbE+ly ].</p><p>(2)</p><p>The initial condition is</p><p>x() = 1() 0, y() = 2() 0, (, 0].</p><p>Delay in the system (2) is based on the assumption that inthe absence of predators the prey satisfies the delayed logisticequation [20]. The objective is to study the dynamic behaviorof the system (2) and determine the direction and stabilityof the bifurcating periodic solutions, if any.</p><p>The organization of the paper is as follows: In SectionII, we perform the local stability of the model system.Section III deals with the direction and stability of bifurcatingperiodic solutions. Rigorous numerical simulations of themodel system are performed in the Section IV. Finally, abrief discussion is presented in the Section 4.</p><p>II. LOCAL STABILITY</p><p>Coexistence of species are given utmost importance in apredator-prey interaction and we are, therefore, only inter-ested in the coexistence equilibrium of the system (2). Thecoexistence equilibrium point of the above system is given byE(x, y), where the equilibrium prey density x is given</p><p>Engineering Letters, 20:3, EL_20_3_05</p><p>(Revised online publication: 30 August 2012)</p><p> ______________________________________________________________________________________ </p></li><li><p>TABLE IPARAMETER VALUES USED IN THE MODEL SIMULATION</p><p>Parameter Default value</p><p>r 1</p><p>k 100</p><p> 0.7</p><p>a 10</p><p>b0 0.45</p><p>d0 0.04</p><p>q 1</p><p>b 0.1</p><p>E 0.01</p><p>l 0.04</p><p>by the positive root of the cubic equation</p><p>p1x3 + p2x</p><p>2 + p3x+ p4 = 0</p><p>with</p><p>p1 =ab0lr</p><p>2</p><p>k2,</p><p>p2 =rk</p><p>[</p><p>(b0 d0)l ab0r</p><p>(</p><p>2l + abEk</p><p>)]</p><p>,</p><p>p3 = r</p><p>[</p><p>2abb0Ek</p><p>(ar ) + abd0Erk</p><p>+ b0l(ar )</p><p>+ d0l +aqE</p><p>k</p><p>]</p><p>and</p><p>p4 = E</p><p>[</p><p> abb0r(ar 2) q(ar )</p><p> bd0(ar ) 2bb0</p><p>]</p><p>.</p><p>Observe that p1 is always positive and p4 is negative if ar &gt;2. Thus, the above cubic equation has at least one positiveroot when ar &gt; 2. The equilibrium predator density isgiven by</p><p>y =</p><p>rx(</p><p>1 x</p><p>k</p><p>)</p><p> ar</p><p>(</p><p>1 x</p><p>k</p><p>) .</p><p>Note that y will be feasible if kar</p><p>(ar ) &lt; x &lt; k.To observe the number of feasible coexistence equilibriumpoints, we plot the prey and predator isoclines of the system(1) in the Fig. 1. This figure confirms that there existstwo positive equilibrium points E(93.9826, 57.5527) andE1(99.3691, 0.9843) for the parameter set given in the Table1. Eigenvalues evaluated at these equilibrium points showthat the equilibrium point E has two negative eigenvaluesand thus stable. The eigenvalues corresponding to the equi-librium point E1 are of opposite signs and therefore unstable.The rest of the paper deals with the stable equilibrium pointE(x, y).</p><p>Let X(t) = x(t) x and Y (t) = y(t) y are the per-turbed variables. Linearizing the system (2) at E(x, y),we get the linear system as follows:</p><p>dXdt</p><p>= xyX</p><p>(ay+x)2 x2Y</p><p>(ay+x)2 rxX(t)</p><p>k,</p><p>dYdt</p><p>= ab0y2X</p><p>(ay+x)2 +[</p><p>qEly</p><p>(bE+ly)2 ab0x</p><p>y</p><p>(ay+x)2</p><p>]</p><p>Y.(3)</p><p>The characteristic equation of the corresponding variationalmatrix is given by</p><p>2 +A+B + (C +D)e = 0, (4)</p><p>20 40 60 80 100 1200</p><p>20</p><p>40</p><p>60</p><p>x</p><p>y</p><p>E(99.3691, 0.9843)</p><p>1 = 0.9878</p><p> 2= 0.1608</p><p>Predator Isocline</p><p>Prey Isocline</p><p> 1= 0.931</p><p>2= 0.0342</p><p>E*(93.9826, 57.5527)</p><p>Fig. 1. The prey and predator isoclines of the model system (1). Thereexists two equilibrium points (shown by bullets), where prey and predatorisoclines intersect. Equilibrium values and the corresponding eigenvaluesare also mentioned. Parameters are as in the TABLE 1.</p><p>where</p><p>A = [</p><p>(1ab0)xy</p><p>(ay+x)2 +qEly</p><p>(bE+ly)2</p><p>]</p><p>,</p><p>B = xy</p><p>(ay+x)2</p><p>[</p><p>qEly</p><p>(bE+ly)2 ab0x</p><p>y</p><p>(ay+x)2</p><p>]</p><p>+ 2ab0x</p><p>2y2</p><p>(ay+x)4 ,</p><p>C = rx</p><p>kand</p><p>D = rx</p><p>k</p><p>[</p><p>qEly</p><p>(bE+ly)2 ab0x</p><p>y</p><p>(ay+x)2</p><p>]</p><p>.</p><p>We now discuss two following cases.</p><p>Case 1: Instantaneous feedback mechanism ( = 0)</p><p>In the case of instantaneous feedback response of preypopulation, the value of becomes zero. In this case, theequation (4) becomes</p><p>2 + (A+ C) + (B +D) = 0. (5)</p><p>All roots of the equation (5) have negative real parts if andonly if</p><p>(H1) A+ C &gt; 0 and B +D &gt; 0.</p><p>Therefore, the equilibrium point E(x, y) is locallyasymptotically stable when (H1) holds. Then the followingtheorem is true.</p><p>Theorem 1. The coexistence equilibrium E(x, y) of thesystem (1) is locally asymptotically stable in absence ofdelay if H1 holds.</p><p>Case 2: Delayed feedback mechanism ( 6= 0)</p><p>We first reproduce some definitions given by [21], [22].</p><p>Definition 1. The equilibrium E is called asymptoticallystable if there exists a &gt; 0 such that</p><p>sup0[ 1() x , 2() y</p><p> ] &lt; </p><p>Engineering Letters, 20:3, EL_20_3_05</p><p>(Revised online publication: 30 August 2012)</p><p> ______________________________________________________________________________________ </p></li><li><p>implies that</p><p>limt(x(t), y(t)) = (x, y),</p><p>where (x(t), y(t)) is the solution of the system (2) whichsatisfies the prescribed initial condition.</p><p>Definition 2. The equilibrium E is called absolutelystable if it is asymptotically stable for all delays 0 andconditionally stable if it is stable for in some finite interval.</p><p>For the delay-induced system (2), the interior equilibriumE(x, y) will be asymptotically stable if all the roots of thecorresponding characteristic equation (4) have negative realparts. But, the classical Routh-Hurwitz criterion cannot beused to discuss stability of the system since equation (4) is atranscendental equation and has infinitely many eigenvalues.To determine the nature of the stability, we require the signof the real parts of the roots of the equation (4). We startwith the assumption that E(x, y) is asymptotically stablein the case of non-delayed system (Case 1) and then we findconditions for which E(x, y) is still stable in presenceof delay [23]. By Rouches Theorem [24] and the continuityin , the transcendental equation (4) has roots with positivereal parts if and only if it has purely imaginary roots. Fromthis, we shall be able to find conditions for all eigenvaluesto have negative real parts.</p><p>Let() = () + i(),</p><p>where and are real. As the interior equilibriumE(x, y) of the non-delayed system is stable, we have(0) &lt; 0. By continuity, if (&gt; 0) is sufficiently small,we still have () &lt; 0 and E(x, y) is still stable. Thechange of stability will occur at some values of for which() = 0, () 6= 0, i.e., will be purely imaginary.Let 0 be such that (0) = 0 and (0) = 0 6= 0,so that = i(0) = i0. In this case, the steady stateloses stability and eventually becomes unstable when ()becomes positive. Now i0 is a root of (4) if and only if</p><p>02 + iA0 +B + (iC0 +D)e</p><p>00 = 0.</p><p>Equating the real and imaginary parts of both sides, we get</p><p>Dcos00 + C0sin00 = 02 B,C0cos00 Dsin00 = A0.</p><p>(6)</p><p>This leads to</p><p>04 [C2 A2 + 2B]0</p><p>2 +B2 D2 = 0. (7)</p><p>This equation has no positive roots if the following condi-tions are satisfied:</p><p>(H2) A2 C2 2B &gt; 0</p><p>and B2 D2 &gt; 0.</p><p>Hence, all roots of the equation (7) will have negative realparts when [0,) if conditions of the Theorem 1 and(H2) are satisfied.</p><p>Let(H3) B</p><p>2 D2 &lt; 0.</p><p>If Theorem 1 and (H3) hold then the equation (7) has aunique positive root 20. Substituting </p><p>20 into (6), we have</p><p>n =1</p><p>0cos1</p><p>[</p><p>(D AC)20 BD</p><p>D2 + C220</p><p>]</p><p>+2n</p><p>0, n = 0, 1, 2, ..</p><p>Let,</p><p>(H4) C2 A2 + 2B &gt; 0,</p><p>B2 D2 &gt; 0 and(C2 A2 + 2B)2 &gt; 4(B2 D2).</p><p>If Theorem 1 and (H4) hold then (7) has two positive roots2+ and </p><p>2. Substituting into (7), we get</p><p>2k =1</p><p>cos1</p><p>[</p><p>(D AC)2 BD</p><p>D2 + C22</p><p>]</p><p>+2k</p><p>, k = 0, 1, 2, ..</p><p>If () be the root of (4) satisfying Re(n) = 0 (respec-tively, Re(k ) = 0) and Im(n) = 0 (respectively,Im(k ) = ), we get[</p><p>ddRe()</p><p>]</p><p>=0,=0</p><p>= (A2C22B)0</p><p>2+24</p><p>C202+D2&gt; 0. (by (H2))</p><p>Similarly, we can show that[</p><p>ddRe()</p><p>]</p><p>=+k</p><p>,=+&gt; 0 and</p><p>[</p><p>ddRe()</p><p>]</p><p>=k</p><p>,=&lt; 0.</p><p>From Corollary (2.4) in Ruan and Wei [25], we have thefollowing conclusions.</p><p>Theorem 2. Assume 6= 0 and conditions of the Theorem1 are satisfied. Then the following conclusions hold:</p><p>(i) If H2 hold then the equilibrium E(x, y) is asymp-totically stable for all 0.</p><p>(ii) If H3 hold then the equilibrium E(x, y) is condi-tionally stable. It is locally asymptotically stable for &lt; 0and unstable for &gt; 0. Furthermore, the system (1)undergoes a Hopf bifurcation at E(x, y) when = 0,where</p><p>0 =1</p><p>0cos1</p><p>[</p><p>(D AC)20 BD</p><p>D2 + C220</p><p>]</p><p>.</p><p>(iii) If H4 holds then there is a positive integerm such that the equilibrium is stable when [0, +0 )(</p><p>0 , </p><p>+1 ).....(</p><p>m1 , </p><p>+m), and unstable when</p><p> [+0 , 0 ) (</p><p>+1 , </p><p>1 ) ..... (</p><p>+m1, </p><p>m1) (</p><p>+m,).</p><p>Further more, the system undergoes a Hopf bifurcation atE(x, y) when = +k , k = 0, 1, 2, ....</p><p>III. DIRECTION AND STABILITY OF HOPF BIFURCATION</p><p>In the previous section, we obtain the conditions underwhich a family of periodic solutions bifurcate from the co-existence equilibrium E(x, y) at the critical value = 0.As pointed out in Hassard et al. [26], it is interested todetermine the direction, stability and period of the periodicsolutions bifurcating from the positive equilibrium. Follow-ing the ideas of Hassard et al., we derive the explicit formulaefor determining the properties of the Hopf bifurcation at thecritical value 0 by using the normal form and the centermanifold theory. Throughout this section, we always assumethat the system (2) undergoes a Hopf bifurcation at thepositive equilibrium E(x, y) for = 0 and then i0 isthe corresponding purely imaginary roots of the characteristicequation at the positive equilibrium E(x, y).</p><p>Let x1 = xx, x2 = y y, xi = xi(t), = 0 +.Dropping the bars for simplification of notations, system (2)</p><p>Engineering Letters, 20:3, EL_20_3_05</p><p>(Revised online publication: 30 August 2012)</p><p> ______________________________________________________________________________________ </p></li><li><p>is transformed into an functional differential equation (FDE)in C = C([1, 0], R2) as</p><p>x(t) = L(xt) + f(, xt), (8)</p><p>where x(t) = (x1, x2)T R2 and L : C R, f : R C R are given by</p><p>L() = (0 + )</p><p>(</p><p>a11 a12a21 a22</p><p>)(</p><p>1(0)2(0)</p><p>)</p><p>+ (0 + )</p><p>(</p><p>b11 00 0</p><p>)(</p><p>1(1)2(1)</p><p>) (9)</p><p>and</p><p>f(, ) = (0 + )</p><p>(</p><p>c1c2</p><p>)</p><p>, (10)</p><p>where,</p><p>a11 =xy</p><p>(ay+x)2 ,</p><p>a12 = x2</p><p>(ay+x)2 ,</p><p>a21 =ab0y</p><p>2</p><p>(ay+x)2 ,</p><p>a22 =qEly</p><p>(bE+ly)2 ab0x</p><p>y</p><p>(ay+x)2 ,</p><p>b11 = rx</p><p>k,</p><p>c1 = r1(0)1(1)</p><p>k 1(0)2(0)</p><p>a2(0)+1(0)and</p><p>c2 =b01(0)2(0)a2(0)+1(0)</p><p> ql2(0)2</p><p>b2E.</p><p>By Riesz representation theorem, there exits a function(, ) of bounded variation for [1, 0], such that</p><p>L = 0</p><p>1 d(, )() for C. (11)</p><p>In fact, we can choose</p><p>(, ) = (0 + )</p><p>(</p><p>a11 a12a21 a22</p><p>)</p><p>()</p><p> (0 + )</p><p>(</p><p>b11 00 0</p><p>)</p><p>( + 1),(12)</p><p>where is the Dirac delta function. For C1([1, 0], R2),define</p><p>A() =</p><p>d()</p><p>d, [1, 0),</p><p> 0</p><p>1 d(, s)(s), = 0</p><p>and</p><p>R() =</p><p>{</p><p>0, [1, 0)f(, ), = 0.</p><p>Then system (8) is equivalent to</p><p>x(t) = A()xt +R()xt, (13)</p><p>where xt() = x(t + ) for [1, 0]. For C1([0,1], (R2)), define</p><p>A(s) =</p><p>d(s)</p><p>ds, s (0, 1]</p><p> 0</p><p>1dT (t, 0)(t), s = 0</p><p>and a bilinear inner product</p><p>(s), () = (0)(0) </p><p> 0</p><p>1</p><p>=0</p><p>( )d()()d,</p><p>(14)where () = (, 0). Then A(0) and A are adjoint opera-tors. By the discussion in section II, we know that i00are eigenvalues of A(0). Thus, they are also eigenvaluesof A. We first need to compute the eigenvalues of A(0)</p><p>and A corresponding to i00 and i00, respectively.Suppose that q() = (1, )T ei00 is the eigenvector ofA(0) corresponding to i00, then A(0)q() = i00q(). Itfollows fro...</p></li></ul>


View more >