Global dynamics of an additional food provided predator–prey system with constant harvest in predators

Applied Mathematics and Computation 250 (2015) 193–211

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Global dynamics of an additional food provided predator–prey system with constant harvest in predators Moitri Sen a, P.D.N. Srinivasu b, Malay Banerjee a,⇑ a b

Department of Mathematics and Statistics, I.I.T. Kanpur, Kanpur 208016, India Department of Mathematics, Andhra University, Visakhapatnam 530003, India

a r t i c l e

i n f o

Keywords: Stable coexistence Oscillatory coexistence Extinction Local bifurcation Global bifurcation

a b s t r a c t The article aims to study the global dynamics associated with a predator prey system when the predator is provided with additional food and harvested at a constant rate. This study supplements the existing literature on the dynamics of additional food provided predator prey system by focusing on the consequences of harvesting the predators. It presents a comprehensive view on the entire range of bifurcations that take place in the considered system and highlights the dependence of the system dynamics on its vital parameters. This study provides important tools for investigations pertaining to controllability of the system which are essential from the real world applications perspective. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Dynamics of predator–prey systems when the predator is provided with additional food has become a topic of intense study due to its wide applications in biological control and species conservation. Some theoretical studies [2,20,28] conclude that provision of additional food to predators enhances the predator density which can influence the predator in a way that increases target predation resulting in lowering the density of target prey. This reduction in prey density due to presence of alternative food for predators is known to be apparent competition [19]. Whereas empirical studies indicate that provision of additional food to predators need not always increase target predation [16,20,24,37]. This apparent conflict between theory and observations led to in depth mathematical analysis of additional food provided predator–prey systems. Under the assumption that the additional food is uniformly distributed in the habitat and the number of encounters per predator with the additional food is proportional to the density of the additional food, it is observed that additional food serves as a biological control agent with which both prey as well as predators can be controlled. The nutritive value and the quantity of the additional food play vital role in the controllability of the predator–prey system [31,32,35]. By and large, it is observed that provision of additional food of high nutritive value could increase target predation and reduce the prey density and, additional food with low nutritive value could release the prey from predation pressure and decrease predator density. The recent works [27,31–34], present a comprehensive account on the dynamics associated with predator–prey systems in presence of additional food to predators. They also bring out controllability of the said predator–prey system using the quality (nutritive value) and quantity of the additional food as control parameters. These works not only supplement the understanding on pest management using biological means but also find their applicability in biological conservation and ⇑ Corresponding author. E-mail address: [email protected] (M. Banerjee). http://dx.doi.org/10.1016/j.amc.2014.10.085 0096-3003/Ó 2014 Elsevier Inc. All rights reserved.

194

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

resource management problems. Therefore, availability of the studies on effect of additional food on sustainability of the involved species leaves us with an important problem of assessing the considered system dynamics under harvesting which is the focus in the present article. Dynamics of predator–prey models with constant-yield harvesting (harvested biomass is independent of the population size) and constant-effort harvesting (a constant proportion of the population is harvested) are studied extensively by several authors [3,4,6–8,13,21,38]. Predator–prey models with the harvesting of either prey or predators or both the species have been studied as the outcomes of such analysis are helpful for the management of renewable resources [1,5]. Effect of harvesting on the dynamic behavior of the concerned model is not only interesting from theoretical point of view but also important for the management and sustainability of commercially/economically important resources [10–12,14,15,18,21]. Systematic analysis of harvested predator–prey models indicate the chance of extinction of one or both the species due to uncontrolled harvesting. Further, it determines the admissible level of harvesting to ensure the long term survival of the renewable resource. In this paper, we bring forward a systematic study on the dynamics of additional food provided predator–prey system with constant harvesting (also known as ‘constant-yield harvesting’) in predators. Dynamics of predator–prey system with constant harvesting in predators is presented in [39]. Essentially, this work presents complete bifurcation analysis for the considered model. Similar kind of dynamical analysis has been carried out for predator–prey systems in [26,29,36]. This analysis unfolds the rich dynamics exhibited by the predator–prey system when the predator is provided with additional food and also harvested at a constant rate. It is observed that constant harvesting brings in more complexity into the system dynamics when compared with the systems which are (i) free from harvesting; (ii) subjected to constant-effort harvesting [21,22]. This is due to inclusion of more number of equilibrium solutions and involvement of global bifurcation such as homoclinic bifurcation and codimension two local bifurcation namely Bogdanov–Takens bifurcation. Existence of such kind of bifurcations always indicate the over exploitation of resource and leading to the extinction scenario. The section-wise division of this paper is as follows. In the next section we describe the model representing dynamics of a predator–prey system when the predator species is provided with additional food and simultaneously harvested with constant harvesting. Section 3 presents analysis pertaining to the existence of equilibria and their local stability nature. Local bifurcation analysis and global dynamics of the system are presented in Sections 4 and 5 respectively. Section 6 presents an over view of the variety of dynamics of the considered system. Finally, Section 5 presents discussion and conclusions. 2. Basic model Let N and P represent biomass of prey and predator respectively. Let us assume that the predator is provided with additional food of biomass A which is distributed uniformly in the habitat. If h1 ðh2 Þ; e1 ðe2 Þ; n1 ðn2 Þ respectively represent the handling time of the predator per unit quantity of prey (additional food), ability of the predator to detect the prey (additional food) and the nutritive value of the prey (additional food) then the following system describes the predator–prey dynamics in presence of additional food to predators.

  dN N e1 NP  ¼ rN 1  ; dT K 1 þ e 1 h1 N þ e 2 h2 A

ð2:1aÞ

dP n1 e1 NP þ n2 e2 AP ¼  mP; dT 1 þ e1 h1 N þ e2 h2 A

ð2:1bÞ

where r and K respectively stand for intrinsic growth rate and carrying capacity of the prey population, and m is the death or starvation rate of the predator. Defining c ¼ h11 ; b ¼ n1 c; a ¼ e11h1 ; g ¼ nn21 ee21 and a ¼ nn12 hh21 , the system (2.1) takes the form

  dN N cNP  ¼ rN 1  ; dT K a þ agA þ N

ð2:2aÞ

dP bðN þ gAÞP ¼  mP dT a þ agA þ N

ð2:2bÞ

with all the involved parameters being positive. Here the parameters a; b and c stand for half saturation value of the predator in the absence of additional food, maximum birth rate of the predator due to consumption of the food perceived by the predator and maximum rate of capturing of prey by predator, respectively. Observe that the above system assumes Holling type II predator functional response towards its available food and that the number of encounters per predator with the additional food is proportional to additional food biomass. Further, this model reduces to the well known Rosenzweig–MacArthur model [5] in the absence of additional food to predators. The model (2.2) has been analyzed from biological control perspective [31] and studies pertaining to controllability of the system with respect to the parameters a and g have been presented in [32,34]. In the present study we wish to investigate consequences of harvesting the predators at a constant rate where in the predators are provided with additional food. Incorporating constant yield harvesting (q) into the model (2.2) we obtain the following system:

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

195

  dN N cNP ¼ rN 1   ; dT K a þ agA þ N

ð2:3aÞ

dP bðN þ gAÞP ¼  mP  q: dT a þ agA þ N

ð2:3bÞ

From physical considerations we need to note that the constant harvesting can be implemented in the system only if Pð0Þ > q and the harvesting activity ceases to exist beyond a specific time instant T ¼ s > 0 if PðsÞ ¼ 0 (due to non availability of predators in the system). Hence, without loss of generality we assume that Nð0Þ > 0; Pð0Þ > q for the considered system. Now we intend to study the dynamics of the above model and investigate the role of constant harvesting on the dynamics of considered predator–prey system. Before we proceed to analyze the model we non-dimensionalise it using the transformations x ¼ Na, y ¼ cP ; t ¼ rT and obtain the following equivalent form. ar

  dx x xy  ¼x 1 ¼ F 1 ðx; yÞ; dt 1 þ an þ x c

ð2:4aÞ

dy bðx þ nÞy ¼  dy  h ¼ F 2 ðx; yÞ; dt 1 þ an þ x

ð2:4bÞ

where c ¼ Ka ; b ¼ br, d ¼ mr ; n ¼ gaA and h ¼ arcq2 . Observe that the non-dimensional model (2.4) involves three exogenous parameters a; n and h. The first two describe the quality and quantity of the additional food [31,32,34] and the third is a measure for the harvest rate. Note that a is directly proportional to the handling time h2 of the additional food and inversely proportional to its nutritive value. Hence a is inversely related to the quality of the additional food and n represents the quantity of the additional food supplied to the predator. The parameters c; b and d represent the dimensionless versions of carrying capacity of x and maximum growth rate and death rate of y. All the parameters involved with the model system (2.4) are positive. 3. Existence of equilibria In this section we shall discuss the existence of equilibria depending upon the parametric restrictions for the system (2.4). h i bðxþnÞy The nullclines of the system under consideration are given by x 1  cx  1þaynþx ¼ 0 and 1þ anþx ¼ dy þ h. The system (2.4)   hð1þanÞ admits only one axial equilibrium point given by Ey ¼ 0; bndð1þanÞ which is a feasible equilibrium point for

D1  bn  dð1 þ anÞ > 0. The model (2.4) admits at most two interior equilibrium points and we denote them by Eh1 ðxh1 ; yh1 Þ and Eh2 ðxh2 ; yh2 Þ where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ½bn  dð1 þ anÞ  ðb  dÞc  ½bn  dð1 þ anÞ þ ðb  dÞc2  4hðb  dÞc h x1;2 ¼ ; ð3:5Þ 2ðb  dÞ   xh and yhi ¼ 1  ci ð1 þ an þ xhi Þ, ði ¼ 1; 2Þ. Conditions required for the feasible existence of interior equilibrium points are summarized in the following two propositions. For the sake of simplicity in the forthcoming discussion we define D2 and D3 as follows:

D2  bn  dð1 þ anÞ  h;

ð3:6Þ

D3  ½bn  dð1 þ anÞ þ cðb  dÞ2  4hcðb  dÞ:

ð3:7Þ

Proposition 3.1. Assume that b 6 d. (a) If D1 < 0 the system (2.4) has no interior equilibrium point. (b) If D1 > 0 then, (i) for D2 > 0 the system (2.4) has one feasible interior equilibrium point, (ii) for D2 < 0 the system (2.4) does not have any feasible interior equilibrium point. Proposition 3.2. Assume that b > d. (a) If D1 < 0 then, (i) for bn  dð1 þ anÞ þ ðb  dÞc > 0 the system (2.4) has no interior equilibrium point if D3 < 0 and two interior equilibrium points if D3 > 0, (ii) for bn  dð1 þ anÞ þ ðb  dÞc < 0 the system (2.4) has no interior equilibrium point.

196

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

(b) If D1 > 0 then, (i) for D2 > 0 the system (2.4) has unique interior equilibrium point, (ii) for D2 < 0 and bn  dð1 þ anÞ < ðb  dÞc the system (2.4) has no interior equilibrium point if D3 < 0 and two interior equilibrium points if D3 > 0, (iii) for D2 < 0 and bn  dð1 þ anÞ > ðb  dÞc the system (2.4) has no interior equilibrium point. The local asymptotic stability conditions of axial equilibrium point Ey and two interior equilibrium points Eh1 and Eh2 (whenever they exist) can be obtained from the nature of eigenvalues of the Jacobian matrix for the system (2.4) evaluated at the concerned equilibrium point. anÞh anÞ Two eigenvalues of the Jacobian matrix for the system (2.4) evaluated at Ey are bndð1þ and bndð1þ . The second eigenbndð1þanÞ 1þan value is positive whenever Ey is feasible. The local asymptotic stability of Ey can be summarized as follows: Lemma 3.3. Ey is an unstable node if D2 > 0 and is a saddle-point if D2 < 0. The trace and the determinant of the Jacobian matrix for the system evaluated at the interior equilibrium points are given by

Tr½J Eh  ¼

h i xh1;2 c  ð1 þ anÞ  2xh1;2

þ

bðxh1;2 þ nÞ

 d; 1 þ an þ xh1;2 2qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3 ½bn  dð1 þ anÞ þ ðb  dÞc2  4hðb  dÞc 5: Det½J Eh  ¼ xh1;2 4 1;2 cð1 þ an þ xh1;2 Þ 1;2

cð1 þ an þ

xh1;2 Þ

Here it is important to note that Det½J Eh  is always negative and the term Det½J Eh  is always positive. 2

1

The stability properties for the axial equilibrium point and interior equilibrium points are summarized below (we omit the proof as they involve straight forward algebraic calculations): Lemma 3.4. Eh1 is locally asymptotically stable if Tr½J Eh  < 0, is unstable when Tr½J Eh  > 0. 1

1

Lemma 3.5. Eh2 is a saddle-point whenever it exists. 4. Local bifurcation analysis In this section we present the local bifurcation results for the model under consideration. Here, conditions for transcritical, saddle-node, Hopf and Bogdanov–Takens bifurcations are derived. First three bifurcations are of co-dimension one and the fourth one is of co-dimension two [23]. The detailed calculations for the derivation of normal forms (near the bifurcation threshold) are presented in the appendix for the transcritical, saddle-node, Hopf and Bogdanov–Takens bifurcations. These normal forms are required to infer the nature of equilibrium points when the system parameters are at the bifurcation threshold. 4.1. Transcritical bifurcation   hð1þanÞ From the previous section it is clear that the model possesses only one axial equilibrium point given by Ey ¼ 0; bndð1þ anÞ . Also note that under D2 ¼ 0; Ey coincides with Eh1 if bn  dð1 þ anÞ < ðb  dÞc and Ey bn  dð1 þ anÞ > ðb  dÞc. Thus we state the following theorem.

coincides with Eh2

if

Lemma 4.1. If D1 > 0 the model undergoes a transcritical bifurcation when the system parameters satisfy the condition D2 ¼ 0. Proof. Here we need to verify the transversality condition for transcritical bifurcation using the Sotomayor’s theorem [25,30]. If D1 > 0, under the condition D2 ¼ 0 the equilibrium point is given by Ey  ð0; 1 þ anÞ. Now Det½J Ey  ¼ 0 and hence the Jacobian matrix J Ey has a zero eigenvalue. Let v and w be the two eigenvectors corresponding to the zero eigenvalue of the matrices J Ey and J TEy respectively and they are given by,

v¼



v1  ¼ v2

dð1þanÞbn bð1þannÞ

1

! ;

w¼



w1 w2

 ¼

  1 0

:

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

197

Thus we have,

wT F a ðEy ; aTC Þ ¼ 0; n½dð1 þ anÞ  bn – 0; bð1 þ anÞð1 þ an  nÞ     2 2 ðbn þ d þ dnaÞ2 2ðbn þ d þ dnaÞ bn ðn þ cÞb wT D2 FðEy ; aTC Þðv ; v Þ ¼  þ :  –0 d– &d – 2 2 ð1 þ naÞ ð1  na þ cÞ c ð1 þ naÞ b ð1 þ na  nÞ ð1 þ naÞbð1 þ na  nÞ

wT DF aðEy ; aTC Þv ¼

Hence the theorem. 4.2. Saddle-node bifurcation In Section 3 we have summarized the conditions for the existence of two interior equilibrium points (Eh1 and Eh2 ) based upon some restrictions satisfied by the system parameters. These equilibrium points are distinct if D3 > 0, they coincide if D3 ¼ 0 and they cease to exist if D3 < 0. In the case D3 ¼ 0 the components of the unique interior equilibrium point E ðx ; y Þ are given by

x ¼ 

½bn  dð1 þ anÞ  ðb  dÞc ; 2ðb  dÞ

ð4:8aÞ

  x ð1 þ an þ x Þ: y ¼ 1 

ð4:8bÞ

c

The emergence/annihilation of equilibria is due to the occurrence of saddle-node bifurcation for interior equilibrium points which takes place when the harvesting parameter (h) assumes the value

h ¼ hSN ¼

½bn  dð1 þ anÞ þ ðb  dÞc2 : 4ðb  dÞc

The following theorem presents occurrence of saddle-node bifurcation in dynamics of the system (2.4). Lemma 4.2. The system (2.4) undergoes saddle-node bifurcation when the system parameters satisfy the restriction D3 ¼ 0 along with the condition jD1 j < ðb  dÞc and D2 < 0 as given in the Proposition 3.2. Here the saddle-node bifurcation threshold with 2

anÞþðbdÞc respect to h as the bifurcation parameter is given by, h  hSN ¼ ½bndð1þ . 4ðbdÞc

Proof. Now we use Sotomayor’s theorem [25,30] to verify the transversality condition for the occurrence of saddle-node bifurcation at h ¼ hSN . It can be easily verified that Det½J E  ¼ 0 implying that J E admits a zero eigenvalue denoted by k1 . If V and W represent eigenvectors corresponding to the eigenvalue k1 for the matrices J E and J TE respectively, then they are given by

V¼



V1 V2

 ¼

1

!

bnbð1þanÞ cðbdÞ

;

W¼



W1 W2

 ¼

!

1 bndð1þanÞcðbdÞ  ½bndð1þ anÞþcðbdÞðbdÞ

:

Clearly, V and W satisfy the transversality conditions

W T F h ðE ; hSN Þ ¼ W 2 – 0; W T D2 FðE ; hSN Þ ¼ 

4x ðb  dÞ

c½bn  dð1 þ anÞ þ cðb  dÞ

–0

for the occurrence of saddle-node bifurcation at E when h ¼ hSN .

h

4.3. Hopf bifurcation In the previous section we have shown that Eh2 is always a saddle point whenever it exists and presented the condition required for local asymptotic stability of Eh1 . Furthermore, it can be easily concluded that the equilibrium point Eh1 may loose its stability through Hopf-bifurcation under certain parametric restrictions. Considering h as the bifurcation parameter one can detect the threshold magnitude h ¼ hH which satisfies Tr½J Eh jh¼hH ¼ 0 and Det½J Eh jh¼hH > 0. Observe that, 1

Tr½J Eh  ¼ 1

xh1 ½bn

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  bð1 þ anÞ þ ½bn  dð1 þ anÞ þ ðb  dÞc2  4hðb  dÞc

cðb  dÞð1 þ an þ xh1 Þ

1

þ

ðb  dÞxh1 þ bn  dð1 þ anÞ 1 þ an þ xh1

ð4:9Þ

198

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

can assume non positive value only if 1þnan < 1. Otherwise the interior equilibrium Eh1 is always unstable. The Hopf-bifurcation threshold is a positive root of Tr½J Eh  ¼ 0, say h ¼ hH > 0 and stability property of Eh1 changes when h passes through the crit1 ical magnitude h ¼ hH . Thus we summarize our findings in the following theorem. Lemma 4.3. Assume that 1 þ an > n and let system parameters satisfy the conditions for existence of interior equilibrium point given in Proposition 3.2. The interior equilibrium point Eh1 changes its stability through Hopf-bifurcation at the threshold h  hH such that Tr½J Eh jh¼hH ¼ 0. 1

Proof. The system may undergo Hopf-bifurcation only at Eh1 if 1þnan < 1. In order to ensure the change of stability through non-degenerate Hopf-bifurcation, we need to verify the transversality condition for the Hopf-bifurcation [9,17,25] which is given by 

  d TrðJ Eh Þ ¼ 1 dh h¼hH

  ð1 þ b  dÞc  ð1 þ anÞ  4xh1  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ð1 þ an þ xh1 Þ ½bn  dð1 þ anÞ þ ðb  dÞc2  4hðb  dÞc

– 0: h¼hH

The interior equilibrium point Eh1 loses it stability through non-degenerate Hopf-bifurcation when the parametric restriction Tr½J Eh jh¼hH ¼ 0 and the transversality condition mentioned above are satisfied simultaneously. Introducing the perturbations 1 x ¼ x1 þ xh1 jh¼hH ; y ¼ y1 þ yh1 jh¼hH in (2.4) and then expanding in Taylor series, we obtain

x_1 ¼ a10 x1 þ a01 y1 þ a20 x21 þ a11 x1 y1 þ a30 x31 þ a21 x21 y1 þ    ;

ð4:10aÞ

y_1 ¼ b10 x1 þ b01 y1 þ b20 x21 þ b11 x1 y1 þ b30 x31 þ b21 x21 y1 þ    ; where a10 ; a01 ; b10 ; b01 are the elements of the Jacobian matrix evaluated at the equilibrium point a10 þ b01 ¼ 0 and D ¼ a10 b01  a01 b10 > 0. The coefficients aij and bij are determined by

a20

" # 1 @ 2 F 1 ðx; yÞ ¼ 2 @x2

;

a11

ðEh1 ;hH Þ

a30

" # 1 @ 3 F 1 ðx; yÞ ¼ 6 @x3

b20

;

a21

;

b11

ðEh1 ;hH Þ

b30

" # 1 @ 3 F 2 ðx; yÞ ¼ 6 @x3

and h ¼ hH , hence

;

ðEh1 ;hH Þ

ðEh1 ;hH Þ

" # 1 @ 2 F 2 ðx; yÞ ¼ 2 @x2

" # @ 2 F 1 ðx; yÞ ¼ @x@y

ð4:10bÞ Eh1

" # 1 @ 3 F 1 ðx; yÞ ¼ 2 @x2 @y " # @ 2 F 2 ðx; yÞ ¼ @x@y

; ðEh1 ;hH Þ

;

ðEh1 ;hH Þ

;

ðEh1 ;hH Þ

b21

" # 1 @ 3 F 2 ðx; yÞ ¼ 2 @x2 @y

: ðEh1 ;hH Þ

The first Lyapunov number [25] to determine the stability of limit cycle arising through Hopf-bifurcation is given by,

   3p

a10 b10 a211 þ a10 a01 b211 þ a20 b11  2a10 a01 a220  a201 ð2a20 b20 þ b11 b20 Þ  a01 b10  2a210 a11 a20 2a01 D3=2 

 a210 þ a01 b10 ½3a01 a30 þ 2a10 ða21 þ b12 Þ þ ðb10 a12  a01 b21 Þ :

l1 ¼

If l1 < 0 the equilibrium point Eh1 destabilized through a supercritical Hopf-bifurcation, and if l1 > 0 then the Hopf-bifurcation is subcritical. h The expression for the first Lyapunov r number at the Hopf-bifurcation point is cumbersome and hence we are not providing its expression from the sake of brevity. We have calculated the Hopf-bifurcation curve in matcont and observed that the Hopf-bifurcation is always supercritical in nature as the first Lyapunov number is always negative. The absence of Bautin bifurcation point, on the Hopf-bifurcation curve, eliminates the chance of having unstable limit cycle. To support our claim for the stability of Hopf-bifurcating limit cycle we consider a numerical example. Let us choose parameter the values a ¼ 2; n ¼ 1, b ¼ 3; c ¼ 3; d ¼ 1:325 and the Hopf-bifurcation threshold for h is given by hH ¼ 0:3887 (correct up to 4th decimal place). For this choice of parameter values one can calculate the first Lyapunov number as l1 ¼ 1:64 < 0. Hence the Hopf-bifurcation is supercritical. 4.4. Bogdanov–Takens bifurcation Apart from the occurrence of codimension one bifurcations discussed so far, codimension two bifurcations such as Bogdanov–Takens bifurcation is also possible in the considered system (2.4). We shall illustrate the occurrence of this bifurca-

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

199

tion by considering h and d as bifurcation parameters. This co-dimension two bifurcation occurs at a point lying on the saddle-node bifurcation curve and hence we carry on the bifurcation analysis around the unique interior equilibrium point E , under the parametric restriction D3 ¼ 0. Let h ¼ hBT and d ¼ dBT be the threshold magnitude of bifurcation parameters such that DetðJ E Þjðd;hÞ¼ðdBT ;hBT Þ ¼ 0 and TrðJ E Þjðd;hÞ¼ðdBT ;hBT Þ ¼ 0. Here we would like to mention that, for the model under consideration, the expressions for dBT and hBT can not be obtained explicitly in terms of other parameters due to algebraic complexity of the expressions involved. Below we shall derive normal form of the Bogdanov–Takens bifurcation for the considered model (2.4), using the technique discussed in [39]. Lemma 4.4. Under the condition D3 ¼ 0, the unique interior equilibrium point E (see (4.8)) undergoes Bogdanov-Takens bifurcation at the threshold ðhBT ; dBT Þ if DetðJ E ÞjðdBT ;hBT Þ ¼ 0 and TrðJ E ÞjðdBT ;hBT Þ ¼ 0. Moreover if we give small perturbation to the parameters d and h around their BT-bifurcation threshold then the model (2.4) is topologically equivalent to the following model,

x_1 ¼ x2 ;

ð4:11aÞ

d2 x_2 ¼ l1 ðk1 ; k2 Þ þ l2 ðk1 ; k2 Þx2 þ x21 þ pffiffiffiffiffi x1 x2 þ Sðx1 ; x2 ; k1 ; k2 Þ: d1

ð4:11bÞ

where Sðx1 ; x2 ; k1 ; k2 Þ is a power series in x1 and x2 involving terms of the form xi1 xj2 satisfying i þ j P 3. Further, the three bifurcation curves through the BT-point in the k1 k2 plane are given by,

Saddle-node curve: SN ¼ ðk1 ; k2 Þ : l1 ðk1 ; k2 Þ ¼ 0 ,   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hopf bifurcation curve: H ¼ ðk1 ; k2 Þ : l2 ðk1 ; k2 Þ ¼ pd2ffiffiffiffi l1 ðk1 ; k2 Þ; l1 ðk1 ; k2 Þ < 0 , d1   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5d l ðk ;k Þ Homoclinic bifurcation curve: HL ¼ ðk1 ; k2 Þ : l2 ðk1 ; k2 Þ ¼ 2 p1ffiffiffiffi 1 2 ; l1 ðk1 ; k2 Þ < 0 . 7

d1

Proof. To obtain the analytical expressions for the saddle-node bifurcation curve, Hopf-bifurcation curve and homoclinic bifurcation curve passing through the Bogdanov–Takens bifurcation point, in the vicinity of BT point, we have to calculate the normal form. For the derivation of the normal form we perturb the parameters d and h around their BT-bifurcation threshold. Let us consider a perturbation to the parameters d and h around their BT-bifurcation threshold given by d ¼ dBT þ k1 and h ¼ hBT þ k2 where jk1 j; jk2 j 1. Substituting these perturbations in (2.4) we get

  dx x xy  ¼x 1  F 1 ðx; yÞ; dt 1 þ an þ x c

ð4:12aÞ

dy bðx þ nÞy ¼ þ ðd þ k1 Þy  h þ k2  F 2 ðx; y; k1 ; k2 Þ: dt 1 þ an þ x

ð4:12bÞ

Substituting n1 ; n2 for x  x ; y  y in (4.12) we obtain

n_1 ¼ an1 þ bn2 þ p11 n21 þ p12 n1 n2 þ P1 ðn1 ; n2 Þ;

ð4:13aÞ

n_2 ¼ k1 y þ k2 þ cn1 þ dn2 þ q11 n21 þ q12 n1 n2 þ P2 ðn1 ; n2 Þ;

ð4:13bÞ

where a; b; c; d are the elements of the Jacobian matrix evaluated at the equilibrium point E and d ¼ dBT , h ¼ hBT . The coefficients pij and qij are given by

p11 ¼

" # 1 @ 2 F 1 ðx; yÞ 2 @x2

;

p12 ¼

ðE ;dBT ;hBT Þ

p22 ¼

" # 1 @ 2 F 1 ðx; yÞ 2 @y2

q11 ¼

ðE ;dBT ;hBT Þ

q12

;

ðE ;dBT ;hBT Þ

;

" # @ 2 F 2 ðx; y; k1 ; k2 Þ ¼ @x@y

" # @ 2 F 1 ðx; yÞ @x@y

ðE ;dBT ;hBT Þ

" # 1 @ 2 F 2 ðx; y; k1 ; k2 Þ 2 @x2

;

ðE ;dBT ;hBT Þ

;

q22

" # 1 @ 2 F 2 ðx; y; k1 ; k2 Þ ¼ 2 @y2

;

ðE ;dBT ;hBT Þ

and P1 ; P2 are power series in n1 ; n2 with terms ni1 nj2 satisfying i þ j P 3. Using the affine transformation z1 ¼ n1 ; z2 ¼ an1 þ bn2 , the system (4.13) gets modified to

1 1 z_1 ¼ z2 þ n00 ðkÞ þ n10 ðkÞz1 þ n01 ðkÞz2 þ n20 ðkÞz21 þ n11 ðkÞz1 z2 þ n02 ðkÞz22 þ Q 1 ðz1 ; z2 Þ; 2 2

ð4:14aÞ

200

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

1 1 z_2 ¼ g00 ðkÞ þ g10 ðkÞz1 þ g01 ðkÞz2 þ g20 ðkÞz21 þ g11 ðkÞz1 z2 þ g02 ðkÞz22 þ Q 2 ðz1 ; z2 Þ; 2 2

ð4:14bÞ

where,

 ap  p n20 ðkÞ ¼ 2 p11  12 ; n11 ðkÞ ¼ 12 ; b b n02 ðkÞ ¼ 0; g00 ðkÞ ¼ bðk1 y þ k2 Þ; g10 ðkÞ ¼ ak1 ; g01 ðkÞ ¼ k1 ;   ap  a2 p 12 g20 ðkÞ ¼ 2 ap11  12 þ bq11  aq12 ; g11 ðkÞ ¼ þ q12 ; g02 ðkÞ ¼ 0; b b n00 ðkÞ ¼ 0;

n10 ðkÞ ¼ 0;

n01 ðkÞ ¼ 0;

and Q 1 ; Q 2 are power series in z1 and z2 with terms zi1 zj2 satisfying i þ j P 3. Now we use C 1 change of coordinates in a small neighborhood of ð0; 0Þ in (4.14),

v 1 ¼ z1 

p12 2 z ; 2b 1

v 2 ¼ z2 þ

 ap  p11  12 z21 ; b

to obtain

v_1 ¼ v 2 þ R1 ðv 1 ; v 2 Þ;

ð4:15aÞ

v_2 ¼ bðk1 y þ k2 Þ  ak1 v 1 þ k1 v 2 þ



ap11 

  a2 p12 ak1 ap  þ bq11  aq12  p12  k1 p11  12 v 21 b 2b b

  ap þ 2p11  12 þ q12 v 1 v 2 þ R2 ðv 1 ; v 2 Þ; b

ð4:15bÞ

where R1 and R2 are again two power series involving terms like v i1 v j2 satisfying i þ j P 3. Finally we use another C 1 change of coordinates in a small neighborhood of ð0; 0Þ in (4.15)

u1 ¼ v 1 ;

u2 ¼ v 2 þ R2 ðv 1 ; v 2 Þ

and obtain

u_1 ¼ u2 ;

ð4:16aÞ

  ap u_2 ¼ bðk1 y þ k2 Þ  ak1 u1 þ k1 u2 þ 2p11  12 þ q12 u1 u2 þ G1 ðu1 Þ b    a2 p12 ak1 ap  þ bq11  aq12  p12  k1 p11  12 u21 þ u2 G2 ðu1 Þ þ u22 G3 ðu1 ; u2 Þ; þ ap11  b 2b b

ð4:16bÞ

where G1 and G2 are power series involving the terms uk11 ; uk12 respectively with k1 P 3, k2 P 2. G3 is another power series which involves terms like ui1 uj2 such that i þ j P 1. Condition for Bogdanov–Takens bifurcation of co-dimension 2 is given by



a b c d

 – h2 2 ; ðk1 ¼0;k2 ¼0Þ

and the transversality conditions are

  ap 2p11  12 þ q12 – 0; b    a2 p12 ak1 ap  ap11  þ bq11  aq12  – 0: p12  k1 p11  12 b 2b b ðk1 ¼0;k2 ¼0Þ Thus, we have

u_1 ¼ u2 ;

ð4:17aÞ

u_2 ¼ Pðu1 ; kÞ þ u2 M1 ðu1 ; kÞ þ u22 N1 ðu1 ; u2 Þ:

ð4:17bÞ

where P; M 1 ; N 1 2 C 1 and

   a2 p12 ak1 ap  Pðu1 ; kÞ ¼ bðk1 y þ k2 Þ  ak1 u1 þ G1 ðu1 Þ þ ap11  þ bq11  aq12  p12  k1 p11  12 u21 b 2b b   ap12 M 1 ðu1 ; kÞ ¼ k1 þ 2p11  þ q12 u1 þ G2 ðu1 Þ b N1 ðu1 ; u2 Þ ¼ G3 ðu1 ; u2 Þ: 0Þ So, Pð0; k0 Þ ¼ @Pð0;k ¼ 0; @ @u1

2

Pð0;k0 Þ @u21

0Þ ¼ 2d1 – 0; M 1 ð0; k0 Þ ¼ 0; @M1@uð0;k ¼ d2 – 0. Here k0  fk1 ¼ 0; k2 ¼ 0g. 1

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

201

Applying the Malgrange Preparation Theorem on Pðu1 ; kÞ we have

Pðu1 ; kÞ ¼ ð/1 ðkÞ þ /2 ðkÞu1 þ u21 ÞBðu1 ; kÞ; u2 ffi;s ¼ where, /1 ; /2 ; B 2 C 1 and Bð0; k0 Þ ¼ d1 ; /i ðk0 Þ ¼ 0; i ¼ 1; 2. Let y1 ¼ u1 , y2 ¼ pffiffiffiffiffiffiffiffiffiffi Bðu1 ;kÞ (4.17) becomes

R t pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Bðu1 ðsÞ; kÞds. Then the system 0

y_1 ¼ y2 ;

ð4:18aÞ

y_2 ¼ /1 ðkÞ þ /2 ðkÞy1 þ y21 þ y2 ðwðkÞ þ Mðy1 ; kÞÞ þ y22 Nðy1 ; y2 ; kÞ

ð4:18bÞ

M 1 ð0;kÞ M 1 ðu1 ;kÞ ðu1 ;u2 ;kÞ ffiffiffiffiffiffiffiffiffi ; Mðu1 ; kÞ ¼ p ffiffiffiffiffiffiffiffiffiffiffi ; Nðu1 ; u2 ; kÞ ¼ N1Bðu where wðkÞ ¼ p . 1 ;kÞ Bð0;kÞ

Bðu1 ;kÞ

Now, considering the transformation

x1 ¼ y 1 þ

/2 ðkÞ ; 2

x2 ¼ y2

we obtain

x_1 ¼ x2 ;

ð4:19aÞ

d2 x_2 ¼ l1 ðk1 ; k2 Þ þ l2 ðk1 ; k2 Þx2 þ x21 þ pffiffiffiffiffi x1 x2 þ Sðx1 ; x2 ; kÞ: d1

ð4:19bÞ

where

l1 ðk1 ; k2 Þ ¼ /1 ðkÞ 

/22 ðkÞ ; 4

d

l2 ðk1 ; k2 Þ ¼ wðkÞ  p2ffiffiffiffiffi /2 ðkÞ; 2 d1

and Sðx1 ; x2 ; kÞ is a power series in x1 and x2 involving terms of the form xi1 xj2 satisfying i þ j P 3. h i @ðl ;l Þ The system undergoes a Bogdanov–Takens bifurcation if the matrix @ðk11 ;k22Þ is nonsingular. When the model ðk1 ¼0;k2 ¼0Þ

undergoes Bogdanov–Takens bifurcation at k1 ¼ 0 ¼ k2 , we find three bifurcation curves in k1 k2  plane and through the BTpoint and they are given by, Saddle-node curve: SN ¼ fðk1 ; k2 Þ : l1 ðk1 ; k2 Þ ¼ 0g,   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Hopf bifurcation curve: H ¼ ðk1 ; k2 Þ : l2 ðk1 ; k2 Þ ¼ pd2ffiffiffiffi l1 ðk1 ; k2 Þ; l1 ðk1 ; k2 Þ < 0 , d1   pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5d2ffiffiffiffi l1 ðk1 ; k2 Þ; l1 ðk1 ; k2 Þ < 0 . h Homoclinic bifurcation curve: HL ¼ ðk1 ; k2 Þ : l2 ðk1 ; k2 Þ ¼ p 7

d1

Example 4.5. Now we illustrate locations of these bifurcation curves and dynamics of the system at BT bifurcation point with specific choices for the dimensionless parameters numerically. We choose parameter values a ¼ 2; b ¼ 3; c ¼ 3; n ¼ 1, then threshold magnitudes are given by dBT ¼ 1:3043793; hBT ¼ 0:8561257094853. Substituting these parameter values, we get from (4.14)

n00 ðkÞ ¼ 0;

n10 ðkÞ ¼ 0;

n01 ðkÞ ¼ 0;

n20 ðkÞ ¼ :2473147070;

n11 ðkÞ ¼ 0:3555309510;

n02 ðkÞ ¼ 0; g00 ðkÞ ¼ :7258323869k1  :3709720606k2 ; g10 ðkÞ ¼ :4375648878k1 ; g01 ðkÞ ¼ k1 ; g20 ðkÞ ¼ 0:4193519533; g11 ðkÞ ¼ 0:1082162383; g02 ðkÞ ¼ 0; Thus we have,

n20 ð0Þ þ g11 ð0Þ ¼ :13909847 – 0;

and g20 ð0Þ ¼ 0:4193519533 – 0;

and hence the Bogdanov–Takens bifurcation is non-degenerate [23]. Finally, we can compute the quantity r determining the bifurcation structure of the BT point [23],

r ¼ signðg20 ð0Þðn20 ð0Þ þ g11 ð0ÞÞ ¼ signð0:05833121455Þ < 0: Hence we can conclude that the system (2.4) undergoes a supercritical Bogdanov–Takens bifurcation. Now from (4.16) we have,

u_1 ¼ u2 ;

ð4:20aÞ

u_2 ¼ 0:7258323869k1  0:3709720606k2 þ 0:4375648878k1 u1 þ ð0:2096759767 þ 0:2014412838k1 Þu21 þ G1 ðu1 Þ þ k1 u2  0:1390984687u1 u2 þ u2 G2 ðu1 Þ þ u22 G3 ðu1 ; u2 Þ:

ð4:20bÞ

202

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

Pðu1 ; kÞ is given by

Pðu1 ; kÞ ¼ 0:7258323869k1  0:3709720606k2 þ 0:4375648878k1 u1 þ ð0:2096759767 þ 0:2014412838k1 Þu21 þ G1 ðu1 Þ; and the applying the Malgrange Preparation Theorem on Pðu1 ; kÞ we obtain

 Pðu1 ; kÞ ¼ /1 ðkÞ þ /2 ðkÞu1 þ u21 Bðu1 ; kÞ;

where

0:7258323869k1  0:3709720606k2 ; 0:2096759767 þ 0:2014412838k1

/1 ðkÞ ¼

/2 ðkÞ ¼

0:4375648878k1 ; 0:2096759767 þ 0:2014412838k1

Bð0; kÞ ¼ 0:2096759767 þ 0:2014412838k1 : u2 ffi;s ¼ Using the transformations y1 ¼ u1 , y2 ¼ pffiffiffiffiffiffiffiffiffiffi Bðu1 ;kÞ

R t pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Bðu1 ðsÞ; kÞds in (4.20), we find 0

y_1 ¼ y2 ;

ð4:21aÞ

k1 0:1390984687 y_2 ¼ /1 ðkÞ þ /2 ðkÞy1 þ y21 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi y1 y2 þ G4 ðy1 ; y2 ; kÞ; Bðy1 ; kÞ Bðy1 ; kÞ

ð4:21bÞ

where G4 is a power series with terms yi1 yj2 satisfying i þ j P 3 and j P 2. Using the affine transformation x1 ¼ y1 þ /22ðkÞ ; x2 ¼ y2 in (4.21) we obtain

x_1 ¼ x2 ;

ð4:22aÞ

d2 x_2 ¼ l1 ðk1 ; k2 Þ þ l2 ðk1 ; k2 Þx2 þ x21 þ pffiffiffiffiffi x1 x2 þ G5 ðx1 ; x2 ; kÞ: d1

ð4:22bÞ

where

l1 ðk1 ; k2 Þ ¼ /1 ðkÞ 

/22 ðkÞ ; 4

d

l2 ðk1 ; k2 Þ ¼ wðkÞ  p2ffiffiffiffiffi /2 ðkÞ;

ð4:23Þ

2 d1

k1 ffi ¼ wðkÞ. Also Bð0; 0Þ ¼ d1 ¼ 0:2096759767; d2 ¼ 0:1390984687 and G5 is a power series in x1 and x2 involving and pffiffiffiffiffiffiffiffi Bð0;kÞ

terms of the form xi1 xj2 satisfying i þ j P 3 and j P 2. Since

  @ðl1 ; l2 Þ   ¼ 4:424626862 – 0;  @ðk ; k Þ  1 2 ðk1 ¼0;k2 ¼0Þ the rank of the matrix

h

i

@ðl1 ;l2 Þ @ðk1 ;k2 Þ

ðk1 ¼0;k2 ¼0Þ

is 2 and hence the parametric transformation (4.23) is nonsingular. Hence the system

of Eqs. (4.22) can be written as follows 4

3.5

2.5

λ2 →

predator →

3

2

H SN HL

1.5

1

0.5

0 0

0

0.5

1

1.5

2

prey →

2.5

3

3.5

4

λ → 1

Fig. 1. (a) phase portrait when parameter values satisfy the conditions for Bogdanov–Takens bifurcation; (b) the bifurcation diagram in a small neighborhood of the Bogdanov–Takens bifurcation point in k1 k2 -parameter space, the curves labelled with SN; H and HL denote the Saddle-Node, Hopf and Homoclinic bifurcation curves respectively.

203

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

x_1 ¼ x2 ;

ð4:24aÞ

x_2 ¼ l1 ðk1 ; k2 Þ þ l2 ðk1 ; k2 Þx2 þ x21  0:303772189616647x1 x2 þ G5 ðx1 ; x2 ; kÞ:

ð4:24bÞ

Hence the model (2.4) undergoes Bogdanov–Takens bifurcation when k1 ¼ 0 ¼ k2 for chosen parameter values and the bifurcation curves in ðk1 ; k2 Þ plane and passing through the BT-point (within a small neighborhood of ðk1 ; k2 Þ ¼ ð0; 0Þ) are given pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi by, l1 ðk1 ; k2 Þ ¼ 0, (saddle-node bifurcation curve); l2 ðk1 ; k2 Þ ¼ 0:303772189616647 l1 ðk1 ; k2 Þ, (Hopf bifurcation curve); pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi and l2 ðk1 ; k2 Þ ¼ 0:216980135440462 l1 ðk1 ; k2 Þ, (Homoclinic bifurcation curve) where l1 ðk1 ; k2 Þ < 0 for the last two bifurcation curves. These bifurcation curves are shown in Fig. 1 along with the phase portrait.

5. Global dynamics In this section we intend to throw some light on the global behavior of the considered system. In the previous two sections we obtained the local asymptotic stability results for the equilibrium points and derived conditions for the occurrence of local bifurcations. Based upon those analytical results now we can prepare various bifurcation diagrams which enable us to understand the complete variety of dynamics exhibited by the model under consideration for all possible parameter values. We have discussed in Section 4 that the interior equilibrium point Eh2 is a saddle point whenever it is feasible. The stability property of Eh1 remains unaltered under the parametric restriction (H1): unstable node if (H1) holds. The local stability behavior of

Eh1

n 1þan

> 1; the equilibrium point Eh1 is an

can change under the parametric restriction (H2):

n 1þan

< 1.

In this case Eh1 can loose its stability through Hopf-bifurcation and as a consequence we find existence of stable limit cycle encircling Eh1 . The size of the limit cycle generated through the Hopf-bifurcation changes with variation in model parameters and it disappears through homoclinic bifurcation, which is a global bifurcation. Under the parametric restriction (H2) Bogdanov–Takens bifurcation takes place at the point of intersection of saddle-node and Hopf-bifurcation curves. Here we discuss two sets of qualitatively different bifurcation diagrams under the parametric restrictions (H1) (see Fig. 2 and (H2) (see Fig. 3). All the bifurcation diagrams are prepared in dh-parameter space. In each of the bifurcation diagrams the vertical green line represents the line d ¼ 1þbnan, that is the curve D1 ¼ 0. The transcritical and saddle-node bifurcation curves, corresponding to the equations D2 ¼ 0 and D3 ¼ 0, are marked with red and blue colors respectively. First we describe two bifurcation diagrams presented in Fig. 2, the system parameters a and n satisfy the restriction (H1), which attributes unstable nodal nature to Eh1 and hence refrains the system from under going Hopf bifurcation, thereby avoiding the occurrence of homoclinic and BT bifurcations. In the sub-figures of Fig. 2, the yellow line represents the straight line d ¼ b. The curves D1 ¼ 0 and D2 ¼ 0 always intersect at the point ð1þbnan ; 0Þ in dh-parameter space and hence lies on the horizontal axis. Relative magnitudes of c and n are responsible for the generation of two bifurcation diagrams Fig. 2(a) and (b). For c > n, a portion of the D3 ¼ 0 (saddle-node bifurcation) curve is present within the positive quadrant of dh-plane and   bðcnÞ bc½nð1þanÞ (see (Fig. 2(a))). The bifurcation curves D1 ¼ 0, intersect the transcritical bifurcation curve at cð1þ anÞ ; cð1þanÞ

D2 ¼ 0; D3 ¼ 0 and the line d ¼ b divide, the parametric domains presented in Fig. 2(a) and (b) into six and five regions respectively. In the region F the model possesses no equilibrium point. The dynamics in the domains E and G are qualitatively similar. In these regions only the axial equilibrium point Ey exists and it is a saddle point. Again the dynamics in domains D

γ>ξ

γ<ξ

E

E G G

h→

h→

C

F

F D

D

H

H

δ→

(a)

δ→

(b)

Fig. 2. Bifurcation diagram for different values of c under the condition 1 þ an < n. For further details see Table 1.

204

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

E

A

BT

E

h→

h→

A C

3

C1

D

BT

D

C3

C2

C1

C

B

2

B 0

0

δ→

δ→

E BT

E

A C1

A C h→

h→

3

B2

C

BT

3

D

D

C1

C2

B1

B2 0

B1

0

C2

δ→

δ→

E

A

E h→

h→

A

D

C

δ→

B2

BT

D

B2 B3

0

C

B1

B

3

BT

B1

0 0

δ→

E E C1 C2

C

A

A h →

h→

BT

3

C D

B

B

2

3

D

B2

B3

BT

B1

B1 δ→

0 0

δ→

Fig. 3. Bifurcation diagram for different magnitudes of c under the condition 1 þ an > n. For further details see Table 1.

205

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

and H are qualitatively similar. In these regions one interior equilibrium point Eh2 and the axial equilibrium point Ey exist and they are saddle point and unstable node respectively. In the region C (appeared in Fig. 2(a) only) we find two interior equilibrium points Eh1 and Eh2 along with the axial equilibrium point Ey . Eh1 is an unstable node, Eh2 and Ey are saddle points. The model admits at most one interior equilibrium point if c < n and at most two interior equilibrium points if c > n but none of them are stable whenever the parametric restriction (H1) holds. Although the dynamics of the model under consideration is same in the domains E; G and D; H, the differentiation is required to ensure that the point of intersection of the saddle-node and transcritical bifurcation curves should lie within the domain d < b and occurrence of one interior equilibrium point under the parametric restriction (H1). Now we consider the possible bifurcation diagrams under the parametric restrictions (H2). These diagrams are presented in Fig. 3. In all the sub-figures of Fig. 3, the curve D3 ¼ 0 (saddle-node bifurcation curve) divides the domain D1 < 0 into two parts and we label them as A and B respectively. Similarly the curves corresponding to D3 ¼ 0 (saddle-node bifurcation curve) and D2 ¼ 0 (transcritical bifurcation curve) divide the domain D1 > 0 into three parts and we label them as D; C   bðcnÞ bc½nð1þanÞ and E. Also D2 ¼ 0 and D3 ¼ 0 touch each other at cð1þ anÞ ; cð1þanÞ , we mark this point by a black filled circle. This point will be in the interior of the positive quadrant of the hd plane only when c < n and it occurs only in the region D1 > 0. Under the parametric restriction (H2), there will be no feasible equilibrium point for the model in the region A. Only axial equilibrium point exists in the region E which is a saddle point. The region D signifies the region where Eh2 and Ey exist where Eh2 is a saddle point and Ey is an unstable node. Dynamics of the system under consideration is relatively rich in the region B [ C, due to presence of three equilibria and occurrence of various local and global bifurcations. Now we shall concentrate on the system dynamics when the parameters belong to B [ C. In order to visualize the Hopf, homoclinic bifurcation curves and BT bifurcation point we have plotted the curves D1 ¼ 0; D2 ¼ 0; D3 ¼ 0, Tr½J Eh  ¼ 0 and the homoclinic bifurcation curve in the ðdhÞ-parameter space and these 1

h¼hH

curves are always confined within the region B [ C. The Hopf and homoclinic bifurcation curves are presented in cyan and magenta colors and the Bogdanov–Takens bifurcation point by a red star. The Bogdanov–Takens bifurcation point is the point of intersection of the saddle-node and Hopf bifurcation curve. It is important to observe that the curve D1 ¼ 0 forms the boundary between the regions B and C. Also the system admits only two interior equilibria in B where as it admits two interior equilibria and an axial equilibrium Ey within C. Hence, the global dynamics of the system in C significantly differ from the dynamics in B due to the presence of the axial equilibrium, which is of saddle nature. The bifurcation diagrams presented in the sub-figures of Fig. 3 differ from one another due to the relative positions of the Hopf and homoclinic bifurcation curves with respect to the line D1 ¼ 0. Solving the equations TrðJ E Þ ¼ DetðJ E Þ ¼ 0 one can find the coordinates of BT bifurcation point in the ðdhÞ plane and the magnitudes of the components of BT point determine its

Table 1 Summary of the bifurcation diagrams under various possible parametric restrictions. 1 þ an  n

1þan b  1þ ann

Restriction on n

Restriction on c

Bifurcation diagram

>0

>0

0 < n < c2 < c1 < 1 þ an

0
Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig. Fig.

0
Fig. 2(a) Fig. 2(b)

0 < c2 < n < c1 < 1 þ an

0 < c2 < c1 < n < 1 þ an

<0

0 < n < c2 < 1 þ an < c1

0 < c2 < n < 1 þ an < c1

<0

NA

NA

3(a) 3(b) 3(d) 3(f) 3(h) 3(a) 3(c) 3(d) 3(f) 3(h) 3(a) 3(c) 3(e) 3(f) 3(h) 3(a) 3(b) 3(d) 3(g) 3(h) 3(a) 3(c) 3(d) 3(g) 3(h)

206

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211 2

ð1þanÞ position relative to the line D1 ¼ 0. We observe that the BT bifurcation point is located on the line D1 ¼ 0 if c ¼ bð1þ annÞ  c1

(say) and it lies in the region D1 > ð<Þ0 if c < ð>Þc1 . Analyzing TrðJ Eh Þ ¼ 0 we observe that when c < c1 the Hopf bifurcation 1

2

ð1þanÞ curve intersects the D1 ¼ 0 line if c > ð1þbÞð1þ anÞbn  c2 (say). Again the Hopf bifurcation curve intersects the horizontal axis

an the conh ¼ 0 on the right side of the point ð1þbnan ; 0Þ if c > 1 þ an. Otherwise it starts from ð1þbnan ; 0Þ. Also note that if b > 1þ1þann

dition c1 < 1 þ an is automatically satisfied. This suggests us to consider two situations c1 < 1 þ an and c1 > 1 þ an to prepare various bifurcation diagrams. an Starting with b > 1þ1þann (i.e. c1 < 1 þ an), we find four possible division for the domain B [ C by the Hopf and homoclinic bifurcation curves. If (i) 0 < c < c2 the Hopf and homoclinic bifurcation curves are confined within the domain C only and the BT bifurcation point appear on the boundary between the regions C and E. Here the Hopf and homoclinic bifurcation curves divide the region C into three sub-regions C 1 ; C 2 and C 3 . If (ii) c2 < c < c1 the Hopf bifurcation curve initiates from ð1þbnan ; 0Þ and intersects D1 ¼ 0 line, the BT bifurcation point appear on the boundary between the regions C and E. In this situation the regions B and C are divided into the sub-regions B1 ; B2 and C 1 ; C 2 and C 3 , respectively. If (iii) c1 < c < 1 þ an or (iv) c > 1 þ an the BT bifurcation point appears on the boundary between the regions B and A. Under these parametric restrictions, the Hopf bifurcation curve does not intersect the line D1 ¼ 0. So, the Hopf and homoclinic bifurcation curves stays in the domain B and divide it into three sub-regions B1 ; B2 and B3 . The only visible difference is that the Hopf bifurcation curve starts from the point ð1þbnan ; 0Þ when c1 < c < 1 þ an and it starts from a point on the d axis lying in the region D1 < 0 if c > 1 þ an. Hopf     cð1þanÞ ; 0 when c > 1 þ an and from the point 1þbnan ; 0 otherwise. In these cases the bifurcation curve starts from b½2nþ 1þanþc phase portrait for parameter values lying within the region C is topologically equivalent to the parameter values lying in the region C 3 .

Region E or G

y→

y →

Region A or F

x→

x→

y→

Region D or H

x→

Fig. 4. Phase portraits for parameter values lying within the regions A or F; E or G and D or H.

207

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

an For b < 1þ1þann , (i.e. c1 > 1 þ an), we have another four possible situations depending upon the magnitude of c relative to the thresholds c1 and C2 : (v) 0 < c < c2 , (vi) c2 < c < 1 þ an, (vii) 1 þ an < c < c1 and (viii) c1 < c. In case of (v), (vi) and (viii) the position of the Hopf bifurcation curve and the Homoclinic bifurcation curve will be same as that is in the case (i), (ii) and (iv) respectively. For 1 þ an < c < c1 , the BT point belongs to the boundary between C and E and Hopf bifurcation curve intersects the line D1 ¼ 0. As a result each of the domains B and C are divided into three sub-regions B1 to B3 and C 1 to C 3 . The phase portraits for the parameter values lying in different domains are presented in Figs. 5 and 6. Region B1 is characterized by the existence of one saddle point Eh2 and one stable point Eh1 . Here Eh1 is locally asymptotically stable whose basin of attraction is bounded by the stable manifolds of Eh2 . In the region B2 ; Eh1 losses its stability through Hopf-bifurcation and one stable limit cycle gets generated whose basin of attraction is bounded by the stable manifolds of the saddle point Eh2 . As the parameter values enter the region B3 the stable limit cycle disappears through Homoclinic bifurcation and the region is left with one unstable focus Eh1 and one saddle interior equilibrium point. In the case D1 > 0 when the parameter values lie in the region C along with Eh1 and Eh1 the axial equilibrium point Ey also exists and it is an unstable point. Now in the subdivision C 1 ; Eh2 is a saddle node and Eh1 is locally stable with the stable manifolds of Eh2 being the boundaries for its basin of attraction. In the region C 2 ; Eh1 losses its stability and a stable limit cycle is formed whose basin of attraction is bounded by the stable manifolds of the saddle point Eh2 . Lastly in the region C 3 , limit cycle around Eh1 disappears through homoclinic bifurcation. In all the figures the unstable equilibrium points are marked with black circles, stable equilibrium points with black filled up circles. Stable and unstable manifolds of the equilibrium points are presented with magenta color curves and red closed curve stands for stable limit cycles. All possible parametric domains described so far and related bifurcation diagrams under various parametric restrictions are summarized in Table 1.

6. Over view of the system dynamics In this section we present the influence of provision of additional food to predator and constant harvesting on the dynamics of predator–prey system. Here we recall that the parameters a and n are related to quality and quantity of additional food respectively. From the definitions of these parameters we observe that a is inversely related to the quality and n is directly related to the quantity of additional food and hence we take them as representatives for quality and quantity of the additional food provided to the predator per unit time. Now we define the term 1n þ a as characteristic of the additional food. Given the values of a and n for the additional food, the dynamics of the system can be characterized depending on the relative magRegion B1

Region B2

5.5

4

5

3.5

4

3

3.5

2.5

3

y→

y→

4.5

2.5

2 1.5

2 1.5

1

1 0.5

0.5 0 0

0.5

1

1.5

2 x→

2.5

3

3.5

0 0

4

0.5

1

1.5

x→

2

2.5

3

3.5

Region C3/C

Region B

3 6

8 7

5

6 4 y→

y→

5 4 3

3

2

2 1 1 0 0

0.5

1

1.5

x→

2

2.5

3

3.5

0 0

0.5

1

1.5 x→

2

2.5

3

Fig. 5. Phase portraits corresponding to the domains B1 ; B2 ; B3 and C of the bifurcation diagram presented in Fig. 3(c). Here the values of the parameters are d ¼ 1:4; h ¼ :4 (B1 ), d ¼ 1:325; h ¼ :6 (B2 ), d ¼ 1:1; h ¼ :5 (B3 ), d ¼ :5; h ¼ 2:3 (C). The other parameters a ¼ 2; b ¼ 3; c ¼ 3; n ¼ 1 are fixed for each domain. Phase portraits for parameter values lying in the domain B of Fig. 3(c) and in the domain C of Fig. 3(c) and (d) will be topologically equivalent to phase portraits in the domains B1 and C 3 respectively. Stable equilibrium point is marked with filled black circle, unstable equilibrium points are marked with black circles and stable limit cycle is indicated in red color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

208

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

nitudes of the additional food characteristic and bd with respect to 1 as described bellow. It is pertinent to note that the magnitudes of the characteristic of additional food and bd decide the sign of bn  dð1 þ anÞ. The parameters b and d respectively represent maximum growth rate of the predator in the absence of additional food and their natural death rate. First, let us consider the system with bd < 1, that is additional food independent maximum growth rate of predator is less than the death rate. Under this assumption, the system without provision of additional food does not exhibit coexistence due to domination of the death rate over the maximum growth rate of the predator. As a result the predator gets extinct in a finite time under constant harvest (however small is the harvest rate). With the provision of additional food to the predator three distinct cases arise. They are given by (i) bd < a þ 1n < 1, (ii) bd < 1 < a þ 1n and (iii) a þ 1n < bd < 1. In the first two cases the system does not posses any equilibrium point and the predator goes to extinction even if the predator is provided with additional food (Fig. 4(a)). On the other hand, if the parametric restriction (iii) is satisfied then from Table 1 and Fig. 2 it is clear that irrespective of the environmental carrying capacity of prey (c), the system admit at most one interior equilibrium Eh2 which is a saddle point. In this case if the constant harvest rate, ðhÞ, is smaller than the value of bn  dð1 þ anÞ, the system admits an unstable axial equilibrium point Ey on the predator axis and one saddle interior equilibrium point Eh2 . The stable manifolds of Eh2 divides the first quadrant of the phase space into two disjoint regions with one of its branches connecting the axial equilibrium point Ey to the interior equilibrium point Eh2 which is heteroclinic connection between them. The eventual state of the system depends on its initial state relative to the stable manifold of the interior equilibrium. Although, it is difficult to maintain the state of the system at the interior equilibrium (due to its instability nature) still it is possible to have coexistence of predator and prey by ensuring that the system initiates on the stable manifold of the interior equilibrium. Again if the system initiates above the stable manifold then the predator survive at all future time on the additional food only whereas the prey goes to extinction after a finite time. If the initial value is situated below the stable manifold then the predator gets extinct in a finite time (Fig. 4(c)). As the harvesting rate (h) gets increased the system exhibits transcritical bifurcation for h ¼ bn  dð1 þ anÞ, at which the interior equilibrium collides with the axial equilibrium on the predator axis. For 0 < bn  dð1 þ anÞ < h the system admits only the axial equilibrium Ey on the predator axis. Here also location of the initial state relative to the stable manifold of the axial equilibrium Ey on the predator axis decides the eventual state of the system. Having the initial state above the stable manifold drives the predator to increase monotonically and simultaneously the prey gets driven to extinction asymptotically. Solutions initiating below the stable manifold suffer predator extinction in a finite time (Fig. 4(b)). Next, we consider the case where in absence of additional food the growth rate of predator is greater than the death rate, that is b > d. Here also three different parametric conditions need to be considered namely (iv) a þ 1n < 1 < bd, (v) 1 < bd < a þ 1n and (vi) 1 < a þ 1n < bd. If the additional food characteristic satisfies the condition a þ 1n < 1 < bd then the predator–prey system admits one axial equilibrium and can admit at most one or two interior equilibrium point/s if c < n or c > n respectively. In Region C2

Region C1 20 25

20

15

y→

y→

15 10

10

5 5

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0

0.7

0

0.1

0.2

0.3

0.4

x→

0.5

0.6

0.7

x→

2.4

2.2

2.2

2

y→

y→

2

1.8

1.8 1.6

1.4

1.6

0.4

0.45

0.5

0.55

x→

0.6

0.65

1.2 0.3

0.35

0.4

0.45

0.5

0.55

0.6

0.65

0.7

0.75

x→

Fig. 6. Phase portraits corresponding to the domains C 1 ; C 2 of the bifurcation diagram presented in Fig. 3. Here the values of the parameters are d ¼ :97; h ¼ :65 (C 1 ), d ¼ :9595; h ¼ :65 (C 2 ). The other parameter values are a ¼ 2; b ¼ 3; c ¼ 1:2, n ¼ 1.

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

209

case of c < n, the system has one interior equilibrium point Eh2 or no interior equilibrium if h is either below or above the line

D2 ¼ 0 respectively (Fig. 2(b)). If c > n, the system possesses one interior equilibrium point (Eh2 ) or two interior equilibrium points (Eh1 & Eh2 ) or no interior equilibrium point if h lies either below the line D2 ¼ 0 or in between D2 ¼ 0 and the curve D3 ¼ 0 or above the curve D3 ¼ 0 respectively (Fig. 2(a)). Now when h is bellow the line D2 ¼ 0 or above the curve D3 ¼ 0 the dynamics will be same as it was under the parametric restriction (iii). But if h lies between the line D2 ¼ 0 and the curve D3 ¼ 0 the predator–prey system fails to coexist (Region C=C 3 in Fig. 5). Now consider the case where the food characteristic satisfies the parametric restriction 1 < bd < a þ 1n. In this situation the system does not posses any axial equilibrium point. Under the above restriction if the environmental support for the prey (c) satisfies the restriction bn  dð1 þ anÞ þ cðb  dÞ < 0 then the predator goes to extinction in a finite time under constant harvesting as there is no equilibrium point. On the other hand if bn  dð1 þ anÞ þ cðb  dÞ > 0. Presence of two interior equilibria Eh1 and Eh2 for h lying below the curve D3 ¼ 0 and Eh2 being saddle, facilitates occurrence of Hopf, homoclinic and BT bifurcations in the system with variations in the parameters (Fig. 3). In this case, if ðd; hÞ is in B1 then the system coexists locally (Fig. 5), if h is in B2 then the system exhibits an oscillatory coexistence (basin of attraction in both the cases is bounded by the stable manifolds of the saddle point Eh2 (Fig. 5)), if h is in B3 a path of the system reaches the prey axis in a finite time (driven by the stable and unstable manifolds of the saddle interior equilibrium (Fig. 5)). Hence, the eventual state of the system is strongly dependent on the initial state of the system. Stable or oscillatory coexistence of predator and prey or the extinction of predator in a finite time depends on their initial states (Fig. 5)). Increasing the harvesting rate beyond ½bndð1þanÞþcðbdÞ2 4cðbdÞ

results in finite time extinction of the predator (Fig. 4(a)). From Table 1 and Fig. 3 it can be inferred that

as the magnitude of c increases the region of stability decreases. However it is possible to choose an appropriate value for the harvesting rate h (depending upon the magnitude of c) such that h lies in the region B1 or B2 and the predator–prey system coexists locally. Finally we consider the parametric restriction 1 < a þ 1n < bd. In this case the system admits one axial equilibrium and may admit one or two or no interior equilibrium depending upon the harvesting rate. If harvesting rate is comparatively low such that h < bn  dð1 þ anÞ holds, then the dynamics of the system is qualitatively similar to what we have observed under the case of bd < 1 (Fig. 4(c)). If the harvest rate (h) crosses the line D2 ¼ 0, and c satisfies the restriction 0 < bn  dð1 þ anÞ < cðb  dÞ, transcritical bifurcation takes place at the axial equilibrium point on the predator axis when D2 ¼ 0, and the axial equilibrium point Ey turns into saddle point for h lying above the line D2 ¼ 0. If h lies between

D2 ¼ 0 and the curve D3 ¼ 0 the system admits two interior equilibria Eh1 and Eh2 with Eh2 being saddle point. Here Eh1 is stable if ðd; hÞ is in the region C 1 and the coexistence of the prey and the predator depends upon their initial states (Fig. 6), Eh1 looses stability through Hopf-bifurcation and a stable limit cycle appears for d and h chosen in the region C 2 which signifies the oscillatory coexistence of the prey and predator (Fig. 6), the stable limit cycle disappears through homoclinic-bifurcation when the parameter values enter the region C 3 and as a result either the predator or the prey goes to extinction depending upon whether the system initiates from below or above the stable manifold of Eh2 (Fig. 5). In this case also as the magnitude of c gradually increases the region of coexistence gradually decreases and when it crosses a threshold value (c > c1 ) the system fails to exhibit the coexistence scenario (see Table 1) in the dh plane.The qualitative behavior of solutions of the system for 2

anÞþcðbdÞ h > ½bndð1þ with 0 < bn  dð1 þ anÞ < cðb  dÞ or h > bn  dð1 þ anÞ with bn  dð1 þ anÞ > cðb  dÞ are similar to the 4cðbdÞ

case when bd < 1 (Fig. 4(b)). Here, it is important to note the behavior of solutions of the system when the predator is provided with additional food. Irrespective of the productive nature of the system without additional food (i.e., position of bd relative to 1), provision of additional food introduces Ey on the predator axis and this equilibrium plays a crucial role to determine eventual state of the system depending upon the initial state. If the system is originally unproductive, it is possible to choose appropriate additional food, so that the predator grows with positive gradient in spite of harvesting. On the other hand if the system without additional food is productive (bd > 1Þ, but the system without the provision of additional food fails to coexist then appropriate introduction of additional food with 1 þ an P n to the system can bring coexistence into the system. Before concluding this section we provide two preliminary examples to support our claim presented above. Let us consider the system (2.4) without the provision of additional food that is a ¼ 0; n ¼ 0, and the other parameter values are b ¼ 1; c ¼ 1:45; d ¼ 0:4; h ¼ 0:1. For this choice the system does not admit any interior equilibrium point and so the predator goes to extinction in finite time. Now without changing any parameter value, if we choose a ¼ 1 and n ¼ 0:5 that is the system is now provided with additional food then the system parameters satisfy the parametric restrictions 1 þ an  n > 0, an , 0 < n < c2 < 1 þ an < c1 and c2 < c < 1 þ an. Also the restriction given by (v) and the condition b < 1þ1þann 2

anÞþcðbdÞ are satisfied, the parameters ðd; hÞ lies in the region B1 . Consequently the system possesses two interior h < ½bndð1þ 4cðbdÞ

equilibrium points Eh1 ð0:396; 1:378Þ, which is locally asymptotically stable and Eh2 ð1:221; 0:43Þ, a saddle point. There is no axial equilibrium point in this case. Thus depending upon the initial state, the prey and predator either sustain for all future times or the predator goes to extinction in a finite time. an If the additional food is chosen with a ¼ 1:49 and n ¼ 1, then we have 1 þ an  n > 0; b < 1þ1þann , 0 < n < c2 < 1 þ an < c1 ; n < c < c2 . Now ðd; hÞ lies in the region C 1 . Hence the system has one saddle axial equilibrium point Ey ð0; 62:25Þ as bn  dð1 þ anÞ > 0 and the parametric restriction (vi) are satisfied. We find two interior equilibria given by

210

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

Eh1 ð0:184; 2:334Þ and Eh2 ð1:259; 0:494Þ. Here Eh1 is locally asymptotically stable and Eh2 is a saddle point. In this situation either the prey and predator can coexist for all future times or one of them can go to extinction depending upon their initial states. Based upon these simple examples we can understand that the appropriate choice of the quantity and the quality of additional food (which in turn determines the level of additional food characteristic) and the initial state of the prey and predator lead to the most favorable situation where we find eventual coexistence scenario for all future times. 7. Discussion and conclusions In this article, we have studied the global dynamics of predator–prey system which is analogous to the model considered in [31], where predator is provided with additional food as well as harvested at a constant rate. To our knowledge, there exists no literature that presents consequences of harvesting a additional food provided predator–prey system. This is the first attempt to study the effect of constant harvesting on the foresaid predator–prey system. The assumption of constant harvesting of the predator is similar to the assumption made in [6,21,39]. The system is analyzed under various natural conditions on the parameters to identify local and global bifurcations involved in the dynamics. This study is useful in identifying appropriate characteristics for the additional food and the level constant harvesting to guarantee continuous availability of predator. This work presents all possible bifurcation scenarios that can arise in the dynamics of the considered model. We observe that a combination of the quality and quantity of the additional food that is available to the predator plays a key role in determining the amount of harvesting that guarantees continuous supply of predator. Apart from providing right kind of additional food to the predator, the initial state of the predator–prey system is also equally important for the sustenance of the predator. Various bifurcation diagrams presented here are consequences of variations in the magnitude of additional food characteristic with respect to the natural productivity of the prey. The features that are observed through the dynamical analysis are of practical importance and also provide guideline to determine the appropriate amount of additional food to have continuous supply of the predator with which the demand in the society can be met. Important ecological interpretations of our findings are now in order. If the predator goes extinct under harvesting because of its higher death rate in the absence of additional food then it is not possible to maintain the predator at a positive level by providing it with additional food. Essentially, in the absence of additional food if births in the predator by consuming the prey are not compensated by the deaths then providing additional food can not improve the situation. If the system without additional food is such that the birth rate of the predator is higher than its death rate then additional food promotes continuous survival of the predator for appropriate choice of constant harvesting in predators. Here it is important to note that the environmental carrying capacity plays a significant role in determining the additional food such that the system coexists. As we observe that for certain provision of additional food if the magnitude of environmental carrying capacity gradually increases the region of coexistence squishes. Thus the choice of the quantity of additional food will also depend on the environmental carrying capacity to have stable or oscillatory coexistence of both the species although the predators are harvested at a constant rate. Finally, we conclude that any predator–prey system can be turned into a good source of commercially important predator by providing the predator with additional food of desirable amount and characteristic. Coexistence of predator and prey in presence of additional food to predator depends on the magnitude of the additional food characteristic and harvesting intensity. Appropriate choice for the additional food characteristic and/or reasonable harvesting rates promote coexistence of predator and prey. Therefore, systems intending to promote coexistence need to regulate the harvesting activity and the provision of the additional food with desirable nutritional value. Unbounded growth in predator under certain circumstances suggests a necessity to include more realism in the model formulation and we hope that periodic harvesting effort or density dependent linear/nonlinear harvesting can lead to more realistic as well as challenging research problems. Another interesting issue of research will be to design a state feedback controller to eliminate the induced bifurcations. Analysis of such types of reformed models shall be our next endeavor. Acknowledgement We highly appreciate Prof. A. Hastings (University of California, Davis, CA) for his valuable suggestions to prepare the revised version. Authors are thankful to the anonymous reviewers for their careful reading and constructive suggestions. The work of M.S. is supported by the Council of Scientific and Industrial Research, India through their senior research fellowship. References [1] L.J.S. Allen, An Introduction to Mathematical Biology, Springer-Verlag, New York, 2006. [2] M. van Baalen, V. Krivan, P.C.J. van Rijn, M.W. Sabelis, Alternative food, switching predators, and the persistence of predator–prey systems, Am. Nat. 157 (5) (2001) 512–524. [3] J.R. Beddington, R.M. May, Maximum sustainable yields in systems subject to harvesting at more than one trophic level, Math. Biosci. 51 (1980) 261– 281. [4] J.R. Beddington, J.G. Cooke, Harvesting from a prey–predator complex, Ecol. Model. 14 (1982) 155–177. [5] F. Brauer, C. Castilo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer-Verlag, New York, 2012.

M. Sen et al. / Applied Mathematics and Computation 250 (2015) 193–211

211

[6] F. Brauer, A.C. Soudack, Stability regions and transition phenomena for harvested predator–prey systems, J. Math. Biol. 7 (1979) 319–337. [7] F. Brauer, A.C. Soudack, Stability regions in predator–prey systems with constant-rate prey harvesting, J. Math. Biol. 8 (1979) 55–71. [8] F. Brauer, A.C. Soudack, Coexistence properties of some predator–prey systems under constant rate harvesting and stocking, J. Math. Biol. 12 (1981) 101–114. [9] J. Carr, Applications of Centre Manifold Theory, Springer-Verlag, New York, 1981. [10] V. Christensen, Managing fisheries involving predator and prey species, Rev. Fish Biol. Fish. 6 (1996) 417–442. [11] C.W. Clark, Bioeconomic Modelling and Fisheries Management, Wiley, New York, 1985. [12] C.W. Clark, in: Mathematical Bioeconomics, The Optimal Management of Renewable Re-sources, John Wiley & Sons, New York, 1990. [13] G. Dai, M. Tang, Coexistence region and global dynamics of a harvested predator–prey system, SIAM J. Appl. Math. 58 (1998) 193–210. [14] O. Flaaten, On the bioeconomics of predator–prey fishing, Fish. Res. 37 (1998) 179–191. [15] R.P. Gupta, M. Banerjee, P. Chandra, The dynamics of two-species allelopathic competition with optimal harvesting, J. Biol. Dyn. 6 (2012) 674–694. [16] J.D. Harwood, J.J. Obrycki, The role of alternative prey in sustaining predator population, in: M.S. Hoddle (Eds.), Proceedings of Second International Symposium on Biological Control of Arthropods 2, 2005, pp. 453–462. [17] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Application of Hopf-Bifurcation, Cambridge University Press, Cambridge, 1981. [18] S.L. Hill, E.J. Murphy, K. Reid, P.N. Trathan, A.J. Constable, Modelling Southern Ocean ecosystems: krill, the food-web, and the impacts of harvesting, Biol. Rev. 81 (2006) 581–608. [19] R.D. Holt, Predation, apparent competition, and the structure of prey communities, Theor. Popul. Biol. 12 (1977) 197–229. [20] R.D. Holt, J.H. Lawton, The ecological consequences of shared natural enemies, Ann. Rev. Ecol. Syst. 25 (1994) 495–520. [21] J. Huang, Y. Gong, S. Ruan, Bifurcation analysis in a predator–prey model with constant-yield harvesting, DCDS-B 18 (2013) 2101–2121. [22] J. Huang, Y. Gong, J. Chen, Multiple Bifurcations in a predator–prey system of Holling and Leslie type with constant-yield prey harvesting, Int. J. Bifurcation Chaos 23 (2013), http://dx.doi.org/10.1142/S0218127413501642. [23] Y.A. Kuznetsov, Elements Applied Bifurcation Theory, second ed., Springer, 1997. [24] W.W. Murdoch, J. Chesson, P.L. Chesson, Biological control in theory and practice, Am. Nat. 125 (3) (1985) 344–366. [25] L. Perko, Differential Equations and Dynamical Systems, Springer, New York, 2001. [26] G. Peng, Y. Jiang, C. Li, Bifurcations of a Holling-type II predator–prey system with constant rate harvesting, Int. J. Bifurcation Chaos 19 (2009) 2499– 2514. [27] B.S.R.V. Prasad, M. Banerjee, P.D.N. Srinivasu, Dynamics of additional food provided predator–prey system with mutually interfering predators, Math. Biosci. 246 (1) (2013) 176–190. [28] P.C.J. van Rijn, van Houten, M.W. Sabelis, How plants benefit from providing food to predators even when it is also edible to herbivores, Ecology 83 (2002) 2664–2679. [29] S. Ruan, D. Xiao, Global analysis in a predator–prey system with nonmonotonic functional response, SIAM. J. Appl. Math. 61 (2001) 14451472. [30] J. Sotomayor, Generic bifurcations of dynamical systems, in: M.M. Peixoto (Ed.), Dynamical Systems, Academic Press, New York, 1973. [31] P.D.N. Srinivasu, B.S.R.V. Prasad, M. Venkatesulu, Biological control through provision of additional food to predators: a theoretical study, Theor. Popul. Biol. 72 (2007) 111–120. [32] P.D.N. Srinivasu, B.S.R.V. Prasad, Time optimal control of an additional food provided predator–prey system with applications to pest management and biological conservation, J. Math. Biol. 60 (2010) 591–613. [33] P.D.N. Srinivasu, B.S.R.V. Prasad, Erratum to: time optimal control of an additional food provided predator–prey system with applications to pest management and biological conservation, J. Math. Biol. 61 (2010) 591–613. [34] P.D.N. Srinivasu, B.S.R.V. Prasad, Role of quantity of additional food to predators as a control in predator–prey systems with relevance to pest management and biological conservation, Bull. Math. Biol. 73 (2011) 2249–2276. [35] D. Stiefs, G.A.K. van Voorn, B.W. Kooi, U. Feudel, T. Gross, Food quality in producer-grazer models: a generalized analysis, Am. Nat. 176 (2010) 367–380. [36] G.A. van Voorn, D. Stiefs, T. Gross, B.W. Kooi, U. Feudel, S.A. Kooijman, Stabilization due to predator interference: comparison of different analysis approaches, Math. Biosci. Eng. 5 (3) (2008) 567–583. [37] J.T. Wootton, The nature and consequences of indirect effects in ecological communities, Am. Rev. Ecol. Syst. 25 (1994) 443–466. [38] D. Xiao, L.S. Jennings, Bifurcations of a ratio-dependent predator–prey with constant rate harvesting, SIAM J. Appl. Math. 65 (2005) 737–753. [39] D. Xiao, S. Ruan, Bogdanov–Takens bifurcation in Predator–Prey systems with constant rate Harvesting, Fields Inst. Commun. 21 (1999) 493–506.