Dynamics of an additional food provided predator–prey system with prey refuge dependent on both species and constant harvest in predator

Physica A 534 (2019) 122301

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Dynamics of an additional food provided predator–prey system with prey refuge dependent on both species and constant harvest in predator Sudeshna Mondal, G.P. Samanta

∗

Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India

highlights • • • • •

Formulates a predator–prey model with prey refuge proportional to both species in presence of additional food. Dynamical behaviour of the model is analysed. Influence of gestation delay on the proposed model has been discussed. A modified model is proposed to incorporate constant harvest to the predator. Pictorial scenarios have been drawn through MATLAB and MATCONT to verify the analytical results.

article

info

Article history: Received 10 May 2019 Received in revised form 19 June 2019 Available online 9 August 2019 Keywords: Refuge Additional food Global stability Hopf-bifurcation Time delay Constant-yield harvesting

a b s t r a c t In the present work, we investigate on the dynamics of predator–prey model with prey refuge in proportion to both species and the effect of additional food provided to the predator. A Holling type II functional response in presence of additional food is taken to represent the interaction between prey and predator. Positivity and uniform boundedness of the solutions of the system are discussed analytically. Existence of equilibrium points and their stability criteria are analysed. Condition for occurrence of Hopf-bifurcation is derived analytically. It is found that the consumption rate of predator plays an important role in the dynamics of present system. Also, discrete time lag (gestation delay) in the proposed model is incorporated to observe the switching behaviour of the delay parameter through Hopf-bifurcation. Moreover, this work aims to analyse the dynamical behaviour of the considered predator–prey model when the predator is harvested at a constant rate. Numerical simulations are performed to verify all the analytical results using MATLAB and MATCONT. © 2019 Elsevier B.V. All rights reserved.

1. Introduction In both ecology and mathematical biology, the dynamic relationship between predator and prey has long been and will continue to be one of the dominant themes due to its universal existence and importance [1]. At first sight, mathematically these problems may be very simple but often they are very challenging and complicated. In last 50 years, there are huge development on dynamic behaviour of the predator–prey theory [2,3]. One of the classical applications of mathematical biology is to study the interactions between both the species by differential equation models because predator–prey ∗ Corresponding author. E-mail address: [email protected] (G.P. Samanta). https://doi.org/10.1016/j.physa.2019.122301 0378-4371/© 2019 Elsevier B.V. All rights reserved.

2

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

interaction has strong effect in the dynamical system, where predator plays an interesting role to maintain the food web structure. In general, in population system there are two types of predators: one is specialist (for example: small mustelids) and other is generalist (like foxes, common buzzards, cats etc.). Generalist predators must have different food sources and attack aggressively on their favourite food while specialist predators confide only on the specific food for healthy and successful life [4]. Moreover, the most fundamental idea in predator–prey interaction is ‘‘functional response’’ that describes the rate at which the number of prey is attacked by a predator. There are various types of functional responses namely Holling type I, Holling type II, ratio-dependent, Hassel Varley etc. From literature surveys [5,6], it can be mentioned that the most important and effective functional responses are Holling type I (Volterra) and Holling type II. The mathematical forms of Holling type I is f (x) = ax and Holling type II is f (x) = bax , where a is the maximum rate of predation, b is the half +x saturation constant and x is the prey biomass [7–9]. Provision of additional food to predator in the dynamics of predator–prey interaction has been studied intensively due to its wide applications in biological control and species conservation. It is concluded by some theoretical studies [10] that provision of additional food could increase target predator and decrease the attack rate of prey biomass. This reduction of attack rate on prey due to presence of additional food to the predator is known as apparent competition [11]. On the other hand, some empirical studies indicate that provision of additional food to predator need not always enhances predator biomass [12]. This switching nature between theory and observations led to in depth mathematical analysis of inclusion of additional food to predator–prey systems. Moreover, it is noticed that provision of additional food of high nutritive value could increase the predator population whereas additional food with low nutritive value could reduce predator biomass [13,14]. One of the most relevant concept on the dynamics of predator–prey model is to study the impact of spatial refuge which defends a constant proportion of prey from predation [15–17]. So, there exists a significant influence of prey refuge on the coexistence among predator and prey. Therefore, in theoretical ecology, the effect of prey refuge on the dynamics of predator–prey interaction may be analysed as an interesting topic. It has already been concluded that the prey refuge has a stabilizing effect on predator–prey dynamic and prey biomass can be controlled and saved from extinction due to predation [18–25]. Recently, there are few works on the dynamics of predator–prey model with prey refuge proportional to encounters between both prey and predator [26,27]. Consideration of prey refuge proportional to both biomass creates the system closer to reality as there are some biological systems, where prey refuge might be affected by both prey and predator. In this work, we have proposed a predator–prey model incorporating a prey refuge proportional to both the species (prey and predator), i.e., mxy amount of prey are free from predation. So, the predator attack only (x − mxy) amount of prey for their healthy life, where parameter m is the coefficient of refuge, x is prey biomass and y is predator biomass. 1 which assures the acceptable range 0 ≤ (1 − my) ≤ 1 of refuge for realistic It is assumed that (1 − my) ≥ 0, i.e., y ≤ m biological system. Also, predator is provided with additional food. Here, a Holling type II functional response is taken to represent the interaction between prey and predator in presence of additional food. The present work is organized as follows: in Section 2, a mathematical model is formulated based on the assumptions as discussed. In Section 3, some preliminaries such as positivity and uniform boundedness of solutions have been studied. In Section 4, existence and local stability analysis of each equilibria of the system in absence of delay (τ = 0) are analysed. Also, the condition for occurrence of Hopf-bifurcation is derived in this section. It is more realistic to assume that the reproduction of predator after predating the prey will not be an instantaneous biological process but mediated through some time lag required for gestation of the predator. So, in Section 5, we have analysed the dynamical behaviour of the delayed (gestation delay) model (τ ̸ = 0). The switching behaviour of τ is discussed mathematically through Hopfbifurcation. Also in Section 6, we have mentioned a systematic analysis on the dynamics of such predator–prey interaction with constant harvesting (‘‘constant-yield harvesting’’) in predators. Essentially, this work contains local bifurcations (saddle–node, Hopf) considering coefficient of refuge m (for saddle–node bifurcation) and quality of additional food α (for Hopf-bifurcation) as bifurcation parameters for the constant harvested predator–prey system. All analytical findings are verified numerically in Section 7 using MATLAB and MATCONT. It ends with conclusion in Section 8. 2. Model formulation The following assumptions are taken to formulate our basic mathematical model. 1. Predator is provided with additional food of constant biomass A which is distributed uniformly among the habitat. The number of encounters per predator with the additional food is proportional to the biomass of the additional food. The proportionality constant characterizes the ability of the predator to identify the additional food. 2. First, we consider a predator–prey model incorporating a refuge protecting nx (where n ∈ (0, 1]) of the prey: dx dt dy dt

( x) = rx 1 − −

a(1 − n)xy

k b + αηA + (1 − n)x ca {(1 − n)x + ηA} y = − dy b + αηA + (1 − n)x

(1)

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

3

with initial conditions: x(0) > 0, y(0) > 0, where x, y denote the prey and predator biomass respectively at time t and r, k, a, b, c, d are all positive constants. Here r represents intrinsic growth rate of prey; k is the carrying capacity of prey; d is the death rate of predator; a the per capita predator consumption rate; b is the half saturation constant in absence of additional food and refuge; c is the conversion factor denoting the number of newly born predators for each captured prey. The term a(1−n)xy represents the functional response (Holling type II) for the predator, where the term ηA denotes the b+αηA+(1−n)x effectual additional food level and α represents the quality of additional food [28,29]. It is noted that if n = 1, i.e., total prey population is refuged, then predator biomass is grown up in presence of additional food only. All the parameters contained in the model system (1) are positive and we assume that ca α> (2) d First, we will show that system (1) is uniformly bounded. It is evident that every solution of system (1) exists (unique) and positive for all t ≥ 0 (see Theorem 3.1). y Let us take W = x + . Then c dW dt

=

dx

1 dy

+ c dt ( x) a(1 − n)xy a {(1 − n)x + ηA} y dy = rx 1 − − + − k b + αηA + (1 − n)x b + αηA + (1 − n)x c ( dy ay x) − + ≤ rx 1 − k c α dt

= rx −

rx2

−

dy

k rx2

c

+

ay

α (

ca ) y − rx − d − k [α c ] ( ca ) y since lim sup x(t) ≤ k ≤ rk − rx − d − α c t →+∞ { } ca Let ξ = min r , d − α > 0 by using the condition (2).

= 2rx −

Applying the theorem of differential inequality [30], we obtain: 0 < W (x(t), y(t)) <

) rk ( 1 − e−ξ t + e−ξ t W (x(0), y(0)) . ξ

∴ 0 < W (x(t), y(t)) <

rk

ξ

, as t → ∞.

Thus, all solutions of system (1) enter into the region: B = (1) is uniformly bounded. So, lim supt →+∞ y(t) ≤

crk

ξ

{

(x, y) ∈ R2+ : 0 < W (x(t), y(t)) <

rk

ξ

}

. Hence system

. Therefore, both the prey and predator species are bounded

1 above. As a result, there exists a positive constant m such that y ≤ m . Let us assume that for prey: the admissible ranges of proportions of refuge and available prey for attack are respectively 0 ≤ my ≤ 1 and 0 ≤ (1 − my) ≤ 1.

Under these assumptions, the basic mathematical model can be formulated as: dx dt dy dt

( x) = rx 1 − −

a(1 − my)xy

k b + αηA + (1 − my)x ca {(1 − my)x + ηA} y = − dy b + αηA + (1 − my)x

(3)

with initial conditions: x(0) > 0, y(0) > 0.

(4)

3. Positivity and uniform boundedness Theorem 3.1.

Every solution of system (3) with initial conditions (4) exists (unique) and x(t) > 0, y(t) > 0 for all t ≥ 0.

Proof. Since the right hand side of system (3) is completely continuous and locally Lipschitzian in the domain R2+ , the solution of (x(t), y(t)) of (3) with initial conditions (4) exists and unique on [0, ξ ), where 0 < ξ ≤ +∞ [31].

4

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

From system (3) with initial conditions (4) , we have

[∫ t { ( r

x(t) = x(0) exp

[∫0 t { y(t) = y(0) exp

1−

x(θ )

a (1 − my(θ )) y(θ )

) −

k

b + αηA + (1 − my(θ ))x(θ )

ca {(1 − my(θ )) x(θ ) + ηA} b + αηA + (1 − my(θ ))x(θ )

0

}

}

dθ

]

>0

]

− d dθ > 0 □

This shows that system (3) is positively invariant for all t ≥ 0. Theorem 3.2. All solutions of system (3) which start in R2+ are uniformly bounded. y

Proof. Let us take W = x + dW

dx

=

c

. Then

1 dy

+ dt dt c dt ( x) a(1 − my)xy a {(1 − my)x + ηA} y dy = rx 1 − − + − k b + αηA + (1 − my)x b + αηA + (1 − my)x c ( ay dy x) + − ≤ rx 1 − k α c = rx −

rx2

+

k rx2

ay

α

−

dy c

ca ) y − rx − d − k [α c ] ( ca ) y ≤ rk − rx − d − since lim sup x(t) ≤ k α c t →+∞ { } Let ξ = min r , d − ca > 0 using (2). α

= 2rx −

(

Applying the theorem of differential inequality [30], we get: 0 < W (x(t), y(t)) <

) rk ( 1 − e−ξ t + e−ξ t W (x(0), y(0)) . ξ

∴ 0 < W (x(t), y(t)) <

rk

ξ

, as t → ∞.

Thus, { all solutions of system (3) enter }into the region:

B = (x, y) ∈ R2+ : 0 < W (x(t), y(t)) <

rk

. Hence system (3) is uniformly bounded. □

ξ

Remark. The bounded regions of the solutions of both the systems (1) and (3) are identical, i.e., there is no change of bounded regions due to ‘‘nx’’ or ‘‘mxy’’ amount of prey refuge. So, the formulation of model system (3) with (4) is meaningful. 4. Equilibrium points and their stability 4.1. Equilibria The trivial and axial equilibrium points E0 (0, 0) and E1 (k, 0) of system (3) always exist. The interior equilibrium point E ∗ (x∗ , y∗ ) can be obtained by solving the following equations:

(

r 1−

x) k

a(1 − my)y =0 b + αηA + (1 − my)x ca {(1 − my)x + ηA} −d = 0 b + αηA + (1 − my)x

−

Solving these two equations, we get a cubic equation in x.

χ1 x3 + χ2 x2 + χ3 x + χ4 = 0,

(5) (b+αηA)rm(ca−d)2

where p = db + ηA (α d − ca) > 0, χ1 = − − < 0 ( if ca > d), χ2 = rmp(ca − d) + (b + αηA)rm(ca − k k d)2 > 0, χ3 = −pca2 + pad < 0 and χ4 = {abd + aηA(α d − ca)} p > 0 using (2) The Eq. (5) has at least one positive root x∗ (say). mr(ca−d)p

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

5

Using the value of x∗ , we get the value of y∗ as: y∗ =

1 x∗ (ca − d) − db − ηA(α d − ca) x∗ (ca − d)

m

,

provided x∗ (ca − d) > db + ηA(α d − ca), m ̸ = 0 and ca > d. 4.2. Local stability analysis Now we discuss the local stability of the equilibrium points by computing corresponding Jacobian matrices. The Jacobian matrix J0 at E0 (0, 0) is given by:

⎡

r

0

J0 = ⎣ 0

caηA

⎤

b + αηA

⎦. −d caηA b+αηA

The eigen-values corresponding to J0 are r (> 0) and real number so E0 is unstable. The Jacobian matrix J1 at E1 (k, 0) is expressed as:

⎡ ⎢−r J1 = ⎢ ⎣

⎤

ak

−

− d (< 0 using (2)). Since one eigen-value is strictly positive

b + αηA + k ⎥ ⎥. ⎦ caηA −d b + αηA + k

0

caηA

The eigen-values corresponding to the Jacobian matrix J1 are −r (< 0) and b +αηA+k − d (< 0 using (2)). Since both 1 eigen-values are real and negative, so the axial equilibrium point E1 is Locally asymptotically stable (LAS). The Jacobian matrix J ∗ corresponding to interior equilibrium point E ∗ (x∗ , y∗ ) is given by: J∗ =

]

[

a11 a21

a12 , a22

{

where a11 = x∗ − kr + a21 = y∗ ∗

a22 = y

(

ca 1−my∗

{

)

}

a(1−my∗ )2 y∗ , a12 x∗ (b+αηA+(1−my∗ )x∗ )2 ∗ ∗ ∗ ca 1−my x +ηA 1−my

−

∗ )x∗

{ b+αηA+(1−my ∗

−camx b+αηA+(1−my∗ )x∗

+

)

{(

=

{(

)

}

}

−a(1−2my∗ ) b+αηA+(1−my∗ )x∗

} ) and

}(

(b+αηA+(1−my∗ )x∗ )2 camx∗ 1−my∗ x∗ +ηA (b+αηA+(1−my∗ )x∗ )2

{

−

(

)

amx∗ y∗ 1−my∗ (b+αηA+(1−my∗ )x∗ )2

}

,

.

The characteristic equation corresponding to J ∗ is given by

λ2 − A1 λ + B1 = 0,

(6)

where A1 = (a11 + a22 ) and B1 = a11 a22 − a12 a21 . Case-1: If A1 > 0 and B1 > 0, then the roots of Eq. (6) are either both real with positive signs or complex conjugate with positive real part. Therefore, E ∗ is either an unstable node or an unstable spiral . Case-2: If A1 < 0 and B1 > 0, then the roots of Eq. (6) are either both real with negative signs or complex conjugate with negative real part. So, E ∗ is either a stable node or a stable spiral. Case-3: If A1 > 0 or < 0 and B1 < 0, then the roots of Eq. (6) are both real with opposite sign to each other. Therefore, E ∗ is a saddle–node. Case-4: If A1 = 0 and B1 > 0, then the roots of Eq. (6) are purely complex conjugate to each other. So, E ∗ is a centre. Case-5: If A1 = 0 and B1 < 0, then the roots of Eq. (6) are both real numbers with same magnitude but opposite sign. Therefore, E ∗ is a saddle–node. 4.3. Global stability analysis of the interior equilibrium point Theorem 4.1. The interior equilibrium point E ∗ (x∗ , y∗ ) is globally asymptotically stable in the region:

{ R1 =

(x, y) :

y y∗

>

x x∗

} { } y x > 1 ∪ (x, y) : ∗ < ∗ < 1 y

x

Proof. Let us consider the following positive definite Lyapunov function in R about the interior equilibrium point E ∗ : V = x − x∗ − x∗ ln

x x∗

+

1 c

(

y − y∗ − y∗ ln

y y∗

)

.

(7)

6

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Differentiating (7) with respect to time t we get:

) x˙ 1 ( ) y˙ ( = x − x∗ + y − y∗ dt x c y { ( { } } ) ) ( ) ca {(1 − my)x + ηA} x a(1 − my)y 1( ∗ ∗ r 1− = x−x − + y−y −d k b + αηA + (1 − my)x c b + αηA + (1 − my)x

dV

⇒

dV

∗

(

= x−x

{ ) rx∗

+

a (1 − my∗ ) y∗

−

a (1 − my) y

−

}

b + αηA + (1 − my) x } {( } {( a 1 − my) x + ηA a 1 − my∗ ) x∗ + ηA} − + y−y b + αηA + (1 − my) x b + αηA + (1 − my∗ ) x∗ { } ) ) ( 1 r ( 1 dV ∗ ∗ 2 = − x − x + y − y aη A − ⇒ dt k b + αηA + (1 − my) x b + αηA + (1 − my∗ ) x∗ { } ∗ ∗ ) ( 1 − my y 1 − my y ( ) ( ) +a x − x∗ − b + αηA + (1 − my∗ ) x∗ b + αηA + (1 − my) x } { ( ) 1 − my x (1 − my∗ ) x∗ ( ) ∗ +a y − y − b + αηA + (1 − my) x b + αηA + (1 − my∗ ) x∗ ∗ ) {( dV r ( a η A y − y x − x∗ ) + m (xy − x∗ y∗ )} ( ) 2 ⇒ = − x − x∗ − dt k (b + αηA + (1 − my∗ ) x∗ ) (b + αηA + (1 − my) x) ( ∗ ) ( ∗ ) a (1 − my) a (1 − my∗ ) + x y − xy∗ − x y − xy∗ ∗ ∗ b + αηA + (1 − my) x b + αηA + (1 − my ) x dt

) ∗

(

⇒

dV dt

=−

r (

k

b + αηA + (1 − my∗ ) x∗

rx k

{

x − x∗

)2

aηA (y − y∗ ) {(x − x∗ ) + m (xy − x∗ y∗ )}

−

(b + αηA + (1 − my∗ ) x∗ ) (b + αηA + (1 − my) x) a (x∗ y − xy∗ ) χ − , (b + αηA + (1 − my∗ ) x∗ ) (b + αηA + (1 − my) x) k

(8)

where χ = {m (y − y∗ ) (b + αηA) + (x −{x∗ ) + m (xy − x∗ y∗ ) } + m{(1 − my) (x∗ y − xy∗})}. Therefore,

dV dt

< 0 in the region R1 = (x, y) :

y y∗

>

x x∗

> 1 ∪ (x, y) :

y y∗

<

x x∗

< 1 . Hence all the solutions of system

(3) in region R1 ultimately approach to the equilibrium point (x , y ). Therefore, V is a Lyapunov function in R1 so the interior equilibrium point E ∗ is globally asymptotically stable in the region R1 . □ ∗

∗

4.4. Hopf-bifurcation around E ∗ Characteristic equation of system (3) at E ∗ (x∗ , y∗ ) is given by

λ2 + A1 (a)λ + B1 (a) = 0,

(9)

where A1 = −(a11 + a22 ) and B1 = (a11 a22 − a12 a21 ). In order to see the instability of system (3), let us consider a as bifurcation parameter. For this purpose let us first state the following theorem: Theorem 4.2 (Hopf-bifurcation Theorem [8]). If A1 (a) and B1 (a) are the smooth functions of a in an open interval about a∗ ∈ R such that the characteristic equation (9) has a pair of imaginary eigen-values λ = q1 (a) ± iq2 (a) with q1 (a) and q2 (a) ∈ R so that they become purely imaginary at a = a∗ and

dq1 da a=a∗

⏐ ⏐

̸= 0, then a Hopf-bifurcation occurs around E ∗ at a = a∗ (i.e, a stability changes of E ∗ (x∗ , y∗ ) accompanied by the creation of a limit cycle at a = a∗ ). Theorem 4.3. B1 (a) > 0.

The system (3) possesses a Hopf-bifurcation around E ∗ when a passes through a∗ provided A1 (a) = 0 and

Proof. At a = a∗ , the characteristic equation (9) can be written as

λ2 + B1 (a) = 0.

(10)

√

√

The roots of Eq. (10) are λ1 = i B1 and λ2 = −i B1 . Thus a pair of purely imaginary eigen-values has been occurred at E ∗ . Also A1 and B1 are the smooth functions of a. So, the roots of Eq. (10) are of the form λ1 = q1 (a) + iq2 (a) and λ2 = q1 (a) − iq2 (a) in a neighbourhood of a∗ , where qi (a) are real functions for i = 1, 2. Next, let us verify the transversality

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

7

condition: d da

⏐ ⏐

(Reλi (a))⏐⏐

̸= 0, i = 1, 2.

a=a∗

Putting λ(a) = q1 (a) + iq2 (a) into Eq. (6), we get :

(q1 (a) + iq2 (a))2 + A1 (a) (q1 (a) + iq2 (a)) + B1 (a) = 0.

(11)

Differentiating both sides with respect to a, we have: 2 (q1 (a) + iq2 (a)) (q˙1 (a) + iq˙2 (a)) + A1 (a) (q˙1 (a) + iq˙2 (a)) + A˙1 (a) (q1 (a) + iq2 (a)) + B˙1 (a) = 0. Comparing real and imaginary parts from both sides, we obtain: 2q1 q˙1 − 2q2 q˙2 + A1 q˙1 + A˙1 q1 + B˙1 = 0,

(12)

2q1 q˙2 + 2q2 q˙1 + A1 q˙2 + A˙1 q2 = 0.

(13)

The Eq. (12) can be rewritten as: q˙1 (2q1 + A1 ) + q˙2 (−2q2 ) + A˙1 q1 + B˙1 = 0

⇒ q˙1 X1 − q˙2 X2 + X3 = 0

(14)

and Eq. (13) can be rewritten as: q˙1 (2q2 ) + q˙2 (2q1 + A1 ) + A˙1 q2 =

0

⇒ q˙1 X2 + q˙2 X1 + X4 = 0 ,

(15)

where X1 = (2q1 + A1 ) ̸ = 2q1 , X2 = 2q2 , X3 = A˙1 q1 + B˙1 and X4 = A˙1 q2 . Multiplying (14) by X1 and (15) by X2 and, then adding, we get:

) X1 X3 + X2 X4 ) X12 + X22 q˙1 + X1 X3 + X2 X4 = 0 ⇒ q˙1 = − ( 2 X1 + X22

(

At a = a∗ ,

√ Case 1: q1 = 0, q2 =

√

B1 . Then X1 ̸ = 0, X2 = 2 B1 , X3 = B˙1 and X4 = A˙1

√

B1 . ∴ X1 X3 + X2 X4 ̸ = 0. √ √ √ ˙ case 2: q1 = 0, q2 = − B1 . Then X1 ̸ = 0, X2 = −2 B1 , X3 = B1 and X4 = −A˙1 B1 . ∴ X1 X3 + X2 X4 ̸ = 0. Hence Theorem 4.3 is proved using Theorem 4.2.

□

5. Effect of discrete time-delay Many biological processes (both natural and man-made) involve time delays. In 1925, Volterra first proposed a delayed predator–prey model to study fish population under harvesting [32,33]. Generally, delay differential equations have much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable through Hopf-bifurcation. Kuang [34] studied that some time is needed to digest food before further activities have to be made. Since some amounts of energy assimilate from prey into predator’s energy level after predation, the study of discrete time delay is more realistic in ecology. This conversion process is not instantaneous because many biological processes have to be performed for completion of this mechanism. The whole transformation process takes some time which is known as gestation delay. However, predator–prey models with discrete delays can undergo complex dynamical behaviour such as the existence of Hopf-bifurcation, Bogdanov–Takens bifurcation and even chaos (see, [35,36] etc.). In this section, we will study the effect of discrete time-delay on the model system (3). The corresponding delay differential equations can be written as: dx dt dy dt

( x) = rx 1 − − =

a(1 − my)xy

k b + αηA + (1 − my)x ca {(1 − my(t − τ ))x(t − τ ) + ηA} y(t − τ ) b + αηA + (1 − my(t − τ ))x(t − τ )

(16)

− dy

We represent by C the Banach space of continuous functions ϕ : [−τ , 0] → R2 with

∥ϕ∥ =

sup {|ϕ1 (ψ )|, |ϕ2 (ψ )|} −τ ≤ψ≤0

where ϕ = (ϕ1 , ϕ2 ). As usual, the initial conditions of this system are given as follows: x(ψ ) = ϕ1 (ψ ), y(ψ ) = ϕ2 (ψ ), ψ ∈ [−τ , 0]

8

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

The initial function ϕ = (ϕ1 , ϕ2 )T belongs to the Banach space C = C ([−τ , 0], R2 ) of continuous functions mapping the initial [−τ , 0] into R2 . In view of biological reasons, it is assumed that:

ϕi (ψ ) ≥ 0, ψ ∈ [−τ , 0], i = 1, 2. Using elementary theory of functional differential equations [31], it can be proved that there is a unique solution (x(t), y(t)) to this system. It is also assumed that ϕi (0) > 0 for i = 1, 2 for a biological feasibility. Now, let us linearize system (16) by using the transformations: X = x − x∗ and Y = y − y∗ . Then the linearized system is expressed as: dU1 dt

= CU1 (t) + DU1 (t − τ ),

(17)

where U1 = [X , Y ]T ,

[ C =

a11 0

a12 −d

]

and

[

0 D= a21

0

]

a22 + b22

ca (1−my∗ )x∗ +ηA

{

,

}

where b22 = b+αηA+(1−my∗ )x∗ and a11 , a12 , a21 and a22 all are defined earlier. The characteristic equation of system (17) is given by:

λ2 + C1 λ + C2 + e−λτ {D1 λ + D2 } = 0,

(18)

where C1 = −(a11 − d), C2 = −a11 d, D1 = − (a22 + b22 ) and D2 = a11 (a22 + b22 ) − a12 a21 . In presence of delay τ (̸ = 0), system (16) is LAS if Eq. (18) has no purely imaginary roots and it is asymptotically stable for τ = 0. But there exists τ = τ ∗ , where changes of stability occurs. We have also studied that E ∗ is LAS if the conditions A1 = a11 + a22 < 0 and B1 = a11 a22 − a21 a12 > 0 hold simultaneously for τ = 0. Now, we will examine if the real part of the roots of Eq. (18) increases to reach zero and eventually becomes positive as τ varies. Substituting λ = p1 + ip2 in Eq. (18), we get:

(p1 + ip2 )2 + C1 (p1 + ip2 ) + C2 + e−(p1 +ip2 )τ {D1 (p1 + ip2 ) + D2 } = 0 Comparing real and imaginary parts from both sides, we obtain: p21 − p22 + C1 p1 + C2 + (D1 p1 + D2 ) e−p1 τ cos(p2 τ ) + D1 p2 e−p1 τ sin(p2 τ ) = 0

(19)

2p1 p2 + C1 p2 + D1 p2 e−p1 τ cos(p2 τ ) − (D1 p1 + D2 ) e−p1 τ sin(p2 τ ) = 0

(20)

Now, let us observe whether Eq. (18) has purely imaginary roots or not. Putting p1 = 0 in both Eqs. (19) and (20), we have:

−p22 + C2 + D2 cos(p2 τ ) + D1 p2 sin(p2 τ ) = 0

(21)

C1 p2 + D1 p2 cos(p2 τ ) − D2 sin(p2 τ ) = 0

(22)

Eliminating τ from (21) and (22), we get: p42 + R1 p22 + R2 = 0, where R1 = −2C2 + C12 − D21 R2 = C22 − D22 Let p22 = σ , then this equation becomes: R ≡ σ 2 + R1 σ + R2 = 0

(23)

This is a quadratic equation in σ . It is clearly noticed that R(∞) = ∞. So, the equation has at least one positive real root if R(0) < 0, i.e., if C22 < D22 . √ If σ = σ+ be the positive root of Eq. (23), then p2 = ± σ+ . Here we describe a lemma which was proved by Ruan and Wei [37]:

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Lemma 5.1.

9

Consider the exponential polynomial:

) ( (0) (0) P λ, e−λτ1 , e−λτ2 , . . . , e−λτm = λn + p1 λn−1 + · · · + pn−1 λ + p(0) n [ ] (1) n−1 (1) (1) + p1 λ + · · · + pn−1 λ + pn e−λτ1 + .[. . ] (m) n−1 (m) + p(m) λ + · · · + p e−λτm , λ + p n 1 n−1 (i)

where τi ≥ 0 (i = 1, 2, . . . , m) and pj (i = 0, 1, . . . , m; j = 1, 2, . . . , n) are constants. As (τ1 , τ2 , . . . , τm ) vary, the sum of ( ) the orders of zero of P λ, e−λτ1 , e−λτ2 , . . . , e−λτm in the open half plane can change only if a zero appears on or crosses the imaginary axis. Now we discuss the occurrence of Hopf-bifurcation around the interior equilibrium point E ∗ considering τ as bifurcation parameter. Theorem 5.2. Suppose E ∗ exists and is LAS if A1 < 0 and B1 > 0 for system (3) when τ = 0. If C22 < D22 , then there exists τ = τ ∗ such that E ∗ of the delay system (16) is LAS when τ ∈ [0, τ ∗ ) and unstable when τ > τ ∗ , where cos−1

{

D2 {σ+ −C2 }−C1 D1 σ+

} 2nπ

D22 +D21 σ+

τ+(n) =

+√

√ σ+

σ+

, n = 0, 1, 2, . . .

(0)

and τ ∗ = τ+ (minimum value). In other words, the system (16) exhibits a supercritical Hopf-bifurcation around E ∗ for τ = τ ∗ provided XU > YV , where X , Y , U and V are described in the proof. (n)

Proof. If C22 < D22 , then Eq. (23) has one positive root, say σ+ . From Eqs. (21) and (22), we see that τ+ is a function of σ+ for n = 0, 1, 2, . . .; which is represented as: cos−1

{

D2 {σ+ −C2 }−C1 D1 σ+

} 2nπ

D22 +D21 σ+

τ+(n) =

+√

√ σ+

σ+

, n = 0, 1, 2, . . .

Therefore, the interior equilibrium point E ∗ is locally asymptotically stable if A1 < 0 and B1 > 0 for τ = 0. In this case E ∗ will remain stable for τ < τ ∗ (using Butler’s lemma). ] [ Now, let us verify the transversality condition: ddτ Reλ(τ ) τ =τ ∗ > 0. Differentiating (19) and (20) with respect to τ and putting p1 = 0 and τ = τ ∗ , we get:

[ X

d dτ

[ −Y

[Re {λ(τ )}]

]

[ +Y

]τ =τ

d

[Re {λ(τ )}]

dτ

dτ

∗

[ +X

τ =τ ∗

d d dτ

[Im {λ(τ )}] [Im {λ(τ )}]

] ]τ =τ

=U

(24)

= V,

(25)

∗

τ =τ ∗

where X = {C1 + D1 cos(p2 τ ) − D2 τ cos(p2 τ ) − D1 p2 τ sin(p2 τ )}τ =τ ∗ ,

Y = {−2p2 − D2 τ sin(p2 τ ) + D1 sin(p2 τ ) + D1 p2 τ cos(p2 τ )}τ =τ ∗ , U = D2 p2 sin(p2 τ ) − D1 p22 cos(p2 τ )

√ τ =τ ∗ , p2 = σ+

V = D2 p2 cos(p2 τ ) + D1 p22 sin(p2 τ )

√ τ =τ ∗ , p2 = σ+

{

}

{

√

p2 = σ+

√

p2 = σ+

}

Solving Eqs. (24) and (25):

[

d [Re {λ(τ )}] dτ

Now we have

]

[ = τ =τ ∗

[ d[Re{λ(τ )}] ] dτ

τ =τ ∗

XU − YV X2 + Y 2

]

.

> 0, provided XU > YV . Hence the transversality condition is verified and a Hopf-bifurcation

occurs when τ passes through the critical value τ ∗ .

□

6. Mathematical formulation with constant harvesting In ecology, dynamical behaviour of predator–prey interactions with ‘‘constant-yield harvesting’’ (harvested biomass is independent of the population size) and ‘‘constant-effort harvesting’’ (a constant proportion of the population is harvested)

10

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

are discussed briefly by several researchers [38,39]. Most of researchers have observed the dynamics of concerned models with the harvesting of either prey or predator or both the species as the outcomes of such analysis which are more relevant for the management of renewable resources [40,41]. On the dynamics of predator–prey model, effect of harvesting is not only an interesting concept from theoretical point of view but also significant for the management and sustainability of commercially/economically important resources [42,43]. Some systematic studies of harvested predator–prey interactions assure the chance of extinction of one or both the species due to uncontrolled harvesting. Moreover, it determines the acceptable range of harvesting to assure the long term survival of the renewable resource. Now, we wish to investigate the effect of harvesting the predator at a constant rate in presence of additional food to the predator. Considering constant yield harvesting (h) into system (3), we obtain the following system:

( x) = rx 1 − −

dx

a(1 − my)xy

k b + αηA + (1 − my)x ca {(1 − my)x + ηA} y = − dy − h, b + αηA + (1 − my)x

dt dy dt

(26)

with initial conditions: x(0) > 0, y(0) > 0. Remark. From the second equation of system (26), we get: t

∫

exp {−P(ω)} dω,

y(t) = y(0) exp {P(t) − P(0)} − h exp {P(t)} 0

where P(ω) =

ca{(1−my(ω))x(ω)+ηA} b+αηA+(1−my(ω))x(ω)

− d.

Therefore, T

∫

exp {−P(ω)} dω.

y(T ) = 0 ⇒ y(0) = h exp {P(0)}

(27)

0

Therefore, we must stop to harvest y before time T (due to unavailability of predators in the system), where T is given by Eq. (27). Hence, without loss of generality we assume that x(0) > 0 and y(0) > h for system (26). Now we interest to study the dynamics of model (26) and investigate the role of constant harvesting on the dynamics of considered predator–prey system. All the parameters are positive and defined earlier. 6.1. Existence of equilibria The equilibrium points can be obtained by solving the following equations:

(

x)

a(1 − my)xy

=0 k b + αηA + (1 − my)x ca {(1 − my)x + ηA} y − dy − h = 0 b + αηA + (1 − my)x

rx 1 −

−

They are Ey (0, y¯ ) (where y¯ = caηA−d(b+αηA) exists only if caηA > d(b + αηA)), E1∗ (x∗1 , y∗1 ) and E2∗ (x∗2 , y∗2 ), where both x∗i , yi (for i = 1, 2) are real, positive and can be found by using MATHEMATICA. h(b+αηA)

∗

6.2. Local stability The Jacobian matrix Jy at Ey (0, y¯ ) is expressed as:

⎡

a(1 − my¯ )y¯

r− ⎢ b + αηA Jy = ⎢ ⎣ ca(1 − my¯ )y¯ (b + αηA − ηA) (b + αηA)2

⎤ 0 caηA b + αηA

−d

⎥ ⎥. ⎦

Two eigen-values of system (26) corresponding to the Jacobian matrix Jy are r − eigen-value

caηA b+αηA

a(1−my¯ )y¯ b+αηA

caηA b+αηA

− d. The second − d is positive whenever Ey is feasible. The stability nature of Ey can be described as follows:

Lemma 6.1. Ey is an unstable node if r −

a(1−my¯ )y¯ b+αηA

> 0 and is a saddle-point if r −

a(1−my¯ )y¯ b+αηA

and

< 0.

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

11

The Jacobian matrix Ji∗ corresponding to interior equilibrium point Ei∗ (x∗i , y∗i ) is given by:

[

]

b Ji = 11 b21 ∗

b12 , b22

{

where b11 = x∗i − kr + b21 b22

a(1−my∗ )2 y∗ i i

}

, b12 = x∗i

{

−a(1−2my∗i ) b+αηA+(1−my∗ x∗ i ) i

) } {(1−my∗i )x∗i +ηA}(1−my∗i ) = yi b+αηA+ 1−my∗ x∗ − , ( (b+αηA+(1−my∗i )x∗i )2 } i ) i { ∗ ∗ ∗ ∗ ∗ ∗ ca[(1−2my )xi +ηA] camxi yi {(1−myi )xi +ηA} = b+αηA+ 1i−my for i = 1, 2. ∗ x∗ − d + ( (b+αηA+(1−my∗i )x∗i )2 i ) i ∗

{

(

ca 1−my∗ i

)

(

2 b+αηA+(1−my∗ )x∗ i i

−

} ( ) , 2 (b+αηA+(1−my∗i )x∗i ) amx∗ y∗ 1−my∗ i i i

ca

The trace and the determinant of the Jacobian matrix Ji∗ for system (26) calculated at the interior equilibrium points are given by: Tr(Ji∗ ) = b11 + b22 Det(Ji∗ ) = b11 b22 − b12 b21 , i = 1, 2. The stability properties for interior equilibrium points Ei∗ are stated below: Lemma 6.2. The interior equilibrium points Ei∗ are locally asymptotically stable if and only if Tr(Ji∗ ) < 0 and Det(Ji∗ ) > 0 (for i = 1, 2) by Routh Hurwitz criteria. 6.3. Local bifurcation In this sub-section, we analyse the local bifurcation results of system (26). Here, conditions for saddle–node and Hopfbifurcations are derived which are of co-dimension one. The detailed calculations for the derivation of normal forms (near the bifurcation threshold) are described as follows. 6.3.1. Saddle–node bifurcation The following theorem states the restrictions for occurrence of a saddle–node bifurcation for system (26). Theorem 6.3. System (26) exhibits a saddle–node bifurcation when the system parameters satisfy the restriction Det(J1∗ ) = 0 along with the condition Tr(J1∗ ) < 0. Here the saddle–node bifurcation threshold with respect to m (coefficient of refuge) as the bifurcation parameter is presented as m = m[SN ] . Proof. To verify the transversality condition for the occurrence of saddle–node bifurcation, let us apply Sotomayor’s theorem [44] at m = m[SN ] . It can be easily noticed that Det(J1∗ ) = 0 and Tr(J1∗ ) < 0 at m = m[SN ] which indicates that the Jacobian matrix J1∗ admits a zero eigenvalue denoted by λ1 . If V and W represent eigenvectors corresponding to the zero eigen-value for the matrices J1∗ and (J1∗ )T respectively, then they are given by V = (v1 , v2 )T and W = (w1 , w2 )T , where [SN ] [SN ] v1 = − cc12 , v2 = 1, w1 = 1, w2 = − cc11 , c11 = b[11SN ] , c12 = b[12SN ] , c21 = b21 and c22 = b22 (b11 , b12 , b21 and b22 all are 11 21 defined earlier). Clearly, V and W satisfy the transversality conditions:

( ) ∆1 = W T .Fm x∗1 , y∗1 ; m[SN ] ̸= 0, [ ( ) ] ∆2 = W T . D2 F x∗1 , y∗1 ; m[SN ] (V , V ) ̸= 0, where F = (F1 , F2 )T and F1 , F2 are given by

(

F1 = rx 1 −

F2 =

x) k

−

a(1 − my)xy b + αηA + (1 − my)x

ca {(1 − my)x + ηA} y b + αηA + (1 − my)x

and

− dy − h. □

6.3.2. Hopf-bifurcation around E1∗ Characteristic equation of system (26) at E1∗ (x∗1 , y∗1 ) is given by

λ2 + K (α )λ + L(α ) = 0,

(28)

where K = −Tr(J1∗ ) and L = Det(J1∗ ) are defined earlier. In order to see the instability of system (26), let us consider α as bifurcation parameter. For this purpose let us first state the following theorem:

12

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Theorem 6.4 (Hopf-bifurcation Theorem [8]). If K and L are the smooth functions of α in an open interval about α ∗ ∈ R such that the characteristic equation (28) has a pair of ⏐ imaginary eigen-values λ = ξ (α ) ± iκ (α ) with ξ (α ) and κ (α ) ∈ R so that dξ they become purely imaginary at α = α ∗ and dα ⏐α=α ∗ ̸ = 0, then a Hopf-bifurcation occur around E1∗ at α = α ∗ (i.e, a stability changes of E1∗ (x∗1 , y∗1 ) accompanied by the creation of a limit cycle at α = α ∗ ). Theorem 6.5. System (26) exhibits a Hopf-bifurcation around E1∗ when α passes through α ∗ provided K (α ) = 0 and L(α ) > 0. Proof. At α = α ∗ , the characteristic equation (28) can be written as

λ2 + L(α ) = 0.

(29)

√

√

The roots of Eq. (28) are λ1 = i L and λ2 = −i L. Thus a pair of purely imaginary eigen-values has been occurred at E1∗ . Also K and L are the smooth functions of α . So, the roots of Eq. (28) are of the form λ1 = ξ (α ) + iκ (α ) and λ2 = ξ (α ) − iκ (α ) in a neighbourhood of α ∗ , where ξ (α ) and κ (α ) are real functions. Next, let us verify the transversality condition:

⏐ ⏐ (Reλi (a))⏐⏐ dα d

α=α ∗

̸= 0, i = 1, 2.

Putting λ(a) = ξ (α ) + iκ (α ) into Eq. (28), we get :

(ξ (α ) + iκ (α ))2 + K (α ) (ξ (α ) + iκ (α )) + L(α ) = 0.

(30)

Differentiating both sides with respect to α , we have: 2 (ξ (α ) + iκ (α )) ξ˙ (α ) + iκ˙ (α ) + K (α ) ξ˙ (α ) + iκ˙ (α )

(

)

(

)

˙ α ) = 0. +K˙ (α ) (ξ (α ) + iκ (α )) + L( Comparing real and imaginary parts from both sides, we obtain: 2ξ ξ˙ − 2κ κ˙ + K ξ˙ + K˙ ξ + L˙ = 0,

(31)

2ξ κ˙ + 2κ ξ˙ + K κ˙ + K˙ κ = 0.

(32)

The Eq. (31) can be rewritten as:

ξ˙ (2ξ + K ) + κ˙ (−2κ ) + K˙ ξ + L˙ = 0 ⇒ ξ˙ X1 − κ˙ X2 + X3 = 0

(33)

and Eq. (32) can be rewritten as:

ξ˙ (2κ ) + κ˙ (2ξ + K ) + K˙ κ = 0 ⇒ ξ˙ X2 + κ˙ X1 + X4 = 0 ,

(34)

where X1 = (2ξ + K ) ̸ = 2ξ , X2 = 2κ , X3 = K˙ ξ + L˙ and X4 = K˙ κ . Multiplying (33) by X1 and (34) by X2 and, then adding, we get:

) X1 X3 + X2 X4 ) X12 + X22 ξ˙ + X1 X3 + X2 X4 = 0 ⇒ ξ˙ = − ( 2 X1 + X22

(

At α = α ∗ ,

√

√

Case 1: ξ = 0, κ = L. Then X1 ̸ = 0, X2 = 2 L, X3 = L˙ and X4 = K˙ ∴ X1 X3 + X2 X4 ̸ = 0.

√

√ L.

√

case 2: ξ = 0, κ = − L. Then X1 ̸ = 0, X2 = −2 L, X3 = L˙ and X4 = −K˙ ∴ X1 X3 + X2 X4 ̸ = 0. Hence Theorem 6.5 is proved using Theorem 6.4. □

√ L.

Now, we will discuss the nature of Hopf-bifurcation. For this purpose, we introduce the perturbation x = x1 + xα1 and ∗ y = y1 + yα1 in system (26). Expanding in Taylor series, we obtain: ∗

x˙1 = β10 x1 + β01 y1 + β20 x21 + β11 x1 y1 + β30 x31 + β21 x21 y1 + · · · , y˙1 = γ10 x1 + γ01 y1 + γ

2 20 x1

+ γ11 x1 y1 + γ

3 30 x1

+γ

2 21 x1 y1

(35)

+ ···,

(36)

where β10 , β01 , γ10 and γ01 are the entries of the Jacobian matrix J1 evaluated at the equilibrium point E1 at α = α ∗ . Therefore β10 + γ01 = 0 and β10 γ01 − β10 γ10 > 0. Other coefficients are determined by: ∗

∗

β20 =

∂ 2 F1 2 ∂ x2 1

[

] α=α ∗

, β11 =

[

∂ 2 F1 ∂ x∂ y

] α=α ∗

, β30 =

∂ 3 F1 6 ∂ x3

1

[

] α=α ∗

, β21 =

∂ 3 F1 2 ∂ x2 ∂ y 1

[

] α=α ∗

,

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

γ20 =

∂ 2 F2 2 ∂ x2 1

[

] α=α ∗

, γ11 =

[

∂ 2 F2 ∂ x∂ y

] α=α ∗

, γ30 =

∂ 3 F2 6 ∂ x3

1

[

] α=α ∗

and γ21 =

[

∂ 3 F2 ∂ x2 ∂ y

13

] α=α ∗

,

where

(

F1 = rx 1 −

x) k

−

a(1 − my)xy b + αηA + (1 − my)x

,

ca {(1 − my)x + ηA} y

− dy − h. b + αηA + (1 − my)x The first Lyapunov number [44] to determine the nature (stable or unstable) of limit cycle arising through Hopfbifurcation is given by: F2 =

σ =

−3π 2β01 (β10 γ01 − β01 γ10 )

{[ 3 2

) ( 2 2 2 + β20 γ11 − 2β10 β01 β20 + β10 β01 γ11 β10 γ10 β11

) ) ] ( 2 ( 2 2 + β01 γ10 β11 β20 − β10 −β01 (2β20 γ20 + γ11 γ20 ) − β01 γ10 − 2β10 [−3β01 β30 + 2β10 (β21 + γ12 ) + (γ10 β12 − β01 γ21 )]} If σ < 0, the equilibrium point E1∗ is destabilized through a supercritical Hopf-bifurcation and if σ > 0, then the Hopf-bifurcation is subcritical. We have plotted the Hopf-bifurcation diagram in MATLAB and noticed that the nature of the Hopf-bifurcation is subcritical. To verify our claim we consider a numerical example. Let us choose a set of parametric values {r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.03, η = 0.9, A = 1.9, h = 0.14} and the Hopf-bifurcation threshold for α comes as α ∗ = 1.389. For these choice of parametric values one can calculate the first Lyapunov number as σ = 0.0047 > 0. Hence the Hopf-bifurcation is subcritical. 7. Numerical simulation Numerical simulations of system (3) are performed to verify our analytical findings and to find some other features. We have selected the set of parameters as: {r = 1.5, k = 5.5, d = 10, a = 2.0, b = 1.5, c = 0.8, m = 0.5, α = 0.2, η = 0.5, A = 2}. For these parameter values, the axial equilibrium point E1 (k, 0) ≡ (5.5, 0) is locally asymptotically stable (LAS) and stability nature of E1 with respect to time t is depicted in Fig. 1. Also, Fig. 2 shows the global attractor of E1 . As per the following set: {r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.01, α = 1.7, η = 0.9, A = 1.9}, E ∗ (x∗ , y∗ ) ≡ E ∗ (1.8837, 5.2397) which is LAS (imaginary eigen-values with negative real parts). The stability behaviour with respect to time t and stable phase trajectory of E ∗ are shown in Fig. 3. Fig. 4 depicts that different phase trajectories start from different initial points (0.5, 1.2), (1, 2.5), (3, 8.2) and (3.5, 10) but ultimately all trajectories go to a fixed point E ∗ (1.8837, 5.2397). So, E ∗ (1.8837, 5.2397) is globally asymptotically stable in the region R1 (defined earlier). It has been seen that if we take a = 1.4 and all the parameters are same as in Fig. 3, then E ∗ is unstable (imaginary eigen-values with positive real parts). The oscillating behaviour with respect to time t and phase diagram of limit cycle (isolated closed trajectory) for a = 1.4 are demonstrated in Fig. 5. So, the parameter ‘‘a’’ (consumption rate of predator) is of great importance due to its switching behaviour, i.e., when a(= 1.2) < a∗ (= 1.27), E ∗ is stable and E ∗ is unstable for a(= 1.4) > a∗ (= 1.27). Since the vector fields for a < a∗ and a > a∗ are qualitatively different, a Hopf-bifurcation is occurred at a = a∗ . Bifurcation diagrams (Fig. 6) shows that system (3) has coexistence stable dynamics for 1 ≤ a < 1.27 and oscillatory coexistence behaviour for a > 1.27. Moreover, the bifurcation diagrams of system (3) with respect to the consumption rate a (0.85 ≤ a ≤ 1.2) are depicted in Fig. 7 in absence of additional food (α = 0.0, η = 0.0, A = 0.0) and other parameters are r = 1.5, k = 6.5, d = 0.6, b = 0.9, c = 0.8, m = 0.01. It is observed from Fig. 7 that system (3) has coexistence stable dynamics for 0.85 ≤ a < 0.997 and oscillatory coexistence nature for a > 0.997. In this case, the dynamical behaviour switches at the critical point a∗ = 0.997. Also, from Figs. 6 and 7, one can easily see that the interior equilibrium point changes its stability nature for comparatively high level of consumption rate in presence of additional food. Also, it has been observed from Fig. 8 that when m lies in (0, 0.13], i.e., when coefficient of refuge is very small, then the growth of predator is greater than that of prey but when m ∈ (0.13, 1.5), then gradual increase of the coefficient of refuge indicates more protection of prey species, as a result, the growth of prey species increases and that of predator decreases. From Figs. 8 and 9, it has been observed that the growth of predator in presence of additional food is large enough to that of without additional food for very small m. Therefore, provision of additional food gives more realistic dynamics in ecology. On the other hand, it has to be noted that we always preserve the restriction my∗ ≤ 1 and α > ca d throughout the numerical computations. Now, we perform numerical simulations to confirm the analytical results of delayed system (16) (τ ̸ = 0). As per the following set: {r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.01, α = 1.7, η = 0.9, A = 1.9}, Eq. (18) has complex roots with negative real parts for τ = 4.0 and the roots of that equation will be purely imaginary when τ = τ ∗ = 4.2434. It is to be noted that the real parts of the roots of (18) is positive when τ = 4.5 > τ ∗ . As a result, it is

14

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Fig. 1. The stability behaviour of the axial equilibrium point E1 (k, 0) with respect to time t for the set of parameter values r = 1.5, k = 5.5, d = 10, a = 2.0, b = 1.5, c = 0.8, m = 0.5, α = 0.2, η = 0.5, A = 2.

Fig. 2. Global attractor of E1 (5.5, 0).

Fig. 3. E ∗ (1.8837, 5.2397): (a) stability nature with respect to time t and (b) stable phase portrait for the set of parameter values r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.01, α = 1.7, η = 0.9, A = 1.9.

Fig. 4. Phase trajectory of E ∗ (1.8837, 5.2397) depicts that E ∗ is globally asymptotically stable.

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

15

Fig. 5. E ∗ (1.8837, 5.2397): (a) Oscillation behaviour with respect to time t and (b) phase diagram of limit cycle for ‘‘a = 1.4’’ and other parameters are fixed as in Fig. 3.

Fig. 6. Bifurcation diagram of (a) prey (x) and (b) predator (y) as the consumption rate a crosses through its critical value ‘‘a = a∗ = 1.27’’ and other parameters are fixed as in Fig. 3.

concluded that E ∗ of system (16) is LAS when τ < τ ∗ and unstable when τ > τ ∗ . Figs. 11 and 12 represent (a) stability and oscillation behaviour with respect to time t and (b) stable phase trajectories and phase diagram of limit cycle respectively. So, system (16) exhibits a Hopf-bifurcation around E ∗ considering τ as bifurcation parameter (τ ∗ = 4.2434), i.e., a periodic solution occurs around E ∗ . As per the following set: {r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.03, α = 1.7, η = 0.9, A = 1.9, h = 0.14}, system (26) has two interior equilibrium points E1∗ (x∗1 , y∗1 ) ≡ (2.6488, 5.3163) and E2∗ (x∗2 , y∗2 ) ≡ (6.0163, 0.9239). There is no axial equilibrium point in this case. It has been observed that (2.6488, 5.3163) is a stable spiral and (6.0163, 0.9239) is a saddle point. The stable behaviour of the interior equilibrium point (2.6488, 5.3163) is presented in Fig. 13. If the quality of additional food is chosen as α = 1.2 and all other parameters are fixed as in Fig. 13, interior equilibrium point (2.6488, 5.3163) is unstable which are shown in Figs. 14 and 15. From here, it is concluded that system (26) undergoes a Hopf-bifurcation around (2.6488, 5.3163) considering α = (α ∗ = 1.389) as bifurcation parameter. The Hopf-bifurcation diagrams are depicted in Fig. 16. For different set of parameter values: {r = 4.5, k = 9.5, d = 0.4, a = 1.3, b = 0.9, c = 0.8, m = 0.07, α = 0.7, η = 0.4, A = 1.2, h = 0.14} of system (26), we get two interior equilibrium points E1∗ (x∗1 , y∗1 ) ≡ (9.4614, 14.2681) and E2∗ (x∗2 , y∗2 ) ≡ (9.4364, 0.2477) and one axial equilibrium point Ey (0, y¯ ) ≡ (0, 36.0500) (correct up to 4 decimal places). The eigen-values corresponding to (9.4614, 14.2681) are −4.9483 and −4.2967 which are both negative. So, (9.4614, 14.2681) is locally asymptotically stable (LAS) which is depicted in Fig. 17. For other interior equilibrium point (9.4364, 0.2477), eigen-values are −4.4429 and 0.5636. Since one eigen-value is positive and other is negative, (9.4364, 0.2477) is a saddle point. The eigen-values corresponding to the axial equilibrium point (0, 36.0500) are 62.2660 and 0.0039 which are positive. So, the axial equilibrium point (0, 36.0500) is an unstable node. But as per the following set: {r = 0.02, k = 1.5, d = 0.1, a = 1.5, b = 0.7, c = 0.8, m = 0.01, α = 0.3, η = 0.6, A = 1.9, h = 0.1}, Ey (0, y¯ ) ≡ (0, 0.0824) which is a saddle point because the eigen-values are −0.0986 < 0 and 1.2129 > 0. Also, for {r = 4.5, k = 9.5, d = 0.4, a = 1.3, ∗ (b = 0.9, c = 0.8, m) = 0.07, α = 0.7, η = 0.4, A = 1.2, h = 0.14}: a saddle–node bifurcation occurs around E1 taking m m[SN ] = 0.014678 as bifurcation parameter and bifurcation diagram is presented in Fig. 18. 8. Concluding remarks In this work, we have investigated the dynamical behaviour of a predator–prey model incorporating prey refuge proportional to both species. Also, we have studied the impact of additional food provided to the predator. Here, the

16

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Fig. 7. Bifurcation diagram of (a) prey (x) and (b) predator (y) as the consumption rate a crosses through its critical value ‘‘a = a∗ = 0.997’’ and other parameters are r = 1.5, k = 6.5, d = 0.6, b = 0.9, c = 0.8, m = 0.01, α = 0.0, η = 0.0, A = 0.0.

Fig. 8. Growth of predator–prey population in presence of additional food when the coefficient of refuge m increases and other parameters are fixed as in Fig. 3.

Fig. 9. Growth of population with respect to m in absence of additional food.

interaction between the predator and the prey is governed by Holling type II response function in presence of additional food. The innovation of our work is to analyse such kind of refuge function with the effect of additional food for a more realistic phenomenon in ecosystem. Existence criteria of biological meaningful equilibrium points have been obtained and their stability analysis have also been carried out. The effect of the parameter m (coefficient of refuge) over the growth of both the species has been plotted in Fig. 8. Also, the parameter a (consumption rate of predator) has an effective role in presence of additional food because system (3) exhibits a Hopf-bifurcation around E ∗ taking a as a bifurcation parameter. On the other hand, some time is needed to convert the consumed prey and additional food into predator’s growth which is termed as gestation period (τ ). So, in Section 5, we have studied the dynamical behaviour of delayed (gestation delay) predator–prey model. It has been noticed that the delay parameter τ plays a crucial role because there exists a

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

17

Fig. 10. Effect of m over the growth of population biomass for both the systems (3) and (1) corresponding to the set of parameters: {r = 1.2, k = 6.5, d = 0.1, a = 1.1, b = 0.4, c = 0.15, α = 0.6, η = 0.2, A = 1.2} and m ∈ [0.183, 0.5]. In system (1): n = m.

Fig. 11. E ∗ (1.8837, 5.2397): (a) stability criteria with respect to time t and (b) stable phase portrait when ‘‘τ = 4.0 < τ ∗ = 4.2434’’ and all parameters are fixed as in Fig. 3.

Fig. 12. E ∗ (1.8837, 5.2397): (a) Oscillation behaviour with respect to time t and (b) phase diagram with limit cycle when ‘‘τ = 4.5 > τ ∗ = 4.2434’’ and all parameters are fixed as in Fig. 3.

Fig. 13. Stable dynamics ((a) time series, (b) phase diagram) of the interior equilibrium point (2.6488, 5.3163) for the set of parameter values {r = 1.5, k = 6.5, d = 0.6, a = 1.2, b = 0.9, c = 0.8, m = 0.03, α = 1.7, η = 0.9, A = 1.9, h = 0.14}.

threshold value τ ∗ for which stability behaviour of the interior equilibrium point E ∗ switches to oscillatory behaviour. Hence, a supercritical Hopf-bifurcation occurs around E ∗ considering τ as a bifurcation parameter.

18

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Fig. 14. Oscillation behaviour of the population around the interior equilibrium point (2.6488, 5.3163) for α = 1.2 and all other parameters are fixed as in Fig. 13.

Fig. 15. Limit cycle around the interior equilibrium point (2.6488, 5.3163) for α = 1.2.

Fig. 16. Bifurcation diagrams of (a) prey (x) and (b) predator (y) as the quality of additional food α crosses through its critical value ‘‘α = α ∗ = 1.389’’ and other parameters are fixed as in Fig. 13.

To our knowledge, there does not exist any literature to represents the consequences of constant-yield harvesting in a predator–prey system with prey refuge proportional to both the species in presence of additional food provided to the predator. This is the first attempt to investigate the impact of constant-yield harvesting on the considered predator–prey model (26). So, in Section 6, we have modified the considered system (3) incorporating ‘‘constant-yield harvesting’’ in predator, i.e., harvested biomass is independent of predator population. We have derived the feasibility conditions and stability behaviour of all the equilibrium points under two different sets of parameters. Moreover, this work presents all possible bifurcation (saddle–node and Hopf-bifurcation) scenarios of co-dimension 1 for two distinct parameter values that can arise in the dynamics of system (26). Mathematically, it has been observed that constant-yield harvesting system brings more complicated dynamics than a system which is free from harvesting. Because harvesting may have various effects on ecosystems in general.

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

19

Fig. 17. The stability behaviour of the interior equilibrium point (9.4614, 14.2681) with respect to time t for the set of parameter values {r = 4.5, k = 9.5, d = 0.4, a = 1.3, b = 0.9, c = 0.8, m = 0.07, α = 0.7, η = 0.4, A = 1.2, h = 0.14}.

Fig. 18. Saddle–node bifurcation at m = m[SN ] = 0.014678 and all other parameters are fixed as in Fig. 17.

In Section 7, numerical simulation results have been shown briefly to verify our analytical results. Also, we have compared numerically the influence of both the states, with and without additional food (see Figs. 6, 7, 8, 9). The impact of the parameter m over the growth of the population biomass for systems (3) and (1), in presence of additional food, can be noticed in Fig. 10. From here, it can be concluded that additional food has negative effect in predator biomass, i.e., additional food does not always enhance the growth of predator. From Fig. 10, it has been also observed that the growth of predator of system (3) is comparatively large than that of system (1). But one interesting phenomenon is that predator’s growth declines quickly in system (3) than in system (1). So, the study of system (3) is very meaningful. In future work, a mathematical model can be formulated by considering nonlinear harvesting efforts for both the species. Moreover, a realistic model can be proposed to explore the effects of spatial diffusion in the pattern formation through diffusion-driven instability. Also, one can include stochastic noise in the model through the environment dependent parameters. Modification Let us modify the underlined system (3) with (4) as follows:

(

dx dt dy dt

(

= rx 1 − ca

x) k

{(

1−

=

a 1−

−

my e+y

)

(

b + αηA + 1 − my e+y

)

(

x + ηA y

b + αηA + 1 −

with initial conditions: x(0) > 0, y(0) > 0.

}

my e+y

)

− dy x

xy my e+y

)

x (37)

20

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

Fig. 19. Impact of m over the growth of population for the set of parameters: {r = 1.2, k = 6.5, d = 0.1, a = 1.1, b = 0.4, c = 0.15, α = 0.6,

η = 0.2, A = 1.2} and m ∈ [1, 2].

Here, φ (y) = e+y is a nonlinear prey refuge function, where e > 0 is half saturation constant and m is the maximum value of coefficient of refuge. Moreover, all parameters are positive and defined earlier. It is evident that system (37) is well posed. Numerical simulations provide the following important observations depicted in Fig. 19. my

Acknowledgements The authors are grateful to the anonymous referees and Dr. Eugene Stanley (Editor) for their valuable comments and helpful suggestions, which have helped them to improve the presentation of this work significantly. References [1] A.A. Berryman, The orgins and evolution of predator-prey theory, Ecology 73 (5) (1992) 1530–1535. [2] E. Beretta, Y. Kuang, Convergence results in a well-known delayed predator-prey system, J. Math. Anal. Appl. 204 (3) (1996) 840–853, http://dx.doi.org/10.1006/jmaa.1996.0471. [3] Y. Kuang, Global stability of gause-type predator-prey systems, J. Math. Biol. 28 (4) (1990) 463–474, http://dx.doi.org/10.1007/BF00178329. [4] M. Andersson, S. Erlinge, Influence of predation on rodent populations, Oikos 29 (1977) 591, http://dx.doi.org/10.2307/3543597. [5] M.E. Gilpin, M.L. Rosenzweig, Enriched predator-prey systems: Theoretical stability, Science 177 (4052) (1972) 902–904, http://dx.doi.org/10. 1126/science.177.4052.902. [6] Y. Kuznetsov, S. Rinaldi, Remarks on food chain dynamics, Math. Biosci. 134 (1) (1996) 1–33, http://dx.doi.org/10.1016/0025-5564(95)00104-2. [7] M. Kot, Elements of Mathematical Ecology, Cambridge University Press, Cambridge, 2001, http://dx.doi.org/10.1017/CBO9780511608520. [8] J.D. Murray, Mathematical Biology, Springer, New York, 1993. [9] C.S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Can. 97 (S45) (1965) 5–60, http://dx.doi.org/10.4039/entm9745fv. [10] M. Baalen, V. Krivan, P.C.J. van Rijn, M.W. Sabelis, Alternative food, switching predators, and the persistence of predatorprey systems, Amer. Nat. 157 (5) (2001) 512–524. [11] R.D. Holt, Predation, apparent competition, and the structure of prey communities, Theor. Popul. Biol. 12 (2) (1977) 197–229, http://dx.doi.org/ 10.1016/0040-5809(77)90042-9. [12] J.D. Harwood, J.J. Obrycki, in: M.s. hoddle (Ed.), The Role of Alternative Prey in Sustaining Predator Population, Vol. 2, 2005, pp. 453–462. [13] 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 (10) (2011) 2249–2276. [14] S. Mondal, A. Maiti, G.P. Sananta, Effects of fear and additional food in a delayed predator-prey model, Biophys. Rev. Lett. 13 (2018) 157–177. [15] M.P. Hassell, The Dynamics of Arthropod Predator-Prey Systems, Princeton University Press, Princeton, 1978. [16] M.P. Hassell, R.M. May, Stability in insect host-parasite models, J. Anim. Ecol. 42 (3) (1973) 693–726. [17] M. Smith, Models in Ecology, Cambridge University Press, Cambridge, 1974. [18] J.B. Collings, Bifurcation and stability analysis of a temperature-dependent mite predator-prey interaction model incorporating a prey refuge, Bull. Math. Biol. 57 (1) (1995) 63–76, http://dx.doi.org/10.1016/0092-8240(94)00024-7. [19] T.K. Kar, Stability analysis of a prey–predator model incorporating a prey refuge, Commun. Nonlinear Sci. Numer. Simul. 10 (6) (2005) 681–691, http://dx.doi.org/10.1016/j.cnsns.2003.08.006. [20] V. Krivan, Effects of optimal antipredator behavior of prey on predator–prey dynamics: The role of refuges, Theor. Popul. Biol. 53 (2) (1998) 131–142, http://dx.doi.org/10.1006/tpbi.1998.1351. [21] M.E. Hochberg, R.D. Holt, Refuge evolution and the population dynamics of coupled host-parasitoid associations, Evol. Ecol. 9 (6) (1995) 633–661, http://dx.doi.org/10.1007/BF01237660. [22] Y. Huang, F. Chen, L. Zhong, Stability analysis of a prey–predator model with holling type iii response function incorporating a prey refuge, Appl. Math. Comput. 182 (1) (2006) 672–683, http://dx.doi.org/10.1016/j.amc.2006.04.030. [23] J. Ghosh, B. Sahoo, S. Poria, Prey-predator dynamics with prey refuge providing additional food to predator, Chaos Solitons Fractals 96 (2017) 110–119, http://dx.doi.org/10.1016/j.chaos.2017.01.010. [24] S. Saha, A. Maiti, G.P. Samanta, A michaelis–menten predator–prey model with strong allee effect and disease in prey incorporating prey refuge, Int. J. Bifurcation Chaos 28 (6) (2018) 1850073, http://dx.doi.org/10.1142/S0218127418500736. [25] S. Saha, G.P. Samanta, Analysis of a predator-prey model with herd behavior and disease in prey incorporating prey refuge, Int. J. Biomath. 12 (1) (2019) 1950007. [26] E. González-Olivares, B. González-Yaez, R. Becerra-Klix, R. Ramos-Jiliberto, Multiple stable states in a model based on predator-induced defenses, Ecol. Complex. 32 (2017) 111–120, http://dx.doi.org/10.1016/j.ecocom.2017.10.004. [27] M. Haque, S. Sarwardi, Dynamics of a harvested prey–predator model with prey refuge dependent on both species, Int. J. Bifurcation Chaos 28 (12) (2018) 1830040.

S. Mondal and G.P. Samanta / Physica A 534 (2019) 122301

21

[28] A. Das, G. Samanta, Stochastic prey–predator model with additional food for predator, Physica A 512 (2018) 121–141, http://dx.doi.org/10.1016/ j.physa.2018.08.138. [29] A. Das, G.P. Samanta, Modeling the fear effect on a stochastic prey–predator system with additional food for the predator, J. Phys. A 51 (46) (2018) 465601, http://dx.doi.org/10.1088/1751-8121/aae4c6. [30] S.L. Ross, Differential Equations, John Wiley & Sons, U. K., 1984. [31] J.K. Hale, Theory of Functional Differential Equations, Springer-Verlag, New York, 1977. [32] V. Volterra, Variazionie fluttuazioni del numbero dindividui in specie animali conviventi, Mem. Acad. Lincei. 2 (1926) 31–113. [33] V. Volterra, Lecons sur la Thorie Mathematique de la lutte pour la vie, Paris, 1931. [34] Y. Kuang, Delay Differential Equation with Application in Population Dynamics, Academic Press, New York, 1993. [35] S. Nakaoka, Y. Saito, Y. Takeuchi, Stability, delay, and chaotic behavior in a lotka-volterra predator-prey system, Math. Biosci. Eng. 3 (2006) 173–187. [36] D. Xiao, S. Ruan, Multiple bifurcations in a delayed predator–prey system with nonmonotonic functional response, J. Differential Equations 176 (2) (2001) 494–510, http://dx.doi.org/10.1006/jdeq.2000.3982. [37] S. Ruan, J. Wei, On the zeros of transcendental functions with applications to stability of delay differential equations with two delays, 10 (2003) 863–874. [38] F. Brauer, A.C. Soudack, Stability regions and transition phenomena for harvested predator-prey systems, J. Math. Biol. 7 (1979) 319–337. [39] F. Brauer, A.C. Soudack, Stability regions in predator-prey systems with constant-rate prey harvesting, J. Math. Biol. 8 (1979) 55–71. [40] F. Brauer, C. Castilo-Chavez, Mathematical Models in Population Biology and Epidemiology, Springer-Verlag, New York, 2012. [41] M. Sen, P. Srinivasu, M. Banerjee, Global dynamics of an additional food provided predator–prey system with constant harvest in predators, Appl. Math. Comput. 250 (2015) 193–211, http://dx.doi.org/10.1016/j.amc.2014.10.085. [42] J. Huang, Y. Gong, S. Ruan, Bifurcation analysis in a predator-prey model with constant-yield predator harvesting, Discrete Contin. Dyn. Syst. Ser. B 18 (2013) 2101–2121, http://dx.doi.org/10.3934/dcdsb.2013.18.2101. [43] 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. [44] L. Perko, Differential Equations and Dynamical Systems, Springer-Verlag, New York, 2001.