Stability and Bifurcation of a Prey-Predator System with Additional Food and Two Discrete Delays

Ankit Kumar and Balram Dubey

Department of Mathematics,BITS Pilani,Pilani Campus,Pilani,333031,India

ABSTRACT In this paper,the impact of additional food and two discrete delays on the dynamics of a prey-predator model is investigated.The interaction between prey and predator is considered as Holling Type-II functional response.The additional food is provided to the predator to reduce its dependency on the prey.One delay is the gestation delay in predator while the other delay is the delay in supplying the additional food topredators.The positivity,boundedness and persistence of the solutions of the system are studied to show the system as biologically well-behaved.The existence of steady states,their local and global asymptotic behavior for the non-delayed system are investigated.It is shown that(i)predator’s dependency factor on additional food induces a periodic solution in the system,and(ii)the two delays considered in the system are capable to change the status of the stability behavior of the system.The existence of periodic solutions via Hopf-bifurcation is shown with respect to both the delays.Our analysis shows that both delay parameters play an important role in governing the dynamics of the system.The direction and stability of Hopf-bifurcation are also investigated through the normal form theory and the center manifold theorem.Numerical experiments are also conducted to validate the theoretical results.

KEYWORDS Additional food;Gestation delay;Hopf-bifurcation;Prey-predator

1 Introduction

Living organisms on the surface of the earth adopt only the way that boosts their survival possibilities so that they can pass their genes to the next generation.There are several fundamental instincts in ecological communities and,predation is one of them that constitutes the building blocks for multispecies food webs.Initially,Lotka [1] and Volterra [2] studied the model for preypredator interaction and observed the uniform fluctuations in the time series of the system.On later the fluctuations were removed from the system by taking logistic growth of prey population [3,4].Many researchers have widely studied prey-predator interactions for the last century[5-12].They have considered several essential concepts over time that play a vital role in the dynamics of the system like functional response,time delay,harvesting and conservation policies of species,stage structure,fear induced by predators,etc.The idea of functional response was proposed by Holling [5].It is defined as the consumption rate of prey by predators.Holling considered it nonlinear function of prey species that saturates at a level.Further,it was considered a function of prey and predator both by several authors [6,10,13,14].

In last few decades,many authors have studied the qualitative dynamics of prey-predator systems in the presence of additional food resources for predators [15-20].Additional food is an important component for predators like coccinellid which shapes the life history of many predator species [17].Ghosh et al.[20] investigated the impact of additional food for predator on the dynamics of prey-predator model with prey refuge and they observed that predator extinction possibility in high prey refuge may be removed by providing additional food to predators.Again,to study the role of additional food in an eco-epidemiological system,a model was proposed and studied by Sahoo [21].The author found that the system becomes disease free in presence of suitable additional food provided to predator.Recently,a prey-predator model with harvesting and additional food is analyzed by Rani et al.[22] and they have shown some local and global bifurcation with respect to different parameters.To incorporate the additional food into the model,they modified the Holling type II functional response.

Delayed models exhibit much more realistic dynamics than non delayed models [4,23].In preypredator system,the impact of consumed prey individuals into predator population does not appear immediately after the predation,there is some time lag that is gestation delay [24].We incorporate the effect of time delay into the model with delay differential equations.A delay differential equation demonstrates much more complex character than ordinary differential equation.On the other hand predators do not consume the additional food as soon as it is provided.They take some time to consume and digest the food.Delayed models are widely studied by researchers [10,14,25-33].A delayed prey and predator density dependent system is investigated by Li et al.[10].The authors analyzed stability,Hopf-bifurcation and its qualititative properties by using Poincare normal form and the formulae given in Hassard et al.[34].Sahoo et al.[35]examined prey-predator model with effects of supplying additional food to predators in a gestation delay induced prey-predator system and habitat complexity.They have pointed out that Hopfbifurcation occurs in the system when delay crosses a threshold value that strongly depends on quality and quantity of supplied additional food.The effect of additional food along with fear induced by predators and gestation delay is discussed by Mondal et al.[36].There are several studies carried out with multiple delays [37-40].Li et al.[37] have done stability and Hopf-bifurcation analysis of a prey-predator model with two maturation delays.Gakkhar et al.[30] explored the complex dynamics of a prey-predator system with multiple delays.They established the presence of periodic orbits via Hopf-bifurcation with respect to both delays.Recently,Kundu et al.[40]have discussed about the dynamics of two prey and one predator system with cooperation among preys against predators incorporating three discrete delays.The authors have found that all delays are capable to destabilize the system.

To the best of our knowledge,an ecological model including (i) Effect of additional food supplies to predators,(ii) Dependency factor of supplied additional food,(iii) Holling Type II functional response,(iv) Gestation delay in predator have not been considered.Inspired by this,we establish three dimensional non delayed and delayed models in Section 2.We analyze the dynamics of non delayed model and validate it via some numerical simulations in Section 3.In Section 4,we analyze the dynamics of delayed model through Hopf-bifurcation.Direction and stability of Hopf-bifurcation are carried out in Section 5.Section 6 is devoted to the numerical simulations for delayed model.Conclusions and significance of this work are discussed briefly in Section 7.

2 Proposed Mathematical Model

We consider a habitat where two biological populations,prey population and predator population are surviving and interacting with each other.It is assumed that prey population grows logistically and the interaction between prey and predator follows Holling Type II functional response.We assume that the density of the additional food supplied to the predators is directly proportional to the density of predators present in the habitat.Keeping these in view,the dynamics of the system can be governed by the following system of differential equations:

In the above modelx(t),y(t)are number of prey and predator individuals at timetandAis quantity of additional food provided to predators.A0is dependency factor of predators on provided additional food resources.IfA0=1,then predators depend completely on additional food and prey population grows logistically.IfA0=0,then predators depend only on the prey population and in such a case additional food is not required.λis maximum supply rate of additional food resources.

In real situations,each organism needs an amount of time to reproduce their progeny.Due to this fact the increment in predators does not appear immediately after consuming prey.It is assumed that a predator individual takesτ1time for gestation.Therefore,it seems reasonable to incorporate a gestation delay in the system.Thus,the delayτ1is considered in the numeric response only.Again it is assumed that the additional food is provided to predators with another delayτ2.The generalized model involving these two discrete delays takes the following form

subject to the non negative conditionsx(s)=φ1(s)≥0,A(s)=φ2(s)≥0,y(s)=φ3(s)≥0,s∈[?τ,0],whereτ=max{τ1,τ2}andφi(s)∈C([?τ,0]→R+),(i=1,2,3).

The biological meaning of all parameters and variables in above models is provided in Tab.1.

3 Dynamics of Non-delayed Model

First of all,we examine the boundedness and persistence of the system (1).

3.1 Boundedness and Persistence of the Solution

Theorem 3.1.The set

is a positive invariant set for all the solutions of model (1),initiating in the interior of the positive octant,whereδ=min{r,β,d}.

ProofThe model system (1) can be written in the matrix form

whereX=(x1,x2,x3)T=(x,A,y)T∈R3,andG(X)is given by

SinceG:R3→R3+is locally Lipschitz-continuous inΩandX(0)=X0∈R3+,the fundamental theorem of ordinary differential equation guarantees the local existence and uniqueness of the solution.Since [Gi(X)]xi(t)=0,X∈R3+≥0,it follows thatX(t)≥0 for allt≥0.In fact,from the first equation of model (1),it can easily be seen that|x=0≥0,|y=0≥0 and hencex(t)≥0,y(t)≥0 for allt≥0.Secondly,|A=0=λA0y≥0 for allt≥0 (asy(t)≥0 for allt≥0.) and henceA(t)≥0 for allt≥0.

From the first equation of model (1),we can write

which yields

Now,suppose

W(t)=c1x(t)+c2A(t)+y(t),

then we have

whereδ=min{r,β,d}.

Hence,it follows that

We also note that ifx(t)≥KandthenThis shows that all solutions of system (1) are bounded and remains inΩfor allt＞0 if(x(0),A(0),y(0))∈Ω.

Theorem 3.2.Let the following inequalities are satisfied:

Then model system (1) is uniformly persistence,where,xais defined in the proof.

Proof:System (1) is said to be permanence or uniform persistence if there are positive constantsM1andM2such that each positive solutionX(t)=(x(t),A(t),y(t))of the system with positive initial conditions satisfies

Keeping the above in view,if we define

then from Theorem 3.1,it follows that

This also shows that for any sufficiently smallε＞0,there exists aT＞0 such that for allt≥T,the following holds:

Now from the first equation of model system (1),for allt≥T,we can write

Hence,it follows that

which is true for everyε＞0,thus

wherer＞α(1?A0)ys.

Now from the third equation of model system (1),we obtain

which implies

which is true for everyε＞0,thus

for persistence,we must have

Second equation of model system (1) yields

Hence

TakingM1=min{xa,Aa,ya},the theorem follows.

Remark.Theorem 3.2 shows that threshold values for the persistence of the system are dependent on the parameterA0.

3.2 Equilibrium Points and Their Stability Behavior

System (1) has four equilibrium points,trivial equilibriumE0(0,0,0),axial equilibriumE1(K,0,0),prey free equilibriumE2(0,,) and interior equilibriumE?(x?,A?,y?).E0andE1always exist.

? Existence ofE2(0,,):The prey free equilibriumE2is positive solution of the following system:

From the second equation of above system,we have

Putting the value ofAin the first equation of system (3),we get

Above equation has two positive roots if

System (1) has two prey free equilibrium under conditions given in (5):andAgain,Ifc2φλA0<φd+βe,then Eq.(4) does not have any positive root.Therefore,E2does not exist in this case.

? Existence of interior equilibriumE?(x?,A?,y?):It may be seen thatx?,A?andy?are the positive solution of the following system of algebraic equations:

From the second equation of system (6),we have

Putting this into the first and third equation of system (6),we obtain the following system:

We note the following points from Eq.(7):

1.Wheny=0,thenx=<0 orx=K＞0.

2.Whenx=0 theny=＞0.

Similarly,from Eq.(8),we note the following:

1.Wheny=0,thenx=

Hence,we can state the following theorem.

Theorem 3.3.The system (1) has a unique positive equilibriumE?(x?,A?,y?) if (9)and (10) hold.

Remark.The number of positive equilibrium for the system (1) depends on values of parameters,which we have chosen.Several possibilities are depicted in Fig.1.

The local behavior of a system in the vicinity of any existing equilibrium is very close to the behavior of its Jacobian system.So,we compute the Jacobian matrix to see the local behavior of the system around its equilibrium and we observe that

? The trivial equilibriumE0(0,0,0) is always a saddle point having stable manifold along the

Aandy-axes and unstable manifold along thex-axis.

? The axial equilibriumE1(K,0,0) is locally asymptotically stable ifIfthenE1is a saddle point having stable manifold along thexandA-axes and unstable manifold along they-axis.

? The Jacobian matrix evaluated at prey free equilibriumis given by

Characteristic equation is given by

The roots of Eq.(11) have negative real part if

Eq.(11) have at least one positive and one negative root if

Remark.By replacingbyandby,similar analysis holds good for the stability behavior of

Figure 1:Four possibilities of the prey and predator zero growth isoclines.(a) Interior equilibrium does not exist for the parametric values a=0.08,d=0.01,(b) interior equilibrium exists uniquely for the values of parameters a=0.1,d=0.235,(c) two interior equilibria for parameter values a=0.1,d=0.137,(d) three interior equilibria for parameter values a=0.105,d=0.1.Rest of the parameters are same as that in (25)

? In order to analyze the local stability of unique interior equilibriumE?(x?,A?,y?),we evaluate the Jacobian matrix atE?and it is given by

Characteristic equation corresponding to above matrix is given by

where

Now using the Routh-Hurwitz criterion,all eigenvalues ofJ|E?have negative real part iff

Thus we can state the following theorem.

Theorem 3.4The system (1) is stable in the neighborhood of its positive equilibrium iff inequalities in (15) hold.

It is also noted that inequalities in (15) hold if

Infect,the above two conditions imply thatA1＞0 andA3＞0.The third conditionA1A2＞A3is also satisfied.

Remark.The system (1) is stable around its positive equilibriumE?if inequalities in (16) hold.

In the following theorem we give a criterion for global asymptotic stability of interior equilibriumE?(x?,A?,y?) of the system (1).

Theorem 3.5.The interior equilibriumE?(x?,A?,y?) of the system (1) is globally asymptotically stable under the following conditions:

Proof.We Choose a suitable Lyapunov function aboutE?as

whereγ1andγ2are positive constants,to be specified later.Now,differentiatingVwith respect totalong the solutions of system (1),we get

Choosingγ2=andγ1=we get

Applying Sylvester criterion,is negative definite if conditions in (17) hold.HenceE?is globally stable under conditions in (17).

3.3 Hopf-Bifurcation and Its Properties

Hopf-bifurcation is a local phenomenon where a system’s stability switches and a periodic solution arises around its equilibrium point by varying a parameter.In system (1),the parameterA0seems crucial,therefore we analyze the Hopf-bifurcation by takingA0as bifurcation parameter,then we have someA0=A?0.The necessary and sufficient conditions for occurrence Hopf-bifurcation atA0=A?0are

(a)A1|A?0＞0,A3|A?0＞0,

(b)f(A?0)≡(A1A2?A3)|A?0=0,

FromA1A2?A3=0,we get an equation inA0and assume that it has at least one positive rootA?0.Then for someε＞0,there is an interval containingA?0,(A?0?ε,A?0+ε)such thatA?0?ε＞0 andA2＞0 forA0∈(A?0?ε,A?0+ε).Thus,Eq.(14) cannot have any real positive root forA0∈(A?0?ε,A?0+ε).

Therefore,atA0=A?0,Eq.(14) becomes

(Λ+A1)(Λ2+A2)=0,

this gives us three roots

Λ1,2=±iρ,Λ3=μ,

whereρ=andμ=?A1.

ForA0∈(A?0?ε,A?0+ε),roots can be taken as

Λ1,2=k1(A0)±ik2(A0),Λ3=?A1(A0).

Now,we have to verify the transversality condition.Differentiating Eq.(14) with respect to the bifurcation parameterA0,we obtain

whereR=A1A2?A3and,i=1,2,3 denote the derivative ofAiwith respect to time.Thus

Thus,we can state the following theorem.

Theorem 3.6.The system undergoes Hopf-bifurcation near interior equilibriumE?(x?,A?,y?)under the necessary and sufficient conditions (a),(b) and (c).Critical value of bifurcation parameterA0is given by the equationf(A?0)=0.

In order to see the stability and direction of Hopf-bifurcation,we use center manifold theorem [34] and some concepts used in [41].Now,consider the following transformation

x1=x?x?,x2=A?A?,x3=y?y?.

Using this transformation,system (1) takes the following form

whereX=(x1,x2,x3)T,

Letv1andv2be the eigenvectors corresponding to eigenvaluesiρandμofE?atA0=A?0.Thenv1andv2are given by

and

Let

where

Δ=p11(p22p33?p23p32)+p12(p23p31?p21p33)+p13(p21p32?p22p31)/=0,

q11=p22p33?p23p32,q12=p21p33?p23p31,q13=p21p32?p22p31,

q21=p12p33?p13p32,q22=p11p33?p13p31,q23=p11p32?p31p12,

q31=p12p23?p13p22,q32=p11p23?p13p21,q33=p11p22?p12p21.

Now letX=HYorY=H?1X,whereY= (y2,y2,y3)T.Using this transformation,system (18) can be written as

where

So,we can write system (19) as

Thus,system (20) takes the following form

whereU=(y1,y2)T,V=(y3),C=(μ),f=(f1,f2) andg=(f3).The eigenvalues ofBandCmay have zero real part and negative real parts,respectively.f,gvanish along with their first partial derivative at the origin.

Since the center manifold is tangent toWC(y=0) we can represent it as a graph

WC={(U,V):V=h(U)}:h(0)=h′(0)=0,

whereh:U→R2is defined on some vicinityU?R2of the origin [42,43].

We consider the projection of vector field onV=h(U)ontoWC:

Now we state the following theorem to approximate the center manifold.

Theorem 3.7.LetΦbe aC1mapping of a neighborhood of the origin inR2intoRwithΦ(0)=0 andΦ′(0)=0.If for someq＞1,(NΦ)(U)=o(|U|q) asU→0,thenh(U)=Φ(U)+o(|U|q) asU→0,where(NΦ)(U)=Φ′(U)[BU+f(U,Φ(U))]?CΦ(U)?g(U,Φ(U)).

In order to approximateh(U),we consider

where h.o.t.stands for high order terms.Using (23),we get from (22)After simplification,we get

where

Equating both the sides of Eq.(24),we get

Using Crammer’s rule,

We can find the behavior of the solution of system (21) from the following theorem.

Theorem 3.8.If the zero solution of (22) is stable (asymptotically stable/unstable),then the zero solution of (21) is also stable (asymptotically stable/unstable).

Now from Eq.(22),we have

where

We determine the direction and stability of bifurcation periodic orbit of the system (21) by the following formula [44]

In above expression,ifν＞0(<0),then the Hopf-bifurcation is supercritical (subcritical)and bifurcation periodic solution exists forA0=A?0.The bifurcating periodic solution is stable(unstable) if

The bifurcating direction of periodic solution of the system (1) is same as the system (21).

3.4 Numerical Simulation

To validate our theoretical findings of model (1),we perform some numerical simulations using MATLAB R2018b.We have chosen the following dataset

For the above set of parameters,condition for existence of prey free equilibrium (5) and conditions for existence and uniqueness of interior equilibrium (9) and (10) are satisfied.Therefore,the system (1) has five equilibrium points.The behavior of these equilibrium points are given in Tab.2.

Table 2:Existing equilibria and their stability nature

The eigenvalues of the Jacobian matrix atE0andE1are (3,?0.32,?0.3) and(?3,?0.32,0.4113),respectively.ThereforeE0andE1both are saddle points.SimilarlyE2and ?E2are also saddle points.Again all the inequalities in (15) are satisfied.So according to Theorem 3.4,the interior equilibriumE?is locally asymptotically stable.The stability of system in the vicinity of the positive equilibriumE?is illustrated by Fig.2.In Fig.2a,time evolution of species is shown and it is noted that they converge to their equilibrium levels after some oscillations.In Fig.2b,phase diagram is drawn inxAy-space which shows the asymptotic stability behavior of positive equilibriumE?.

In this study,we found that predators dependency factorA0on additional food plays an important role in the dynamics of the system.If it is less than a threshold value then it can be the cause of destabilizing the system.The threshold value can be calculated by solvingf(A?0)=0(Theorem 3.6).By our computer simulation we obtain it asA?0=0.482.All the conditions of Theorem 3.6 are satisfied,so the system undergoes a Hopf-bifurcation atA?0=0.482.If we keep the value of parameterA0below its threshold value,then the system (1) always remains unstable.The instable behavior of solutions and presence of stable limit cycle atA0=0.45

In the model (1),consumption rate of additional foodφis also a vital parameter.We have noted that if system is stable for parameterA0(A0∈[0.482,1]) then it is stable for all range of parameterφ.But ifA0∈[0.2361,0.482) then system undergoes a Hopf-bifurcation with respect to parameterφ.In Fig.5,we have shown the bifurcation diagram whenA0=0.4 and other parameters are same as given in (25).The Hopf-bifurcation point isφ?=0.02847.

Figure 2:Time series evolution (a) and phase portrait (b) of species for the set of parameters chosen in (25).Positive equilibrium E?is locally asympeoeically stable

As the system (1) shows Hopf-bifurcation with respect to parametersA0andφ,and direction of Hopf-bifurcation is opposite for both the parameters.Therefore,we can divide theA0φ-plane into two regions

Region of stability (green)S1={(A0,φ):system (1) is locally asymptotically stable},

Region of instability (white)S2={(A0,φ):system (1) is unstable}.

Both the regions are drawn in Fig.6.The curve which separates both the regions is called Hopf-bifurcation curve.

The number of interior equilibrium points depend on the values of parameters.In the table below,we have shown dependence of total number of interior equilibrium on parametersaanddand the nature of their stability.It is observed that whena=0.105 andd=0.1 (other parameters are as in (25)),then three interior equilibrium exist for the system (1),E?1(0.5099,1.398,20.9174),E?2(8.2355,1.409,32.9068) andE?3(42.309,1.4136,43.0591).E?1andE?3are locally asymptotically stable andE?2is unstable.Since there are two locally asymptotically stable equilibrium in the system,so it shows bistability.Bistability is a phenomenon where a system converges to two different equilibrium points for the same parametric values based on the variation of the initial conditions.In Fig.7,we initiated two trajectories from two nearby points and they converse to different interior equilibria.The black dotted curve is separatrix,which divides thexy-plane into two regions in such a way that if a solution is initiated from the left of the separatrix,it converses toE?1and if a solution is initiated from the right of the separatrix,it converses toE?3.In other words,left region is region of attraction forE?1and right region is region of attraction forE?3.

Remark.For the best representation of bistability phenomenon and separatrix curve,the Fig.7 is drawn in thexy-plane.But initial conditions and interior equilibrium points are written as they are.

Figure 3:Instable behavior of solutions and existence of stable limit cycle for A0=0.45(

4 Analysis of Delayed Model

In this section,we discuss the local stability and Hopf-bifurcation phenomenon for the delayed system (2).The introduction of time delay does not affects the equilibria of the system.So,all the equilibria remain same as the non-delayed system (1).To see the effect of delay on the dynamical behavior of the interior equilibriumE?,we rewrite the delayed system (2) as

where

U(t)=[x(t),A(t),y(t)]T,U(t?τ1)=[x(t?τ1),A(t?τ1),y(t?τ1)]T,

U(t?τ2)=[x(t?τ2),A(t?τ2),y(t?τ2)]T.

Figure 4:Bifurcation diagram of the prey and predator population with respect to parameter A0

Figure 5:Bifurcation diagram of the prey and predator population with respect to parameter φ

Now we linearize the system (26) by using the following transformations:

Figure 6:Region of stability and instability for system (1) in A0φ-plane

Table 3:Dependence of total number of interior equilibria and their stability on parameters a and d.Rest of the parameters are same as in (25)

where

Thus,the Jacobian matrix of the system (2) atE?is given by

Figure 7:Trajectories initiated from region of attraction of both the locally asymptotically stable equilibrium points,system (1) shows bistability

where

The characteristic equation corresponding to the above Jacobian matrix is

where

b1=?(a1+a3+a6),b2=a3a6+a1a6+a1a3,b3=?a1a3a6,b4=?c1a2,

b5=c1a2(a1+a3+a5),b6=?c1a2a3(a1+a2),b7=?c2φA?,b8=c2φ(a1A?+a3A??a4y?),b9=c2a1φ(?a3A?+a4y?).

Remark.Whenτ1=τ2=0,then the characteristic Eq.(27) is same as the characteristic Eq.(14) for non-delayed system.

Case (1):τ1＞0,τ2=0.Then Eq.(27) becomes

where

d1=b1+b7,d2=b2+b8,d3=b3+b9.

For the delayed system (2),the positive equilibrium is locally asymptotically stable if and only if all the roots of the Eq.(28) have negative real parts.For switching of the stability,the root of the Eq.(28) must cross the imaginary axis.Therefore letiω(ω＞0) be a root of Eq.(28),then it follows that

From the above set of equations,we can obtain

where

h1=d21?b24?2d2,h2=d22?b25?2d1d3+2b4b6,h3=d23?b26.

If we putω2=z,then Eq.(30) becomes

Theorem 4.1.If Eq.(31) has no positive root,then there is no change in the stability behavior ofE?for allτ1≥0.

Corollary.If inequalities in (15) hold and Eq.(31) has no positive root,thenE?is locally asymptotically stable for allτ1≥0.

Corollary.If inequalities in (15) do not hold and Eq.(31) has no positive root,thenE?is unstable for allτ1≥0.

Now let inequalities in (15) hold and Eq.(31) has at least one positive root,sayz1=ω21.Substitutingω1into Eq.(29),we obtain

(H1):g′(ω21)＞0.

Letξ(τ1i)=±iω1be the root of Eq.(28),a little calculation yields

Hence,the transversality condition can be obtained under (H1)

Thus,we can state the following theorem.

Theorem 4.2.For system (2),withτ2=0 and assuming that (H1) holds,there exists a positive numberτ10such that the equilibriumE?is locally asymptotically stable whenτ1<τ10and unstable whenτ1＞τ10.Furthermore system (2) undergoes a Hopf-bifurcation atE?whenτ1=τ10.

Case (2):τ1=0,τ2＞0.Then Eq.(27) becomes

where

e1=b1+b4,e2=b2+b5,e3=b3+b6.

Under an analysis similar to Case (1),one can easily deduce the following theorem.

Theorem 4.3.Forτ1=0,the interior equilibrium point is locally asymptotically stable forτ2<τ20,unstable forτ2＞τ20,and it undergoes Hopf-bifurcation atτ2=τ20given by

whereiω2is root of characteristic Eq.(33).

Case(3):τ1is fixed in the interval (0,τ10) and assumingτ2as a variable parameter.We consider Eq.(27) withτ1as fixed in its stable interval (0,τ10) andτ2as a variable.Letiω (ω＞0) be a root of characteristic Eq.(27).Then separating real and imaginary parts,we obtain

Squaring and then adding (34) and (35) to eliminateτ2,we obtain

Eq.(36) is a transcendental equation in complex form.So,it is not easy to predict the nature of roots.Without going detailed analysis with (36),it is assumed that there exist at least one positive rootω0.Eqs.(34) and (35) can be re-written as

where

D1=?b1ω20+b3+(?b4ω20+b6)cos(ω0τ1)+b5ω0sin(ω0τ1),

D2=?ω30+b2ω0?(?b4ω20+b6)sin(ω0τ1)+b5ω0cos(ω0τ1).

Eqs.(37) and (38) lead to

Now,to verify the transversality condition of Hopf-bifurcation,differentiating equation (34)and (35) with respect toτ2and substituteτ2=τ′20,we obtain

where

(H2):PR?QS/=0.

Theorem 4.4.For system (2),withτ1∈(0,τ10) and assuming that (H2) holds,there exists a positive numberτ′20such thatE?is locally asymptotically stable whenτ2<τ′20and unstable whenτ2＞τ′20.Furthermore,system (2) undergoes a Hopf-bifurcation atE?whereτ2=τ′20.

Case(4):τ2is fixed in the interval (0,τ20) and assumingτ1as a variable parameter Under an analysis similar to Case (3),one can easily prove the following theorem.

Theorem 4.5.Forτ2∈(0,τ20),the interior equilibrium point is locally asymptotically stable forτ1<τ′10and it undergoes Hopf-bifurcation atτ1=τ′10,given by

where

D3=?b1ω2?+b3+(?b7ω2?+b9)cos(ω?τ2)+b8ω?sin(ω?τ2),

D4=?ω3?+b2ω??(?b7ω2?+b9)sin(ω?τ2)+b8ω?cos(ω?τ2),

andiω?is characteristic root of Eq.(27).

5 Direction and Stability of Hopf-Bifurcation

Now with the help of center manifold theory and normal form concept (see [34] for details),we shall study direction and stability of the bifurcated periodic solutions atτ1=τ′10.

Without loss of generality,we assume thatτ?2<τ′10,whereτ?2∈(0,τ20).Let

x1(t)=x(t)?x?,A1(t)=A(t)?A?,y1(t)=y(t)?y?,

and still denotex1(t),A1(t),y1(t)byx(t),A(t),y(t).Letτ1=τ′10+μ,μ∈Rso that Hopfbifurcation occurs atμ=0.We normalize the delay with scalingthen system (2) can be re-written as

whereU(t)=(x(t),A(t),y(t))T,

The nonlinear termf1,f2andf2are given by

The linearization of Eq.(41) around the origin is given by

Forχ=(χ1,χ2,χ3)T∈C([?1,0],R3),define

By the Riesz representation theorem,there exists a 3×3 matrixη(θ,μ),(?1≤θ≤0) whose element are of bounded variation function such that

In fact,we can obtain

Then Eq.(42) is satisfied.

Forχ∈C1([?1,0],R3),define the operatorH(μ)as

and

where

Then system (2) is equivalent to the following operator equation

whereUt=U(t+θ)forθ∈[?1,0].

Forψ∈C1([0,1],(R3)?),define

and a bilinear form

whereη(θ)=η(θ,0),H=H(0) andH?are adjoint operators.From the discussion in previous section,we know that ±are the eigenvalues ofH(0) and therefore they are also eigenvalues ofH?.It is not difficult to verify that the vectors(θ∈[?1,0]) andq?(s)=are the eigenvectors ofH(0) andH?corresponding to the eigenvalueandrespectively,where

Following the algorithms explained in Hassard et al.[34] and using a computation process similar to that in Song et al.[26],which is used to obtain the properties of Hopf-bifurcation,we obtain

Figure 8:Time series evolution and phase portrait of species for the set of parameters in (25) and τ1=0.2<τ10=0.2889 when τ2=0.System is locally asymptotically stable around the positive equilibrium E?

where

Figure 9:System (2) is unstable when τ1=0.35＞τ10=0.2889 and τ2=0.Hopf-bifurcation occurs and stable limit cycle arises in the system

E1= (E1(1),E1(2),E1(3))T∈R3andE2= (E2(1),E2(2),E2(3))T∈R3are constant vectors,computed as:

Figure 10:Bifurcation diagram of the prey and predator population with respect to delay parameter τ1 when τ2=0

Figure 11:Stable time series solutions and phase diagram of system (2) for τ2=0.7<τ20=0.9618 and τ1=0.Other parameters are same as in (25)

Figure 12:Instable behavior and existence of periodic solutions of system (2) around the positive equilibrium E?at τ2=1.2＞τ20=0.9618 and τ1=0

Figure 13:Bifurcation diagram of the prey and predator population with respect to delay parameter τ2 and τ1=0

Figure 14:E?is locally asymptotically stable when τ1=0.12 is fixed in its stable range (0,τ10) and τ2=0.4<τ′20=0.4731

Consequently,gijcan be expressed by the parameters and delaysτ′10andτ?2.Thus,these standard results can be computed as:

These expressions give a description of the bifurcating periodic solution in the center manifold of system (2) at critical valuesτ1=τ10which can be stated in the form of following theorem:

Theorem 5.1

?μ2determines the direction of Hopf-bifurcation.Ifμ2＞0(<0) then the Hopf-bifurcation is supercritical (subcritical).

Figure 15:E?is unstable when τ1=0.12 is fixed in its range of stability (0,τ10) and τ2=0.6＞τ′20=0.4731.Time series solution of species and existence limit cycle

?β2determines the stability of bifurcated periodic solution.Ifβ2＞0(<0) then the bifurcated periodic solutions are unstable (stable).

?T2determines the period of bifurcating periodic solution.The period increases (decreases)ifT2＞0(<0).

Remark.Whenτ1＞0 andτ2=0 orτ1=0 andτ2＞0,then under an analysis similar to Section 5,the corresponding values ofμ2,β2andT2can be computed.Depending upon the sign ofμ2,β2andT2,the corresponding results can also be deduced.

6 Numerical Simulation of Delayed Model

In order to validate our theoretical findings,obtained in previous sections,we perform some simulations by taking the same values of parameters in (25).We consider all four cases on delay parametersτ1andτ2.

Figure 16:Bifurcation diagram of the prey and predator species with respect to parameter τ2 when τ1=0.12 is fixed in its range of stability (0,τ10)

Figure 17:E?is locally asymptotically stable when τ2=0.42 is fixed in its stable range (0,τ′20)and τ1=0.1<τ10=0.1336

Case(I):Whenτ2=0 andτ1＞0,then we see that condition (H1) holds.Since the transversality condition is satisfied,therefore Hopf-bifurcation occurs in the system.To evaluate the critical value of delay parameter,takingi=0 in Eqs.(31) and (32),we obtain

ω1=0.3688,τ10=0.2889.

Thus,the positive equilibrium is locally asymptotically stable forτ1<τ10=0.2889,which is shown in Fig.8.Whenτ1=τ10,system undergoes a Hopf-bifurcation and periodic solution occurs aroundE?.The time series analysis and periodic solution have been shown in Fig.9.If we starts a trajectory from an initial point then it approaches to the periodic solution (Fig.9).This shows that the periodic solution is stable.In Fig.10,we made the bifurcation diagram for both the populations.The blue (red) curve represents the maximum (minimum) values of population at sufficiently large time.It is easy to see that Hopf-bifurcation occurs atτ1=τ10=0.2889.

Figure 18:E?is unstable when τ2=0.42 is fixed in its stable range (0,τ20) and τ1=0.2＞τ′10=0.1336.Time series solution of species and existence limit cycle

Case(II):Whenτ1=0 andτ2＞0.In this case,the transversality condition is satisfied,so the system will show Hopf-bifurcation at a critical value of delay parameterτ2.By some computation,we obtain

ω2=0.317,τ20=0.9618.

Therefore,according to our theoretical analysis,the system (2) is locally asymptotically stable forτ2<τ20.In Fig.11,we draw the time series of both the species forτ2=0.7<τ20=0.9618.From the figure,it can be seen that system is stable around the positive equilibriumE?.Atτ2=τ20,the system goes through a Hopf-bifurcation and forτ2＞τ20,system becomes unstable and limit cycle produces.This behavior is depicted in Fig.12.Again bifurcation diagram with respect to delayτ2for both the species is drawn in Fig.13,which helps us to understand the Hopfbifurcation phenomenon in the system.

Figure 19:Bifurcation diagram of the prey and predator species with respect to parameter τ1 when τ2=0.42 is fixed in its stable range (0,τ20)

Figure 20:Region of stability and instability for system (2) in τ1τ2-plane

Case (III):Whenτ1=0.12 (fixed in the interval (0,τ10)) andτ2as a parameter,then we observe that the condition (H2) holds true.Therefore according to Theorem 4.4 system (2)undergoes a Hopf-bifurcation.Eqs.(36) and (39) give us the values ofω0andτ′20as

ω0=0.445,τ′20=0.4731.

Thus the equilibrium pointE?is locally asymptotically stable forτ2<τ′20=0.4731 which is shown in Fig.14 and unstable forτ2＞τ′20(Fig.15).Whenτ2=τ′20,system undergoes a Hopfbifurcation aroundE?and periodic solution arises in the system.Bifurcation diagram is also presented in Fig.16 with respect toτ2for both the species whenτ1=0.12 (fixed).

Case (IV):Whenτ2=0.42 (fixed in the interval (0,τ20)) andτ1as a parameter,then our computer simulation yieldsω?=0.4491,τ′10=0.1336.

Forτ1=0.1 ∈(0,τ′10),the system is locally asymptotically stable (Fig.17).But forτ1=0.2＞τ′10,the system becomes unstable (Fig.18).Thus the model is stable forτ1<τ′10.Asτ1passes throughτ′10,it loses the stability and a Hopf-bifurcation occurs in the system.Fig.18 shows the existence of periodic solution (closed trajectory).The trajectory started from an initial point,approaches to the closed trajectory.This shows that the closed trajectory is stable.In Fig.19,we present the bifurcation diagram of both the species with respect toτ1whenτ2=0.42 (fixed).

As the system (2) shows Hopf-bifurcation with respect to both the delay parametersτ1andτ2.Therefore,we can bisect theτ1τ2—plane into two regions,which are separated by Hopf-bifurcation curve.

Region of stability (sky blue)S3={(τ1,τ2):system (2) is locally asymptotically stable},

Region of instability (white)S4={(τ1,τ2):system (2) is unstable}.

Both the regions are drown in Fig.20.

7 Conclusion

In this study,we have considered a habitat where two biological populations,prey populationxand predator populationyare surviving and interacting with each other.It is assumed that prey population follows logistic growth in the absence of predator and in the presence of predator,the interaction between them follows Holling type II functional response.We have shown the positivity,boundedness and persistence of the system,which implies that the proposed model is ecologically wellposed.We have defined a parameterA0(0≤A0≤1) which denotes the dependency of predators on supplied additional food.Our system has four kinds of equilibria,trivial equilibriumE0(0,0,0),axial equilibriumE1(K,0,0),two prey free equilibriaandunder condition (5) and unique positive equilibriumE?under conditions (9) and (10).Local and global stability of the positive equilibrium are shown under several conditions which are dependent upon the parameterA0.The parameterA0is crucial,so we have studied its effect via Hopf-bifurcation analysis which is also condescend by the numerical illustration.For a chosen set of parameters we calculated the threshold value of parameterA0,that isA0=0.482,where Hopfbifurcation occurs and system stabilizes.It is also observed that after stabilization of system if predators are more dependent on additional food then prey population increase whereas predators remain in their range.We also have studied the Hopf-bifurcation with respect to consumption rate of additional foodφ.Threshold value ofφis obtained asφ=0.02847.In Tab.3,we have shown the different number of positive equilibrium points by varying the parametric values,whena=0.105 andd=0.1 (other parameters are same as in (25)) then our system has two stable equilibrium together,therefore system shows the phenomenon of bistability,which is depicted in Fig.7.

Models with delay show comparatively more realistic dynamics than non delayed models.When a predator consumes a prey individual,then its effect does not come immediately,it takes some time i.e.,time lag for gestation.Again,predators also take some time to consume and digest the supplied additional food to them.Therefore,to make our model ecologically more realistic,we incorporated two delays;one for gestation delay and other for consuming and digesting the supplied additional food.

For the delayed model,we have analyzed Hopf-bifurcation via local stability taking delay as a bifurcation parameter.We investigated the Hopf-bifurcation phenomenon for all combinations of both delays.We obtained the sufficient conditions for the stability of the positive equilibrium point and existence of Hopf-bifurcation for Case(1):τ1＞0,τ2=0,Case(2):τ1=0,τ2＞0,Case(3):τ1is fixed in the interval (0,τ10) andτ2as a variable parameter,Case(4):τ2is fixed in the interval(0,τ20) andτ1as a variable parameter.Our system undergoes Hopf-bifurcation in the vicinity of the interior equilibrium point with respect to both the delay parameters when they cross their critical values.The qualitative properties of Hopf-bifurcation are studied by using the Normal form theory and the formulae given in Hassard et al.[34].

We have performed some numerical simulations to illustrate our theoretical results.For a biologically feasible set of parameters,the system is stable initially,then we introduce delay and system remains stable till its critical value.If we increase the delay parameter over the critical value,then system goes through Hopf-bifurcation and becomes unstable.Bifurcation diagrams(Figs.10,13,16,19) with respect to different delays depict the dynamical behavior of the system.

Our study is important to conserve the prey population through providing additional food to predators and to establish their balance.Here we have also shown the significance of delay parameters.We hope that this study will help to perceive the dynamics of an ecological system with additional food and two discrete delays.

Acknowledgement:The author (Ankit Kumar) acknowledges the Junior Research Fellowship received from University Grant Commission,New Delhi,India.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.