Bifurcation analysis of a Holling type II predator-prey model with refuge

Bifurcation analysis of a Holling type II predator-prey model with refuge

Journal Pre-proof

Bifurcation analysis of a Holling type II predator-prey model with refuge Debasis Mukherjee, Chandan Maji PII: DOI: Reference:

S0577-9073(20)30032-0 https://doi.org/10.1016/j.cjph.2020.02.012 CJPH 1075

To appear in:

Chinese Journal of Physics

Received date: Revised date: Accepted date:

19 August 2019 27 November 2019 14 February 2020

Please cite this article as: Debasis Mukherjee, Chandan Maji, Bifurcation analysis of a Holling type II predator-prey model with refuge, Chinese Journal of Physics (2020), doi: https://doi.org/10.1016/j.cjph.2020.02.012

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier B.V. on behalf of The Physical Society of the Republic of China (Taiwan).

Highlights • A predator-prey model with prey refuge is considered. • The amount of refuge is proportional to the number of encounters between prey and predator. • Different kind of bifurcation analysis is discussed. • Some numerical simulations are performed to illustrate the results.

1

Bifurcation analysis of a Holling type II predator-prey model with refuge Debasis Mukherjeea,∗, Chandan Majia a

Department of Mathematics, Vivekananda College, Thakurpukur, Kolkata-700063, India

Abstract In this article we discuss bifurcation analysis of a predator-prey model allowing for prey refuge. Here the amount of prey in refugia is proportional to encounters between prey and predator. The equilibrium points and their local stability have been analyzed. The local bifurcation (Hopf, saddle-node, transcritical, Bogdanov-Takens) have been discussed. Stability of the limit cycle emerging through Hopf-bifurcation is given. Persistence criterion of the model is developed. Some numerical simulations are performed to validate the obtained results. Keywords: predator-prey; refuge; stability; bifurcation; persistence. 1. Introduction The study of predator-prey interaction becomes an important area of research due to its universal existence. A significant numbers of papers have been published to understand the dynamics on the predator-prey systems. It is worth mentioning that prey refuges is an important ecological factor that influence systems dynamics. The risk of predation allows the prey population to take refuge. Prey’s response to predator attack using refuge is observed in a wild range of systems [1, 2, 3, 4, 5, 6, 7, 8]. To understand the dynamics of predator-prey systems with prey refuge mathematical models are formulated and analyzed [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27] . The effect of prey refuges on the dynamics of interacting populations is complex in nature. It has mainly two roles. Firstly, it affects positively the growth of prey and negatively that of predators because prey mortality is reduced due to decrease in predation success. The second one may be trade-offs and bi-products of the prey refuging sometimes could be useful or harmful for the involved populations. For example, reduce foraging activity in response to predation risk can lower the prey growth rate, reproduction and ultimately population size [2, 28, 29]. Also prey in refuges give more effort to obtain resources [30, 31] where as the opposite hold outside the ∗

Corresponding author.Ph.:+919007985012 Email addresses: [email protected] (Debasis Mukherjee), [email protected] (Chandan Maji)

Preprint submitted to Chinese Journal of Physics

refuge [32]. In prey-predator model refuge is taken either in constant proportion or in constant numbers [11, 14, 16, 18, 19, 20, 25]. These studies revealed that refugia protecting a constant number of prey lead to a stronger stabilizing effect on population dynamics, in comparison to refugia protecting a constant proportion of prey. It is noted that the dynamics of predator-prey systems may be affected by the preys using spatial refuges. Such type of refuges are observed where environmental heterogeneity provides lessaccessible sides for predators. By gathering group prey can also avoid predation. In modeling purpose two types of refuges are usually taken. The quantity of prey in cover xr is proportion to prey density x that is xr = mx and the other one for which the refuge population is a fixed quantity xr = m . Ruxton [24] considered xr = my where y is predator density. In the present article we are interested to investigate the situation where prey in refugia is proportional to encounters between prey and predators that is xr = mxy. In section 2 we will first propose a two dimensional predator-prey model incorporating prey refugia which depend on the encounters between the prey and predator. We give a thorough analysis of the model to study bifurcations and related dynamics, including Hopf, Saddle-Node, Transcritical, Bogdanov-Takens bifurcations which exhibit the complex dynamics. 2. Model Formulation We choose a habitat where prey and predator living together.To avoid predation, prey population takes refuge. Let x(t), y(t) denote for the prey and predator population. Here we assume that a fraction of prey population making refuge use for avoiding predation is proportional to encounters of prey and predators, that is, xr = mxy with m > 0. The growth rate of prey is of logistic type. Predation process follows Holling type II response function.     x c(1 − my)y dx = x r 1− − dt k 1 + ax(1 − my)   dy pc(1 − my)x = y −d+ dt 1 + ax(1 − my) with initial conditions x(0) > 0, y(0) > 0.

3

(1)

Parameters r k m c p a d

Table 1. Parameters Description intrinsic growth rate of the prey population carrying capacity of the environment prey refuge constant predation coefficient conversion efficiency of the predator half saturation constant death rate of the predator

2.1. Existence and uniqueness For t > 0, we take X ≡ (x, y)T , f : R2 → R2 , f = (f1 , f2 )T , system (1) can be = f (X). Here fi ∈ C ∞ (R2 ) for i = 1, 2 where f1 (x, y) = x{r(1 − xk ) − written as dX dt c(1−my)y pc(1−my)x }, f2 (x, y) = y{−d + 1+ax(1−my) }. Throughout this paper we assume that 1+ax(1−my) my < 1. Clearly, the vector function f is a smooth function of the variables (x, y) in 2 so, the quantity {1 + ax(1 − my)} can not be zero. Thus the positive quadrant R+ local existence and uniqueness of the solution hold. 3. Positivity and Boundedness of Solutions In this section, we shall first show positivity and boundedness of solutions of system (1). These are very important so far as the validity of the model is related. We first study the positivity. 2 Lemma 1. All solutions (x(t), y(t)) of system (1) with initial value (x0 , y0 ) ∈ R+ , remains positive for all t > 0.

Proof. Since right hand side of system (1) is continuous and Lipschitzian on C, thus solution of system (1) with initial conditions exists and unique on [0, ρ) where 0 < ρ < ∞. Integrating the first equation of system (1), we obtain, Rt x(t) = x0 exp{ 0 [r(1 −

x(s) ) k

−

c(1−my(s))y(s) ]ds}, 1+ax(s)(1−my(s))

which is always positive as x0 > 0. Similarly from second equation of system (1) we obtain, Rt y(t) = y0 exp{ 0 [−d +

pc(1−my(s))x(s) ]ds}, 1+ax(s)(1−my(s))

2 which is always positive as y0 > 0. Hence the interior of R+ is an invariant set of system (1).

2 Lemma 2. The set D = {(x, y) ∈ R+ : 0 < W = px + y ≤ µη } in a region of the attractors for all solutions initiating in the interior of the positive quadrant, where 2 0 < η < d and µ = kp(r+η) . 4r

4

Proof. Let W = px + y and η > 0 be a constant. Then   dW x + ηW (t) = px r(1 − ) + η + (η − d)y dt k

(2) 

Now choose η such that 0 < η < d. As, maximum value of px r(1 − kp(r+η)2 . 4r dW dt

Then eq. (2) implies that

+ ηW ≤

kp(r+η)2 4r

x ) k

+η



is

= µ.

Therefore, using the differential inequality [33], we obtain 0 < W (t) ≤

µ(1−e−ηt ) η

+ W (0)e−ηt .

Taking limit when t → ∞, we have, 0 < W (t) ≤ the Lemma.

µ η

for sufficiently large, proving

4. Equilibria The system of equations (1) has the following equilibria : (i) The trivial equilibrium E0 (0, 0), which is always exist. (ii) The predator free equilibrium E1 (k, 0) which is always feasible. (iii) The coexistence equilibrium E2 (x∗ , y ∗ ) where y ∗ = rpx∗ (k − x∗ )(kd)−1 and x∗ is the positive root of the equation mrp(pc − ad)x3 − mrpk(pc − ad)x2 + kd(pc − ad)x − d2 k = 0

(3)

E2 is feasible when x∗ < k.

The eqn. (3) can be written in the form as a0 x3 + 3a1 x2 + 3a2 x + a3 = 0 , a2 = where, a0 = mrp(pc − ad), a1 = − mrpk(pc−ad) 3

kd(pc−ad) 3

and a3 = −d2 k.

Now, define G = a20 a3 − 3a0 a1 a2 + 2a31 , H = a0 a1 − a21 . Again, the eigenvalues of the Jacobian matrix of system (1) evaluated at E0 are r, −d establishing that it is a saddle, unstable along the direction of x-axis as the eigenvalue in this direction is positive. The eigenvalues of the Jacobian matrix of pck pck system (1) at E1 are −r, 1+ak − d and hence it is stable if 1+ak < d. 5

4.1. The structure of the coexistence equilibria Theorem 1. Assume that pc − ad > 0. Equation (3) has (i) a unique positive root x∗ if G2 + 4H 3 > 0. (ii) Three positive roots, denoted by x∗1 , x∗2 , x∗3 if G2 + 4H 3 < 0. (iii) Two positive roots of which two are equal and denoted by x∗11 if G2 + 4H 3 = 0. (iv) Three equal roots denoted by x∗111 if G = H = 0. where 2 2 2 2 G = − m r p k(pc−ad) {27d2 + k(pc − ad)(2mrpk − 9d)}, 27 2 H = mrpk(pc−ad)9 (3d−mrpk) . Proof. The cubic eqn. (3) has always a positive root as the last term of the cubic is negative. Since pc − ad > 0, by Descarte’s rule of signs, eqn. (3) has no negative root. When G2 + 4H 3 > 0, eqn.(3) has two imaginary roots. In that case, eqn. (3) has also a unique positive root x∗ . Again, if G2 + 4H 3 < 0, eqn. (3) has three positive roots. When G2 + 4H 3 = 0, eqn. (3) has two positive roots and one root is zero. Lastly, when G = H = 0, eqn. (3) has three equal positive roots. We note that H > 0 if 3d > mrpk. So one can find unique positive root x∗ if pck pck < d < 1+ak . The inequality d < 1+ak implies that x∗ < k. We also observed that if pc ≤ ad then x∗ cannot be positive. So, for existence of positive equilibrium it should be necessary pc > ad. mrpk 3

4.2. Local bifurcation analysis Theorem 2. System (1) undergoes a transcritical bifurcation around E1 with respect pck to the parameter d if d = 1+ak . Proof. The Jacobian matrix J of system (1) around E1 is given by   ck −r − 1+ak J(E1 ) = pck 0 −d + 1+ak Clearly, one of the eigenvalues of J(E1 ) is negative and the other will be zero if pck ¯ If v = (v1 , v2 )T and w = (w1 , w2 )T denotes the eigenvectors corre= d. d = 1+ak sponding to the eigenvalue zero of the matrix J(E1 ) and J(E1 )T respectively, then ck v2 = 0 and w = (0, 1)T . we get v = (v1 , v2 )T where rv1 + 1+ak ¯ = 0. Note that wT [fd (E1 , d)] Hence system (1) does not attain any saddle-node bifurcation around E1 . ¯ = −v2 6= 0. Again wT [Dfd (E1 , d)v] ¯ Further wT [D2 f (E1 , d)(v, v)] = −

2pcv22 {ck(1+rk)+rkm} r(1+ak)4

6

6= 0.

So, according to Sotomayor’s theorem [34], system (1) has a transcritical bifurcation at E1 with respect to the parameter d. Theorem 3. System (1) undergoes a saddle-node bifurcation around E2 with recay ∗ (1−my ∗ )2 ∗ ∗ ∗ 2 spect to the parameter m if kr = [1+ax ∗ (1−my ∗ )]2 and 2my = 1 + ax (1 − my ) . Proof. The Jacobian matrix J of system (1) around E2 is given by J(E2 ) =

cay ∗ (1−my ∗ )2 } [1+ax(1−my ∗ )]2 pc(1−my ∗ )y ∗ [1+ax∗ (1−my ∗ )]2

x∗ {− kr +

−x

∗ c{1−2my ∗ +ax∗ (1−my ∗ )}

[1+ax∗ (1−my ∗ )]2 ∗ y∗ − [1+axmpcx ∗ (1−my ∗ )]2

!

Let m ¯ be the value of m such that the matrix J(E2 ) has a simple zero eigenvalue at m = m. ¯ By the assumption of the theorem, one of the eigenvalues of J(E2 ) is negative and the other will be zero. If v = (v1 , v2 )T and w = (w1 , w2 )T denotes the eigenvectors corresponding to the eigenvalue zero of the matrix J(E2 ) and J(E2 )T respectively, then we get v = (v1 , v2 )T where (1 − my ∗ )v1 = mv2 x∗ and w = (1, 0)T . Note that wT [fm (E2 , m)] ¯ =

cx∗ y ∗ [1+ax∗ (1−my ∗ )]2

6= 0.

Hence system (1) may admit a saddle-node bifurcation around E2 . Again wT [Dfm (E2 , m)v] ¯ 6= 0, by the assumption of the theorem. Thus system (1) has a saddle-node bifurcation around E2 with respect to the parameter m.

4.3. Numerical example for saddle-node bifurcation 7 Choose r = 10 , c = 8, k = 73 , p = 12 , d = 0.2, a = 15, m = 45 . Since real data is not available to us at present so we have chosen the values hypothetically. Then system (1) has a interior equilibrium point (1, 1). With these choice of parameters, all the conditions of Theorem 3 are satisfied. So system (1) admits a saddle-node bifurcation around the equilibrium point (1, 1) with respect to the parameter m (see Fig. 1.).

7

predator

1.5

1

0.5

0 0

0.5

1 m

1.5

2

Fig. 1. Saddle-node bifurcation occurs around interior equilibrium point E2

4.4. Stability and Hopf-bifurcation of coexistence equilibria : Theorem 4. (i) The equilibrium point E2 is locally stable if 1 + ax∗ (1 − my ∗ )2

r k

>

cay ∗ (1−my ∗ )2 [1+ax∗ (1−my ∗ )]2

and 2my ∗ <

(ii) The system undergoes a Hopf-bifurcation with respect to the bifurcation param∗ {a(1−my ∗ )2 −mp} eter m around the equilibrium point E2 if kr = cy[1+ax and (1 − my ∗ ){1 − ∗ (1−my ∗ )]2 2my ∗ + ax∗ (1 − my ∗ )2 } > m2 px∗ y ∗ . Proof. (i) Assumptions of the theorem imply that T race J(E2 ) is negative and det J(E2 ) is positive. Therefore the eigenvalues of the matrix J(E2 ) have negative real parts. Hence the equilibrium point E2 is locally asymptotically stable. (ii) Now Trace of the matrix J(E2 ) is given by T race (J(E2 )) = x∗ {− kr +

cy ∗ {a(1−my ∗ )2 −mp} } [1+ax∗ (1−my ∗ )]2

which is zero by the assumption of the theorem. Further det J(E2 ) > 0 follows from the assumption. Now we calculate Here

d (T race dm

Again,

dx∗ dm

=

d (T race dm

J(E2 )). ∗

J(E2 )) = − dx {1 + dm k

rpx∗2 (k−x∗ ) 3mrpx∗2 −2mrpkx∗ +kd

ad kpc

+

(pc−ad)2 } rpc

> 0 if m ∈ (0, m∗ ) 8

where m∗ = Thus

3d . rpk

d (T race dm

J(E2 )) 6= 0.

This guarantees the existence of Hopf bifurcation around E2 . The condition kr > cay ∗ (1−my ∗ )2 indicates that the rate of change of prey’s specific growth rate with re[1+ax∗ (1−my ∗ )]2 spect to prey density decreases and the condition 2my ∗ < 1 + ax∗ (1 − my ∗ )2 indicates that the rate of change of prey’s specific growth rate with respect to predator density increases. The stability of the limit cycle is shown in next section. 5. Stability of limit cycle Now we investigate the nature of the limit cycle by calculating the first Lyapunov number σ at the positive equilibrium point E2 of system (1). Let x = u − x∗ , y = v − y ∗ . Then system (1) reduces to du dt

= a10 u + a01 v + a20 u2 + a11 uv + a02 v 2 + a30 u3 + a21 u2 v + a12 uv 2 + a03 v 3 + P (u, v),

dv dt

= b10 u + b01 v + b20 u2 + b11 uv + b02 v 2 + b30 u3 + b21 u2 v + b12 uv 2 + b03 v 3 + Q(u, v),

where a10 = −x∗ { kr − a20 = − kr +

cy ∗ (1−my ∗ )2 }, [1+ax∗ (1−my ∗ )]2

∗

∗

cay ∗ (1−my ∗ )2 , [1+ax∗ (1−my ∗ )]3

∗ )2

∗

∗

},

(1−my )−2my } a11 = − c{1+ax , [1+ax∗ (1−my ∗ )]3

a02 =

cmx∗ (1+ax∗ ) , [1+ax∗ (1−my ∗ )]3

a21 =

ac{1+ax∗ (1−my ∗ )2 −4my ∗ +3m2 y ∗2 } , [1+ax∗ (1−my ∗ )]4

a03 =

ax∗2 m2 c(1+ax∗ ) , [1+ax∗ (1−my ∗ )]4 ∗

∗

+ax (1−my a01 = −x∗ c{ 1−2my [1+ax∗ (1−my ∗ )]2

∗ 3 ∗

2

ca (1−my ) y a30 = − [1+ax ∗ (1−my ∗ )]3 ,

b10 =

∗

∗

∗

∗

{ax (1−my )−2my } a12 = − mcax[1+ax , ∗ (1−my ∗ )]3 pcy ∗ (1−my ∗ ) , [1+ax∗ (1−my ∗ )]2 ∗

∗

∗ 2

pcay (1−my ) b20 = − [1+ax ∗ (1−my ∗ )]3 ,

mpy b01 = − [1+axcx∗ (1−my ∗ )]2 ,

b11 =

pc{1+ax∗ (1−my ∗ )−2my ∗ } , [1+ax∗ (1−my ∗ )]3

pcx m(1+ax ) b02 = − [1+ax ∗ (1−my ∗ )]3 ,

b30 =

pca2 (1−my ∗ )3 , [1+ax∗ (1−my ∗ )]4

(1−my ) −4my +3m b21 = − pca{1+ax [1+ax ∗ (1−my ∗ )]4

b12 =

pcmax∗ {ax∗ (1−my ∗ )−2my ∗ } , [1+ax∗ (1−my ∗ )]3

P (u, v) =

P4

i+j

∗

∗

∗

∗2

2

∗ 2

∗

2 y∗ }

,

∗

ax m c(1+ax ) b03 = − [1+ax ∗ (1−my ∗ )]4 ,

aij ui v j and Q(u, v) =

P4

i+j bij u

i j

v .

Hence the first Lyapunov number σ for a planar system (as in described in [34]) 9

is given by σ=−

3π

3

2a01 ∆ 2

{[a10 b10 (a211 + a11 b02 + a02 b11 ) + a10 a01 (b211 + a20 b11 + a11 b02 ) + b210 (a11 a02 +

2a02 b02 ) − 2a10 b10 (b202 − a20 a02 ) − 2a10 a01 (a220 − b20 b02 ) − a201 (2a20 b20 + b11 b20 ) + (a01 b10 − 2a210 )(b11 b02 − a11 a20 )] − (a210 + a01 b10 )[3(b10 b03 − a01 a30 ) + 2a10 (a21 + b12 ) + (b10 a12 − a01 b21 )]}, where ∆ = a10 b01 − a01 b10 . As the expression for Lyapunov number σ is quite complicated, we are unable to say the sign of σ. Thus we present a numerical example. 5.1. Numerical example for Hopf-bifurcation Choose r = 23 , c = 8, k = 23 , p = 81 , d = 14 , a = 2, m = 21 . There exists a limit 8 7 cycle around (1,1). Note that the Lyapunov number σ = 0.0253 > 0. This implies that system (1) possesses a unstable limit cycle. It is to be noted here that positive equilibrium is unique. (See Fig. 2)

2

1.8

1.8

1.6

1.4

1.6

1.2

Predator

Predator

1.4

1.2

1

0.8

1 0.6 0.8

0.4

0.6

0.4

0.2

0

0.5

1

1.5

2

2.5

Prey

(a) m = 0.45

0

0

0.5

1

1.5 Prey

(b) m = 0.49

10

2

2.5

3

1.8

2

1.6

1.8

1.4 1.6

Predator

Predator

1.2

1

1.4

1.2

0.8 1 0.6

0.8

0.4

0.2

0

0.5

1

1.5 Prey

2

2.5

3

0.4

(c) m = 0.50

0.6

0.8

1

1.2

1.4 Prey

1.6

1.8

2

(d) m = 0.51

Fig.2. Phase portrait of system (1) with different values of m.

3.5 3

prey population

2.5 2 1.5 1 0.5 0 0

0.5

1 m

1.5

2

Fig. 3. Hopf-bifurcation diagram of system (1) with m as bifurcation parameter. 5 Now we choose r = 16 , c = 1, k = 5, p = 12 , d = 18 , a = 2, m = 21 . There exist two positive equilibrium points (1, 1) and (2, 1.5). We observe that there exist a limit

11

2.2

2.4

cycle which contain the point (1, 1) and the other point (2, 1.5) is locally stable. If we decrease the value of m the multiple equilibria disappear (see Fig. 4.). 2.5

2 1.8

2

1.6 1.4

Predator

1.2

Predator

1.5

1

1 0.8 0.6

0.5

0.4 0.2

0

0

0.5

1

1.5

2

2.5 Prey

3

3.5

4

4.5

0

5

0

0.5

1

1.8

1.8

1.7

1.7

1.6

1.6

1.5

1.5

1.4

1.4

1.3

1.2

1.1

1.1

1

1

0.9

0.9

1

1.2

1.4

1.6

2.5

3

3.5

4

1.3

1.2

0.8 0.8

2 Prey

(b) m = 0.40

Predator

Predator

(a) m = 0.20

1.5

1.8 Prey

2

2.2

2.4

2.6

2.8

0.8

1

1.5

2

2.5

3

Prey

(c) m = 0.50

(d) m = 0.60

Fig.4. Phase portrait of system (1) with different values of m.

6. Bogdanov - Takens bifurcation Conditions of Theorem 3 imply that system (1) undergoes a saddle-node bifurcation around E2 , if it exists and Trace J(E2 ) 6= 0. Now, we consider the case Trace J(E2 ) = 0. So one can find double zero eigenvalue of the Jacobian matrix J(E2 ) even it is not a null matrix. So, there is a possibility of occurring co-dimension 2 bifurca12

3.5

tion. One can easily show that the equilibrium point E2 is a cusp of co-dimension 2 ∗ )2 −my ∗ } whenever kr = cy{a(1−my and (1 − my ∗ ){1 − 2my ∗ + ax∗ (1 − my ∗ )2 } = m2 px∗ y ∗ . [1+ax∗ (1−my ∗ )]2 We apply the algorithm discussed in [35] to prove the non-degeneracy conditions of Bogdanov-Takens bifurcation. We have selected c and d as the bifurcation parameters as they are important from ecological point of view. We note that Bogdanov-Takens point (in short, BT-point) in the parameter space is the intersection of the saddlenode bifurcation curve and the Hopf- bifurcation curve. Let (c, d) = (c[BT ] + λ1 , d[BT ] + λ2 ) be a neighboring point of the B − T point (c[BT ] , d[BT ] ) with λi , i = 1, 2 sufficiently small. Then system (1) at (c, d) = (c[BT ] + λ1 , d[BT ] + λ2 ) is given by     x (c + λ1 )(1 − my)y dx = x r 1− − = f (x, y, λ1 ) dt k 1 + ax(1 − my)   dy p(c + λ1 )(1 − my)x = y − (d + λ2 ) + = g(x, y, λ2 ) dt 1 + ax(1 − my)

(4)

The system (4) is C ∞ smooth with respect to the variables x, y in a small neighbourhood of (c[BT ] , d[BT ] ). Define x1 = x − x∗ , x2 = y − y ∗ then system (4) reduces to dx1 = α00 + α10 x1 + α01 x2 + α20 x21 + α11 x1 x2 + α02 x22 + P1 (x1 , x2 ) dt dx2 = β00 + β10 x1 + β01 x2 + β20 x21 + β11 x1 x2 + β02 x22 + P2 (x1 , x2 ) dt

(5)

where ∗

∗

∗

∗

λ1 (1−my)y 1 x (1−my )y α10 = − [1+ax α00 = − λ1+ax ∗ (1−my ∗ ) , ∗ (1−my ∗ )]2 + ∗

∗

∗ )2 }

+ax (1−my α01 = −x∗ (c + λ1 ) {1−2my [1+ax∗ (1−my ∗ )]2 ∗

∗

(1−my )−2x α11 = − (c+λ1 ){1+ax [1+ax∗ (1−my ∗ )]3

β00 = −λ2 y ∗ + pλ1 x∗ (1−my ∗ ) , 1+ax∗ (1−my ∗ )

pλ1 (1−my ∗ )x∗ , 1+ax∗ (1−my ∗ )

∗

∗ 2

∗ y∗ }

, α20 = − kr +

, α02 =

β10 =

1 )ay (1−my ) β20 = − p(c+λ , β11 = [1+ax∗ (1−my ∗ )]3

cayx∗ (1−my ∗ )2 [1+ax∗ (1−my ∗ )]2

− kr x∗ ,

(c+λ1 )ay ∗ (1−my ∗ )2 , [1+ax∗ (1−my ∗ )]3

(c+λ1 )x∗ m(1+ax∗ ) , [1+ax∗ (1−my ∗ )]3

p(c+λ1 )y ∗ (1−my ∗ ) , [1+ax∗ (1−my ∗ )]2

β01 = −λ2 −

p(c+λ1 ){1+ax∗ (1−my ∗ )−2my ∗ } , [1+ax∗ (1−my ∗ )]3

(c+λ1 )x∗ mpy ∗ [1+ax∗ (1−my ∗ )]2

∗

Let us apply the affine transformation y1 = x1 , y2 = a10 x1 + a01 x2 we get = ξ00 (λ) + ξ01 (λ)y2 + ξ20 (λ)y12 + ξ11 (λ)y1 y2 + ξ02 (λ)y22 + Q1 (y1 , y2 ) 13

∗

1 )x m(1+ax ) β02 = − p(c+λ [1+ax∗ (1−my ∗ )]3

and P1 , P2 are the power series in (x1 , x2 ) with powers xi1 xj2 satisfying i + j ≥ 3.

dy1 dt

+

dy2 dt

= η00 (λ) + η10 (λ)y1 + η01 (λ)y2 + η20 (λ)y12 + η11 (λ)y1 y2 + η02 (λ)y22 + Q2 (y1 , y2 )

where ξ00 (λ) = α00 , ξ01 (λ) = α01 , ξ20 (λ) = ξ02 (λ) = αα02 2 , 01 α10 α01 + β11 ,

(α20 α201 −α11 α10 α01 +α02 α210 ) , a201

η00 (λ) = α10 α00 + α01 β00 ,

a11 −2a02 a10 , a01

η10 (λ) = β10 α01 − β01 α10 ,

η20 (λ) =

{α10 (α20 α201 −α11 α10 α01 +α02 α210 )+α201 (α01 β20 −β11 α10 )+β02 α210 α01 } , α201

η11 (λ) =

α10 (α11 −2α02 α10 −2β02 ) , α01

η02 (λ) =

ξ11 (λ) =

η01 (λ) =

α10 α02 +β02 α01 α201

and Q1 , Q2 are the power sums in (y1 , y2 ) with powers y1i y2j satisfying i + j ≥ 3. The non-degeneracy conditions of Bogdanov-Takens bifurcations [35] are 2ξ20 (0) + η11 (0) 6= 0, η20 (0) 6= 0. Thus we have the following theorem. Theorem 5. System (1) undergoes a Bogdanov-Takens bifurcation with respect to the bifurcation parameter c and d around the equilibrium point E2 where kr = cy{a(1−my ∗ )2 −my ∗ } and (1 − my ∗ ){1 − 2my ∗ + ax∗ (1 − my ∗ )2 } = m2 px∗ y ∗ and 2ξ20 (0) + [1+ax∗ (1−my ∗ )]2 η11 (0) 6= 0, η20 (0) 6= 0.

6.1. Numerical example for Bogdanov-Takens bifurcation Choose r = 1, k = 3, d[BT ] = 1, a = 21, c[BT ] = 12, p = 2, m = 21 , together with (x∗ , y ∗ ) = (1, 43 ). Here 2ξ20 (0) + η11 (0) = −0.5581 6= 0, η20 (0) = 0.1095 6= 0. For this set of parameter values, the unique equilibrium (1, 43 ) becomes a cusp of co-dimension 2.

14

2

1.6 1.4 1.2

Predator

Predator

1.5

1

1 0.8 0.6

0.5

0.4

0 1

1.5

2

2.5

Prey

3

3.5

0.2 0.5

4

(a) (λ1 , λ2 ) = (0, 0)

1

1.5

2

2.5

Prey

3

(b) (λ1 , λ2 ) = (1.81, 0.2)

2

Predator

1.5

1

0.5

0 0

1

2

Prey

3

4

(c) (λ1 , λ2 ) = (1.86, 0.2) Fig.5. Phase portraits of system (4): (a) when (λ1 = 0, λ2 = 0) the system converges to a stable equilibrium point where prey equilibrium density is low, (b) (λ1 , λ2 ) = (1.81, 0.2) the system converges into a stable equilibrium point but prey equilibrium density increases, (c) (λ1 , λ2 ) = (1.86, 0.2) a stable limit cycle bifurcates through Hopf-bifurcation around the equilibrium point .

7. Uniform persistence Biologically, persistence means the long term survival of all populations whatever may be the initial populations. Mathematically, it means that strictly positive solutions do not have omega limit points on the boundary of the non-negative cone. Now we state a result guaranteeing the uniform persistence of system (4). Theorem 6. Suppose that d <

pkc 1+ak

then system (1) is uniformly persistent.

Proof. To prove the theorem we use the ’average Lyapunov function’ and apply the theorem by Hutson [36]. Consider, the average Lyapunov function of the form : P (X) = xp1 y p2 , where, each 15

3.5

4

2 pi , i = 1, 2 is assumed positive. In the interior of R+ , one has 1 dP (X) P (X) dt

= ψ(X) = px1 dx + dt = p1 {r(1

p2 dy y dt − xk )

−

c(1−my)y } 1+ax(1−my)

+ p2 (−d +

pc(1−my)x }. 1+ax(1−my)

2 We have to show ψ(X) > 0 for all X ∈ bdR+ , for a suitable choice of p1 , p2 > 0 to prove uniform persistence of system (1). That is one has to satisfy the following conditions corresponding to the boundary equilibria E0 and E1 .

E0 : p1 r − p2 d > 0 pkc )>0 E1 : p2 (−d + 1 + ak

(6) (7)

Positivity of (7) follows from the assumption of the theorem. If we choose p1 and p2 in such a way, positivity of (6) will follow. This completes the proof of the theorem.

8. Discussion In this paper we have considered a predator-prey model with refuge. We have assumed here that the amount of prey in refugia is proportional to the encounters between prey and predator. This assumption is appropriate, as it justify a good and tractable approximation to observe behavior regarding inducible defence [37]. Holling type-II predator prey model with refuge have been studied in [11, 14, 16, 18, 19, 20, 25, 38]. All this work do not address the situation when the number of refugia is proportion between the predator and prey that is xr = mxy. Our analysis indicates multiple positive equilibria can exist. Moreover, complicated bifurcation namely Bogdanov-Takens bifurcation can be established. Previous models do not show such complexity. We note that increasing the prey refuges stabilizes the system. When refuge is low, we observed periodic oscillation. The model experiences different type of local bifurcation such as transcritical, saddlenode, Hopf bifurcation and Bogdanov-Takens bifurcation. These are important in ecological point of view in particular transcritical bifurcation convert prey extinction equilibrium into a unstable equilibrium point and a unstable positive equilibrium to a positive equilibrium. We have found an unique positive equilibrium point which is locally asymptotically stable under certain parametric restriction. Existence of Hopf-bifurcation is ensured when refuge parameter crosses a critical value m∗ and stability of the limit cycle has been tested by numerical example by estimating first Lyapunov number. We have also find two distinct positive equilibrium points one of which is locally asymptotically stable and the other one encircle a limit cycle arising from Hopf-bifurcation (see Fig 4. (c)). Coexistence of all populations are shown by invasion analysis. In absence of refuge a predator prey system with Holling type-II 16

functional response usually admits limit cycle behaviour emerging through Hopf bifurcation but if small amount of refuge is incorporated in the system it will stabilize the system (see Fig 2.(d)). Also number of coexistence equilibrium point varies in case of refuge system. For example the system without refuge admits only one coexistence equilibrium point whereas system with refuge may admits two or three equilibrium points. B-T bifurcation is not possible for the system without refuge whereas it is observed in refuge system. Acknowledgement: The authors are grateful to the anonymous reviewers and Editor for their helpful comments and suggestions for improving the paper.

References [1] S.S. Bell, E.D. McCoy, H.R. Mushinsky, Habitat structure: the physical arrangement of objects in space. Chapman and Hall, New York, New York USA, 1991. [2] J.S. Brown, B.P. Kotler, Hazardous duty pay and the foraging cost of predation. Ecology Letters. 7 (2004) 999-1014. [3] T. Caro, Antipredator defenses in birds and mammals. University of Chicago Press, Chicago, Illinois, USA, 2005. [4] W.E. Cooper, Jr., Theory successfully predicts hiding time: new data for the Sceloporus virgatus and a review. Behavioural Ecology. 20 (2009) 585-592. [5] S.L. Lima, Stress and decision making under the risk of predation: recent developments from behavioral, reproductive, and ecological perspective. Advances in the study of Behaviour, 27(1998) 215-290. [6] S.L. Lima, L.M. Dill, Behavioral decisions made under the risk of predation- a review and prospectus. Canadian Journal of Zoology. 68(1990)619-640. [7] Jin Hyuk Choi, Hyunsoo Kim, Rathinasamy Sakthivel, On certain exact solutions of diffusive predator-prey system of fractional order, Chinese Journal of Physics 54 (2016) 135-146. [8] A.E. Allahverdyan, S. G. Babajanyan, C. -K. Hu, Polymorphism in rapidly changing cyclic environment, Phys. Rev. E 100 (2009) 032401 17

[9] L. Chen and F. Chen, Global analysis of a harvested predator-prey model incorporating a constant prey refuge, Int. J. of Biomath. 3 (2010) 205-223. [10] F. Chen, L. Chen and X. Xie, On a Leslie Gower predator-prey model incorporating a prey refuge, Nonlinear. Anal.RWA. 10 (2010) 2905-2908. Nonlinear. Anal.RWA. 11 (2010) 246-252. [11] L. Chen, F. Chen and L. Chen, Qualitative analysis of a predator-prey model with Holling type II functional response incorporating a constant prey refuge, Nonlinear Anal. RWA. 11 (2010) 246-252. [12] R. Cressman, J. Garay, A predator-prey refuge system : Evolutionary stability in ecological systems, Theoret. Pop.Biol. 76 (2009) 248-257. [13] L.M. Dill, R. Houtman, The influence of distance to refuge on flight-initiation distance in the prey squirrel (Sciurus carolinensis), Canadian J. Zoology. 67 (1989) 232-235. [14] E. Gonz´ alez-Olivares, R. Ramos-Jiliberto,Dynamic consequences of prey refuges in a simple model system : More prey,fewer predators and enhanced stability, Ecol.Mod. 166 (2003) 135-146. [15] M. Haque, S. Rahman, E. Venturino, B.L. Li, Effect of a functional responsedependent prey refuge in a predator-prey model. Ecol. Complex. 20 (2014) 248-256. [16] W. Ko, K. Ryu, Qualitative analysis of a predator-prey model with Holling type-II functional response incorporating a prey refuge. J. Differ. Equ. 231 (2006) 534-550. [17] V. Krivan, Effects of optimal antipredator behavior of predator-prey dynamics: The role of refuges, Theor. Popul. Biol. 53 (1998) 131-142. [18] Z.H. Ma, W.L. Li, Y. Zhao, Effects of prey refuges on a predator-prey model with a class of functional responses: the role of refuges. Math. Biosci. 218 (2009) 73-79. [19] J.M. McNair, The effects of refuges on predator-prey interaction : a reconsideration, Theor.Popul.Biol. 29 (1986) 38-63.

18

[20] J.N. McNair, Stability effects of prey refuges with entry exist dynamics. J. Theor. Biol. 125 (1987) 449-464. [21] D. Mukherjee,The effect of prey refuges on a three species food chain model, J. Diff. Eqns. and Dyn. Syst. 22 (2014) 413-426. [22] D. Mukherjee, Study of refuge use on a predator-prey system with a competitor for the prey. Int. J.Biomath. (10)2 (2017) 1750023. DOI:10.1142/S1793524517500231. [23] D. Mukherjee, The effect of refuge and immigration in a predator-prey systems in the presence of a competitor for the prey. Nonlinear Anal., Real World Appl. 31 (2016) 277-287. [24] G.D. Ruxton, Short term refuge use and stability of predator-prey models, Theor. Popul. Biol. 47 (1995), 1-17. [25] A. Sih, Prey refuges and predator-prey stability, Theor. Popul. Biol. 31 (1987) 1-12. [26] G. Tang, S. Tang, R.A. Cheke, Global analysis of a Holling type II predator-prey model with a constant prey refuge. Nonlinear Dyn. 76 (2014) 635-647. [27] H. Wang, W. Morrison, A. Singh and H. Weiss, Modelling inverted biomass pyramids and refuges in ecosystems, Eco.Mod. 220 (2009) 1376-1382. [28] S. Creel, D. Christianson, Relationships between direct predation and risk effects. Trends in Ecology and Evolution. 23(2008) 194-201. [29] E.L. Preisser, D.I. Bolnick, M.F. Benard, Scared to death? The effects of intimidation and consumption in predator-prey interactions. Ecology. 86 (2005) 501-509. [30] J. Martin, P. Lopez, W.E. Cooper, When to come out from a refuge: balancing predation risk and foraging opportunities in an alpine lizard. Ethology. 109(2003) 77-87.

19

[31] L. Persson, P. Eckl¨ ov, Prey refuges affecting interactions between piscivorous perch and juvenile perch and roach. Ecology. 76 (1995) 70-81. [32] A. Sih, Prey uncertanity and the balancing of antipredator behavior and feeding needs. American Naturalist. 139(1992) 1052-1069. [33] G. Birkhoff, G.C. Rota, Ordinary Differential Equation ( Ginn and Co.Boston, 1982) [34] P. Lawrence, Differential Equations and Dynamical Systems, Springer-Verlag, New York, 1991. [35] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer-Verlag, New York, 1995. [36] V . Hutson, A theorem on average Lyapunov function, Monatsh. Math. 98 (1984) 267-275. [37] R. Ramos-Jiliberto and E. Gonz´ alez-Olivares, Relating behavior to population dynamics : a predator-prey metaphysiological model emphasizing zooplankton diel vertical migration as an inducible response, Ecological Modelling. 127 (2000) 221-233. [38] A. Sih, J.W. Petranka and L.B. Kats, The dynamics of prey refuge use : a model and tests with sunfish and salamander larvae, Am.Nat . 132 (1998), 463-483.

Conflict of Interest The authors declare that there is no conflict of interest in publishing this paper.

20