Nonlinear Dyn https://doi.org/10.1007/s11071-018-4079-3

ORIGINAL PAPER

Dynamical behaviour of a delayed three species predator–prey model with cooperation among the prey species Soumen Kundu · Sarit Maitra

Received: 14 May 2017 / Accepted: 14 January 2018 © Springer Science+Business Media B.V., part of Springer Nature 2018

Abstract In this paper we have discussed about the dynamics of three species (two preys and one predator) delayed predator–prey model with cooperation among the preys against predation. We accept that the rate of change of density of population relies on the growth, death and in addition intra-species competition for the predators. The growth rate for preys is thought to be logistic. Delays are taken just in the growth components for each of the species. With this model we have demonstrated that the system has permanence. Taking the delays as the bifurcation parameter, the stability of the interior equilibrium point has been discussed analytically and numerically. Critical value of the delay is obtained where the Hopf-bifurcation happens. In presence of delay by constructing a Lyapunov function local asymptotic stability of the positive equilibrium point is discussed. Keywords Predator–prey logistic model · Delay · Mutualism · Permanence · Bifurcation

1 Introduction The investigations about the dynamical relationship between predators and their preys are amongst the most S. Kundu (B) · S. Maitra Department of Mathematics, National Institute of Technology, Durgapur, Durgapur 713209, India e-mail: [email protected]

essential issues in population ecology. To portray better about the growth of population in a constrained condition, Verhulst [34] presented a special sort of growth equation, named as logistic equation. Logistic equation is very useful for one animal group of population as well as more than one associating populations. Without predator the growth rate of prey population must fulfil logistic growth equation. Numerous researchers have done a lot of work on this model by distinguishing the parameters for various species [18–20,29,31]. The established logistic model is given by  x(t) , x(t) ˙ = r x(t) 1 − K

(1)

where r is the intrinsic growth rate and “K ” is the carrying capacity for x. At some points in down to earth applications like population growth, physiology of breeding, organic insusceptible reaction and so on, the growth rate of the species does not respond immediately, i.e. we need to confront a type of time delay. The presence of delay in a system greatly affects its stability. In some cases for small value of delay the system is stable and for large value of delay the stable system turns into an unstable one [2]. Hutchinson [17], first introduced the delay in a logistic differential equation. He proposed the model as:   x(t − τ ) , x(t) ˙ = r x(t) 1 − K

(2)

123

S. Kundu, S. Maitra

where τ is the time delay and other parameters have the standard significance as that of (1). Lotka and Volterra [1] gave a mathematical model for two interacting populations, where growth rate of one animal type increases in presence of other’s abatement. This model is known as Lotka–Volterra predator–prey model: dx = x(t)(a − by(t)), dt dy = y(t)(cx(t) − d), (3) dt where a, c are the growth rate of x and y, respectively, b is the rate of predation to decrease the prey growth rate, d is the common demise rate of predator species. As per the model (3), prey growth is unbounded in absence of predator, which is an unlikely supposition. Consequently, presenting carrying capacity (K ) to the growth rate of prey equation the model (3) was adjusted into the logistic Lotka–Volterra model as: x(t) dx = r x(t)(1 − ) − bx(t)y(t), dt K dy = y(t)(cx(t) − d). (4) dt There are a lot of papers [6,8,10,11,21,32,35] that have detailed discussion of at least two species for delayed or non-delayed predator–prey model. The examination on predator–prey model has relevance to subjects other than ecology, such as, game theory [16]. Other than the dynamical investigation about the stability, periodic oscillation, bifurcation, chaos, limit cycle and so forth there are numerous other dynamics to look into, for example, permanence. There are a considerable measure of articles that have the discussion about permanence in various cases like: discrete n-species food chain model with delays [7], delayed ratio-dependent predator–prey model with Holling type functional response [11], global attraction of Lotka– Volterra competitive model [5], Kolmogorov system for delayed discrete non-autonomous system [33] and extinction for non-autonomous Lotka–Volterra system [36]. Also, there are many articles [7,9,14,15,22– 24,26,30] that have the investigation about the permanence and the bifurcation for various predator–prey models. Castro et al. [4] have considered a predator– prey model with general logistic growth and exponential fading memory. First, they have demonstrated that the system has a stable solution at interior positive equilibrium. Additionally they have demonstrated that the system has a Hopf-bifurcation. Greenhalgh et al. [13]

123

have considered an eco-epidemiological predator–prey model. They assumed that the preys are influenced by diseases and the predator–prey functional responses for both the preys are taken as ratio-dependent functional responses. The boundedness for their system and stability has been checked at each of the equilibrium points. By taking search rate among suspected and infected preys they have determined the conditions for bifurcation by analytically and numerically both. Here, we take a system that have three interacting populations, two preys and one predator. The preys are picked such that there is no competition among them at any condition, while mutualism can be found among them to spare each other against predation. However in mutualism, it is seen that two species help each other to protect each other’s existence. Mutualism has awesome effects on the dynamics of numerous ecological models [12]. The mutualism between plant and animals: plant producing food for animals and animals dispersing the seeds has been reported for Zoochory plants [3]. The mutualism for growth and survivals of the anemone fish and anemones is discussed in [25]. Goh [12] discussed general nonlinear models of mutualism and prescribed simple tests for the global stability or stability in a finite region. It was pointed out that in comparison with competitions or predator–prey models, mutualistic models are more amendable to mathematical analysis. For a two species mutualism model of Lotka–Volterra type it was shown that local stability implies global stability. The system can be assumed as: dx = xg1 (x, K 1 ) − a1 x z + a2 x yz, dt dy = yg2 (y, K 2 ) − byz + a2 x yz, dt dz = −δz − c1 z 2 + c2 x z + c3 yz, dt

(5)

where K 1 , K 2 are the carrying capacities, g1 (x, K 1 ) and g2 (y, K 2 ) are the growth functions with regard to preys x and y. The coefficients a1 , b are the rate of predation, and a2 is the rate of cooperation for the preys x and y, individually, δ is the intrinsic mortality rate, and c1 is the rate of intra-species competition within the predator species z. The rate of transformation because of predation from preys to predators is c2 x z and c3 yz. Here, we assume that the rate of cooperation is not as much as the rate of predation, i.e. x < ab2 and y < aa21 . Presently we consider the delay logistic growth functions for the preys x

Dynamical behaviour of a delayed three species predator–prey model

and y, i.e. we pick g1 (x, K 1 ) = l1 (1 − y(t−τ2 ) K2 )

x(t−τ1 ) K1 )

and

then above system (5) g2 (y, K 2 ) = l2 (1 − becomes as:   x(t − τ1 ) dx − a1 x z + a2 x yz, = l1 x 1 − dt K1   dy y(t − τ2 ) − byz + a2 x yz, = l2 y 1 − dt K2 dz = −δz − c1 z 2 + c2 x(t − τ3 )z + c3 y(t − τ3 )z, dt (6)

Theorem 2.1 Let M1 , M2 , M3 , m 1 , m 2 , m 3 are positive constants and independent of the initial solution of system (6), if the following conditions hold

subjected to the initial condition x(θ ) = φ1 (θ ) > 0, θ ∈ [−τ1 , 0), φ1 (0) > 0, y(θ ) = φ2 (θ ) > 0, θ ∈ [−τ2 , 0), φ2 (0) > 0, z(θ ) = φ3 (θ ) > 0, θ ∈ [−τ3 , 0), φ3 (0) > 0,

Lemma 2.1 (see [2,35]) If for t ≥ 0 and x(0) ≥ 0 we have x˙ ≥ x(c −dx) where c > 0, d > 0 then c lim inf x(t) ≥ t→∞ d and if for t ≥ 0 and x(0) ≥ 0 we have x˙ ≤ x(c −dx) where c > 0, d > 0 then c lim sup x(t) ≤ . t→∞ d

(7)

where all the parameters are assumed to be positive and l1 , l2 are the intrinsic growth rates for the preys x, y, respectively. With this model we want to see the effect of delay on permanence and stability of the system in presence of cooperation. In this paper, first, in presence of delay it is shown that the system (6) has permanence (Sect. 2). In Sect. 3, the condition for stability and presence of bifurcation has been examined by taking delays as the bifurcation parameter. In Sect. 4, local stability of the system (6) has been discussed by forming a Lyapunov function. Some numerical results have been incorporated for our expository results. The impact of each parameter is additionally indicated numerically in both absence and presence of delay.

2 Permanence

m 1 ≤ lim inf x(t) ≤ lim sup x(t) ≤ M1 , t→∞

t→∞

m 2 ≤ lim inf y(t) ≤ lim sup y(t) ≤ M2 , t→∞

t→∞

m 3 ≤ lim inf z(t) ≤ lim sup z(t) ≤ M3 , t→∞

t→∞

with positive initial condition, i.e. x(0) > 0, y(0) > 0, z(0) > 0, then we say that the system (6) has permanence. Proof With the positive initial condition (x(0), y(0), z(0)), it is easy to see that the solution (x(t), y(t), z(t)) of the system (6) is positive. From the first equation of system (6) we can write dx ≤ l1 x. dt

Integrating both sides of (8) within t − τ1 to t we get  t x(t) ≤ l1 dt = l1 τ1 ln x(t − τ1 ) t−τ1 ⇒ x(t − τ1 ) ≥ x(t)e−l1 τ1 .

One can check that the system (6) has a positive solution with the positive initial condition given in (7). Now we show that the system (6) has permanence. Before proving, we first state the definition of permanence and a lemma which will be used to prove our statement given below in Theorem 2.1. Definition 2.1 (see [2,35]) A system is said to has permanence if for positive constants m and M we have m ≤ lim inf x(t) ≤ lim sup x(t) ≤ M

(9)

From the first equation of (6) and using (9) we get dx l1 −l1 τ1 ≤ x(t){l1 − e x(t)}. dt K1

(10)

Following Lemma 2.1 we can say lim sup x(t) ≤ K 1 el1 τ1 = M1 ,

t→∞

t→∞

(8)

(11)

i.e. for  > 0, ∃ a T1 > 0 s.t. x(t) ≤ M1 + , ∀t > T1 Similarly following the same computation as done for first equation of (6) we get from the second equation of (6) we get

t→∞

for all positive solutions x(t) of the system.

lim sup y(t) ≤ K 2 el2 τ2 = M2 .

t→∞

(12)

123

S. Kundu, S. Maitra

Therefore, y(t) ≤ M2 + , ∀t > T2 and from the third equation of (6) we have lim sup z(t) ≤ M3 ,

(13)

t→∞

lim inf z(t) ≥ m 3 ,

(19)

t→∞

where

1 c2 m 11 + c3 m 12 − δ , c1  1 1 K 1 {l1 − a1 (M3 + )} exp m1 = l1    l1 (M1 + ) − a1 (M3 + ) τ3 , × l1 − K1

where M3 =

Finally from the third equation of (6), we can write

−δ

+ c2 (M11

+ ) + c3 (M21 c1

+ )

m3 =

,

M11 = K 1 el1 τ3 and

M21 = K 2 el2 τ3 thus,

z(t) ≤ M3 + , ∀t > T3 .

and

Again from the first equation of (6) we can write

m 12

  dx l1 ≥ x l1 − (M1 + ) − a1 (M3 + ) , dt K1

(14)

integrating both sides of (14) within [t − τ1 , t] we get   l1 (M1 + ) x(t − τ1 ) ≥ x(t) exp − l1 − K1   (15) − a1 (M3 + ) τ1 . ∴ From first equation of (6) and using (15) we can write     dx l1 l1 ≥ x(t) l1 − exp − l1 − (M1 + ) dt K1 K1     − a1 (M3 + ) τ1 x(t) − a1 (M3 + ) .

 1 K 2 {l2 − b(M3 + )} exp = l2    l2 (M2 + ) − b(M3 + ) τ3 . × l2 − K2

Hence from (11) & (17), (12) & (18), (13) & (19) we have m 1 ≤ lim inf x(t) ≤ lim sup x(t) ≤ M1 , t→∞

t→∞

m 2 ≤ lim inf y(t) ≤ lim sup y(t) ≤ M2 , t→∞

t→∞

m 3 ≤ lim inf z(t) ≤ lim sup z(t) ≤ M3 , t→∞

t→∞

with positive initial condition, i.e. x(0) > 0, y(0) > 0, z(0) > 0. Thus, we conclude that the system (6) has permanence. 

(16) 3 Stability and bifurcation analysis

Using Lemma 2.1 we can write lim inf x(t) ≥ m 1 ,

t→∞

where m1 =

(17)

 1 K 1 {l1 − a1 (M3 + )} exp l1    l1 (M1 + ) − a1 (M3 + ) τ1 . × l1 − K1

Similarly from second equation of (6) we get lim inf y(t) ≥ m 2 ,

t→∞

where m2 =

(18)

 1 K 2 {l2 − b(M3 + )} exp l2    l2 (M2 + ) − b(M3 + ) τ2 . × l2 − K2

123

The subject of bifurcation is a vast area to do research. In a system of differential equations Hopf-bifurcation happens when the complex conjugate set of eigenvalues of a linearised system become purely imaginary at a fixed point. So with these conditions, a system that contains at least two equations may have Hopfbifurcation. Hopf-bifurcation occurs at a point where a system changes state from stable to unstable, i.e. it is a local bifurcation in which a fixed point of a dynamical system loses stability. In [27, p. 264], it is given that “the Hopf-bifurcation mainly happens when a periodic solution or limit cycle encompassing an equilibrium point emerges or leaves as parameters differs. At the point when a stable limit cycle surrounds an unstable equilibrium point, the bifurcation is called supercritical Hopf-bifurcation. If the limit cycle is unstable and

Dynamical behaviour of a delayed three species predator–prey model

surrounds a stable equilibrium point, at that point the bifurcation is called sub-critical Hopf-bifurcation.” In our previous paper [21], taking τ1 = 0, τ2 = 0 and τ3 = τ the stability analysis for the model (6) has been derived by constructing a suitable Lyapunov functional. Later, by taking τ as bifurcation parameter, numerically we have obtained a critical value of τ where the Hopfbifurcation exists. Also, the effect of the cooperation on the stability of the system had been shown. In this section by taking all the delays τ1 , τ2 , τ3 as bifurcation parameters the necessary conditions for the Hopf-bifurcation are derived. Let the system (6) has a positive interior equilibrium E(x ∗ , y ∗ , z ∗ ). Letting v1 (t) = x − x ∗ , v2 (t) = y − y ∗ , v3 (t) = z − z ∗ and substituting into Eq. (6) we get the linearised form as:

We can say that E(x ∗ , y ∗ , z ∗ ) is locally asymptotically stable (Fig. 1) if

dv1 = − p11 e−λτ1 + p12 v2 + p13 v3 , dt dv2 = p21 v1 − p22 e−λτ2 v2 + p23 v3 , dt dv3 = p31 e−λτ3 + p32 e−λτ3 − p33 v3 , dt

(λ3 + m 11 λ2 + m 12 λ + m 13 )

l1 x ∗ ∗ ∗ K 1 , p12 = a2 x z , p13 = ∗ a1 x ∗ , p21 = a2 y ∗ z ∗ , p22 = l2Ky2 , p23 = by ∗ , p31 = c2 z ∗ , p32 = c3 z ∗ , p33 = c1 z ∗ .

where p11 =

p33 + p22 + p11 > 0, p11 p22 p33 − p12 p21 p33 − p12 p23 p31 − p13 p21 p32 − p11 p23 p32 − p13 p22 p31 > 0 and ( p11 + p22 + p33 )(− p12 p21 + p11 p33 + p22 p33 − p23 p32 − p13 p31 + p11 p22 ) − ( p11 p22 p33 − p12 p21 p33 − p12 p23 p31 − p13 p21 p32 − p11 p23 p32 − p13 p22 p31 ) > 0, i.e. the roots of (22) have negative real part. Case 2 when τ1 = 0, τ2 = 0, τ3 = 0, (21) becomes +e−λτ1 (n 11 λ2 + n 12 λ + n 13 ) = 0,

(20) a2 x ∗ y ∗ − a2 x ∗ y ∗ −

The characteristic equation for the system (20) as given below

where m 11 = p22 + p33 , m 12 = p22 p33 − p12 p21 − p23 p32 − p13 p31 , m 13 = − p12 p21 p33 − p12 p23 p31 − p13 p21 p32 − p13 p22 p31 , n 11 = p11 , n 12 = p11 p33 + p11 p22 , n 13 = p11 p22 p33 − p11 p23 p32 . Let λ = iω(> 0) then from (23) separating real and imaginary part we get m 13 − m 11 ω2 = −n 12 ω sin(ωτ1 ) + (n 13 − n 11 ω2 ) cos(ωτ1 ), m 12 − ω3 = (n 13 − n 11 ω2 ) sin(ωτ1 )

(λ3 + p33 λ2 − λp12 p21 − p12 p21 p33 )

+ n 12 ω cos(ωτ1 ),

+ e−λτ1 (λp11 p33 + p11 λ2 ) + e−λτ2 ( p22 p33 λ + p22 λ2 ) +e

−λτ3

ω6 + (m 211 − 2m 12 + n 211 )ω4

− p13 p21 p32 − p13 p31 λ) +e

( p11 p22 p33 + p11 p22 λ)

+e

−λ(τ1 +τ3 )

(− p11 p23 p32 )

+ e−λ(τ2 +τ3 ) (− p13 p22 p31 ) = 0.

+ (m 212 − 2m 11 m 13 + 2n 11 n 13 − n 212 )ω2 + (m 213 − n 213 ) = 0. (21)

Now we discuss about the bifurcation conditions. Case 1 when all the delays are assumed to be zero, i.e. τ1 = 0, τ2 = 0, τ3 = 0, (21) becomes λ3 + ( p33 + p22 + p11 )λ2 + (− p12 p21 + p11 p33 + p22 p33 − p23 p32 − p13 p31 + p11 p22 )λ + ( p11 p22 p33 − p12 p21 p33 − p12 p23 p31 − p13 p21 p32 − p11 p23 p32 − p13 p22 p31 ) = 0.

(24)

which gives

(−λp23 p32 − p12 p23 p31

−λ(τ1 +τ2 )

(23)

(22)

(25)

Now Eq. (25) has positive root if m 211 − 2m 12 + n 211 > 0 and m 213 − n 213 < 0. Eliminating sin(ωτ1 ) from (24) and substituting ω = ω0 , where ω0 is a positive root of (25), we get  (m 11 ω02 − m 13 )(n 11 ω02 − n 13 ) 1 τ1n = arccos ω0 n 212 ω02 + (n 11 ω02 − n 13 )2

n 12 ω0 (ω03 − m 12 ω0 ) − 2 2 n 12 ω0 + (n 11 ω02 − n 13 )2 + 2nπ } , n = 0, 1, 2, 3, . . .

(26)

123

S. Kundu, S. Maitra

Using Butler’s Lemma, E(x ∗ , y ∗ , z ∗ ) will be stable for τ1 < τ10 (see Fig. 2), where τ10 is obtained from (26) by taking n = 0. Now for τ1 > τ10 , E(x ∗ , y ∗ , z ∗ ) will be unstable (see Fig. 5) and for τ1 = τ10 , it gives a periodic solution (see Figs. 3, 4 ). Now we discuss about the existence of Hopf-bifurcation with respect to the bifurcation parameter τ1 , i.e. if we dλ −1 ) > 0 for ω = ω0 , τ1 = τ1n , then can show Rel( dτ 1 we can say that Hopf-bifurcation occurs for τ1 . Now differentiating (23) with respect to τ1 and then substituting λ = iω0 and simplifying we get  dλ −1 dτ1   1 2ω06 +(m 211 − 2m 12 −n 211 )ω04 +(n 213 −m 213 )

dλ −1 ) > 0 as m 221 − 2m 22 − n 221 > 0 Therefore, Rel( dτ 2 and n 223 − m 223 > 0. Thus, Hopf-bifurcation occurs for τ2 .

Case 4 when τ1 = 0, τ2 = 0, τ3 = 0 following the same computation as done in Case 2 we get τ3n =



Rel =

w02

n 211 ω04 +(n 212 −2n 11 n 13 )ω02 +n 213

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

(27) dλ −1 Therefore, Rel( dτ ) 1 2 2 0 and n 13 − m 13 > 0. for τ1 .

> 0 as − 2m 12 − > Thus Hopf-bifurcation occurs m 211

n 211

Case 3 when τ1 = 0, τ2 = 0, τ3 = 0 following the same computation as done in Case 2 we get  (m 21 ω02 − m 23 )(n 21 ω02 − n 23 ) 1 arccos τ2n = ω0 n 222 ω02 + (n 21 ω02 − n 23 )2

n 22 ω0 (ω03 − m 22 ω0 ) − 2 2 n 22 ω0 + (n 21 ω02 − n 23 )2 + 2nπ } , n = 0, 1, 2, 3, . . .

(28)

where m 21 = p11 + p33 , m 22 = − p12 p21 + p11 p33 − p23 p32 − p13 p31 , m 23 = − p12 p21 p33 − p12 p23 p31 − p13 p21 p32 − p11 p23 p32 , n 21 = p22 , n 22 = p22 p33 + p11 p22 , n 23 = p11 p22 p33 − p13 p22 p31 . Thus, E(x ∗ , y ∗ , z ∗ ) will be stable for τ2 < τ20 (see Fig. 6), where τ20 is obtained from (28) by taking n = 0. Now for τ2 > τ20 , E(x ∗ , y ∗ , z ∗ ) will be unstable (see Fig. 9) and for τ2 = τ20 it gives a periodic solution (see Figs. 7, 8 ). dλ −1 ) > 0 for ω = Also it can be shown that Rel( dτ 2 ω0 , τ2 = τ2n , i.e. Hopf-bifurcation occurs for τ2 , where  dλ −1 dτ2   1 2ω06 +(m 221 −2m 22 −n 221 )ω04 +(n 223 −m 223 )

 Rel =

w02

n 221 ω04 +(n 222 −2n 21 n 23 )ω02 +n 223

.

(29)

123

1 arccos ω0   (m 31 ω02 −m 33 )n 32 +n 31 ω02 (ω02 −m 32 ) +2nπ , × n 232 +n 31 ω02

(30)

where m 31 = p11 + p22 + p33 , m 32 = − p12 p21 + p11 p33 + p22 p33 + p11 p22 , m 33 = p12 p21 p33 + p11 p22 p33 , n 31 = − p23 p32 − p13 p31 , n 32 = − p12 p23 p31 − p13 p21 p32 − p11 p23 p32 − p13 p22 p31 . Thus E(x ∗ , y ∗ , z ∗ ) will be stable for τ3 < τ30 (see Fig. 10), where τ30 is obtained from (30) by taking n = 0. Now for τ3 > τ30 , E(x ∗ , y ∗ , z ∗ ) will be unstable (see Fig. 13) and for τ3 = τ30 it gives a periodic solution (see Figs. 11, 12). dλ −1 ) > 0 for ω = Also it can be shown that Rel( dτ 3 ω0 , τ3 = τ3n , i.e. Hopf-bifurcation occurs for τ3 , where  dλ −1 dτ3   2ω02 +(m 231 −2m 32 )ω04 +(n 232 −m 233 ) 1 = 2 arccos . w0 n 232 +n 231 ω02 

Rel

(31) dλ −1 Therefore, Rel( dτ ) > 0 as m 231 − 2m 32 > 0 and 3 m 233 − n 232 < 0. Thus, Hopf-bifurcation occurs for τ3 .

Case 5 τ1 = 0, τ2 = 0, τ3 = 0 Since all the delays are assumed as bifurcation parameters and in this case they are all nonzero, it is difficult to find the critical value at which Hopf-bifurcation occurs. Therefore, by taking any one of them (delays) as a bifurcation parameter (when others are taken in their stable interval) we will derive the necessary condition(s) for Hopf-bifurcation. First, we take τ1 as bifurcation parameter and τ2 , τ3 are in stable interval. Letting iω as a root of (21) we get

2 2 ω4 + 2d2 + d12 − d41 ω6 + 2d1 ω5 + d21 + (2d3 + 2d21 d22 + 2d1 d2 − 2d41 d42 )ω3

Dynamical behaviour of a delayed three species predator–prey model

2 + 2d21 d23 + 2d1 d3 + d22 + d22 − 2d41 d43

2 2 ω2 − d42 − d31 + (2d22 d23 + 2d2 d3 − 2d42 d43 − 2d31 d32 )ω

2 2 2 + d32 + d23 = 0, − d32 − d43

(32)

where

Now by taking τ1 , τ3 and τ1 , τ2 , respectively, in their stable interval and following the similar calculation as done to get τ1 = τ1∗ , one can easily get the critical values for τ2 = τ2∗ and τ3 = τ3∗ , respectively, where the Hopf-bifurcations occur.

4 Lyapunov stability analysis

d1 = − p22 sin ωτ2 , Here, we have studied the stability of (6) by using a suitable Lyapunov functional as done in [21]. The system (6) has positive equilibrium point E 6 = (x ∗ , y ∗ , z ∗ ). Introducing the new set of variables v1 (t) = x(t) − x ∗ , v2 (t) = y(t) − y ∗ and v3 (t) = z(t) − z ∗ in Eq. (6), we get the following linearised system:

d2 = − p22 p33 cos ωτ2 + p23 p32 cos ωτ3 + p13 p31 cos ωτ3 , d3 = p12 p23 p31 sin ωτ3 + p13 p21 p32 sin ωτ3 − p13 p22 p31 sin ω(τ2 + τ3 ), d21 = p33 + p22 cos ωτ2 ,

l1 v1 + a2 z ∗ v2 + a2 y ∗ v3 , K1 l2 (B(t)) = −bv3 − v2 + a2 z ∗ v1 + a2 x ∗ v3 , K2 (C(t)) = −c1 v3 + c2 v1 + c3 v2 ,

d22 = − p22 p33 sin ωτ2

(A(t)) = −a1 v3 −

+ p23 p32 sin ωτ3 + p13 p31 sin ωτ3 , d23 = p12 p23 p31 cos ωτ3 + p13 p21 p32 cos ωτ3 + p13 p22 p31 cos ω(τ2 + τ3 ), d31 = p11 p33 − p11 p22 cos ωτ2 ,

where

 t v1 l1 − v1 (s)ds, x∗ K 1 t−τ1  t l2 v2 v2 (s)ds B(t) = ∗ − y K 2 t−τ2

d32 = p11 p23 p32 sin ωτ3 − p11 p22 p33 sin ωτ2 ,

A(t) =

d41 = p11 , d42 = p11 p22 sin ωτ2 , d43 = p11 p23 p32 cos ωτ3 − p11 p22 p33 cos ωτ2 .

and Applying Descartes’ rule of signs we can say that (32) has at least one positive root (say ω) if 2 2 2 < d32 + d43 d32 + d23

(33)

holds. Following the earlier steps as done in Case 2 one can easily calculate the critical value for τ1 at which Hopf-bifurcation occurs. Thus, for finite number of positive ω we get finite number of τ1i , i = 1, 2 . . .. Let τ1∗ = min{τ1i , i = 1, 2 . . .}, then   d(Reλ) > 0. dτ1 λ=iω,τ1 =τ ∗

(36)

(34)

v3 C(t) = ∗ + c2 z



t t−τ3

 v1 (s)ds + c3

t

v2 (s)ds.

t−τ3

Now following the steps as in [21,28], we shall check the stability of the system by assuming a suitable Lyapunov function w(v)(t) as follows: w(v)(t) = k1 w1 (v)(t) + k2 w2 (v)(t) + k3 w3 (v)(t) + k4 w4 (v)(t) + k5 w5 (v)(t) + k6 w6 (v)(t), (37) where k1 , k2 , k3 , k4 , k5 , k6 are given in Appendix and w1 (v)(t)

  l1 l1 ∗ ∗ a1 − a2 y − a2 z + = A (t) + K1 K1  t  t v12 (l)dlds, t−τ1 s   l2 l2 b − a2 z ∗ − a2 x ∗ + w2 (v)(t) = B 2 (t) + K2 K2  t  t v22 (l)dlds, 2

(35)

1

Thus, we see that Hopf-bifurcation occurs for τ1 = τ1∗ , when τ2 , τ3 are in stable interval, i.e. τ2 ∈ [0, τ2∗ ] and τ3 ∈ [0, τ3∗ ].

t−τ2

s

123

S. Kundu, S. Maitra

w3 (v)(t) = C 2 (t)

  + (c2 + c3 − c1 ) c2  + c3

t−τ3



t

t−τ3

t

t s

 v22 (l)dlds ,

 s

t

v12 (l)dlds

w4 (v)(t)

 1 l1 (c1 − c2 − c3 ) = A(t)C(t) + 2 K1   ∗ ∗ + a2 c2 z − a1 c2 + a2 c2 y   t  t 1 2 a2 c3 z ∗ v1 (l)dlds + × 2 t−τ1 s  l1 c3 ∗ + a2 c3 y − a1 c3 − K1  t  t v22 (l)dlds, × t−τ3

Fig. 1 Time versus population have been plotted by taking τ1 = 0, τ2 = 0, τ3 = 0 for different values of a2 , which represents a stable solution

s

w5 (v)(t)

l1  b − a2 x ∗ 2K 1  t  t l2 ∗ + − a2 z v12 (l)dlds K2 t−τ1 s   l2 l1 ∗ ∗ a1 − a2 y − a2 z + + 2K 2 K1  t  t v22 (l)dlds, ×

= A(t)B(t) +

t−τ2

s

w6 (v)(t) = B(t)C(t)   l2 c2 1 a2 c2 x ∗ − bc2 − + a2 c2 z ∗ + 2 K2  t  t v12 (l)dlds × t−τ3

s

1 + a2 c3 z ∗ + a2 c3 x ∗ − bc3 2  t  t l2 + (c1 −c2 −c3 ) v22 (l)dlds, K2 t−τ2 s as all the parameters are assumed positive so, k1 , k2 , k3 , k4 , k5 , k6 > 0 and w(v)(t) > 0. Taking the derivative of (37), we get d w(v)(t) ≤ Λ1 v12 + Λ2 v22 + Λ3 v32 , (38) dt where expression for Λ1 , Λ2 , Λ3 are given in Appendix.

123

Fig. 2 Stable solution for the population when τ1 is the bifurcation parameter. τ1 < τ10 (τ2 = 0, τ3 = 0), other parameters are same as that of Fig. 1

Theorem 4.1 If the value of the delays τ1 , τ2 , τ3 satisfies the conditions Λ1 < 0, Λ2 < 0, Λ3 < 0 then the interior equilibrium point E 6 of (6) is locally asymptotically stable. Following the steps as done in [21,28], one can easily prove Theorem 4.1.

5 Numerical analysis In the previous section, the conditions for stability and Hopf-bifurcation of the system (6) have been derived analytically. Now we do some numerical computations of system (6) to understand our results obtained in Sect. 3 by choosing suitable values of the parameters. For different values of delays we have obtained different scenarios with E(x ∗ , y ∗ , z ∗ ) as interior equilibrium point. The values of the parameters are taken as: l1 =

Dynamical behaviour of a delayed three species predator–prey model

Fig. 3 Denote a periodic solution around E(x ∗ , y ∗ , z ∗ ) at τ1 = τ10 (τ2 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

Fig. 4 Stable limit cycle around E(x ∗ , y ∗ , z ∗ ) at τ1 = τ10 (τ2 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

Fig. 5 Unstable solution when τ1 > τ10 (τ2 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

2.0, l2 = 2.0, a1 = 1.3, a2 = 1.2, K 1 = 0.7, K 2 = 0.6, b = 1.5, δ = 0.2, c1 = 0.8, c2 = 1.0, c3 = 1.4. In Fig. 1, a stable solution for the species (6) has been plotted by taking all the delays equal to zero. For increasing value of cooperation (a2 ) we can see that the density of the population also increases. Now if we consider the nonzero value of delay, corresponding results are discussed in what follows:

Fig. 6 Stable solution for the population when τ2 is the bifurcation parameter. τ2 < τ20 (τ1 = 0, τ3 = 0), other parameters are same as that of Fig. 1

Fig. 7 Periodic solution around E(x ∗ , y ∗ , z ∗ ) at τ2 = τ20 (τ1 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

Fig. 8 Stable limit cycle around E(x ∗ , y ∗ , z ∗ ) at τ2 = τ20 (τ1 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

First, by taking τ1 = 0(τ2 = 0, τ3 = 0), τ2 = 0(τ1 = 0, τ3 = 0) and τ3 = 0(τ1 = 0, τ2 = 0), respectively, and with the above set of parameters we obtained the critical value of delays τ10 = 1.32(from 26), τ20 = 1.2449 (from 28) and τ30 = 5.5061 (from 30). The interior equilibrium point E(x ∗ , y ∗ , z ∗ ) is seen to be stable for less values of the delays (Figs. 2, 6, 10). Now with the increased values of delays the oscillation of the population also increases. There-

123

S. Kundu, S. Maitra

Fig. 12 Stable limit cycle around E(x ∗ , y ∗ , z ∗ ) at τ3 = τ30 (τ1 = 0, τ2 = 0), when other parameters are same as that of Fig. 1 Fig. 9 Unstable solution when τ2 > τ20 (τ1 = 0, τ3 = 0), when other parameters are same as that of Fig. 1

Fig. 13 Unstable solution when τ3 > τ30 (τ1 = 0, τ2 = 0), when other parameters are same as that of Fig. 1 Fig. 10 Stable solution for the population when τ3 is the bifurcation parameter. τ3 < τ30 (τ1 = 0, τ2 = 0), other parameters are same as that of Fig. 1

Fig. 14 Stable solution when τ1 < 1.1, τ2 < 1.1, τ3 < 1.1, other parameters are same as that of Fig. 1 Fig. 11 Periodic solution around E(x ∗ , y ∗ , z ∗ ) at τ3 = τ30 (τ1 = 0, τ2 = 0), when other parameters are same as that of Fig. 1

fore, at the critical value of delays we get a stable periodic solution where the Hopf-bifurcation occurs (Figs. 3, 7, 11) and there we get a stable limit cycle around E(x ∗ , y ∗ , z ∗ ) (Figs. 4, 8, 12). For large values of delays (τ10 > 1.32, τ20 > 1.2449, τ30 > 5.5061) the system becomes unstable (Figs. 5, 9, 13). Now we con-

123

sider the case when all the delays are nonzero and they are equal. Here for τ1 < 1.1, τ2 < 1.1, τ3 < 1.1, satisfying the relation (33) we get a stable solution (Fig. 14), and at τ1 = 1.1, τ2 = 1.1, τ3 = 1.1, we get periodic solution (Figs. 15, 16) and for τ1 > 1.1, τ2 > 1.1, τ3 > 1.1 we get an unstable solution (Fig. 17). Now considering system (6) with τ2 , τ3 in its stable interval and τ1 as bifurcation parameter we get a stable solution for τ1 = 0.7, τ2 = 0.5, τ3 = 0.6 (Fig. 18), a periodic solu-

Dynamical behaviour of a delayed three species predator–prey model

Fig. 15 Periodic solution when τ1 = 1.1 = τ2 = τ3 , other parameters are same as that of Fig. 1

Fig. 16 Stable limit cycle around E(x ∗ , y ∗ , z ∗ ) at τ1 = 1.1 = τ2 = τ3 , other parameters are same as that of Fig. 1

Fig. 18 Stable solution for τ1 = 0.7, τ2 = 0.5, τ3 = 0.6, other parameters are same as that of Fig. 1

Fig. 19 Periodic solution around E(x ∗ , y ∗ , z ∗ ) for τ1 = 0.9, τ2 = 0.5, τ3 = 0.6, other parameters are same as that of Fig. 1

Fig. 20 Stable limit cycle for τ1 = 0.9, τ2 = 0.5, τ3 = 0.6, other parameters are same as that of Fig. 1

Fig. 17 Unstable solution for τ1 = τ2 = τ3 > 1.1, other parameters are same as that of Fig. 1

tion for τ1 = in Fig. 25.

12 sin t 10 , τ2

= 1−sin t, τ3 = cos t is shown

tion for τ1 = 0.9, τ2 = 0.5, τ3 = 0.6 (Figs. 19, 20), unstable solution for τ1 = 1.25, τ2 = 0.5, τ3 = 0.6 (Fig. 21). Keeping the value of the parameters same as that of Fig. 18 and taking delay as variable, we get stable solution (Fig. 22) for τ1 = sin2 t , τ2 = 1 − sin t, τ3 = cos t. Here τ1 is changed gradually. t Periodic solution is obtained for τ1 = 9 sin 10 , τ2 = 1 − sin t, τ3 = cos t (Figs. 23, 24), and unstable solu-

6 Discussion of our model In this section we shall discuss about the result obtained by changing the parameter values as well as the effect of predation in presence of cooperation between the prey species. At first we assume the model has no predator and the preys have equal carrying capacities (K 1 = K 2 ) and equal logistic growth rates (l1 = l2 ). In that case

123

S. Kundu, S. Maitra

Fig. 21 Unstable solution for τ1 = 1.25, τ2 = 0.5, τ3 = 0.6, other parameters are same as that of Fig. 1

Fig. 22 Solution for τ1 = sin2 t , τ2 = 1 − sin t, τ3 = cos t, other parameters are same as that of Fig. 1

t Fig. 23 Periodic solution for τ1 = 9 sin 10 , τ2 = 1 − sin t, τ3 = cos t, other parameters are same as that of Fig. 1

the two species increase simultaneously approaching the carrying capacities of the species (Fig. 26). This assumption to have exactly equal K 1 , K 2 and l1 , l2 for

123

t Fig. 24 Stable limit cycle for τ1 = 9 sin 10 , τ2 = 1 − sin t, τ3 = cos t, other parameters are same as that of Fig. 1

sin t Fig. 25 Unstable solution for τ1 = 1210 , τ2 = 1 − sin t, τ3 = cos t, other parameters are same as that of Fig. 1

Fig. 26 Time series plot with equal growth rates, equal carrying capacities between two preys and there is no predator

two species is not reasonable to happen in an ecological system. Now in the presence of predator (equal K 1 , K 2 and equal l1 , l2 for the prey species), but with no delay, we see that all the species decrease and after

Dynamical behaviour of a delayed three species predator–prey model

Fig. 27 Time series plot with equal growth rates, equal carrying capacities between two preys and rate of predation for prey species y is greater than the species x with τ1 = 0, τ2 = 0, τ3 = 0

Fig. 28 Time series plot with equal growth rates, equal carrying capacities between two preys and rate of predation for prey species y is greater than the species x with τ1 = 0, τ2 = 0, τ3 = 0

some time they reach to a stable situation (Fig. 27). But in the presence of delay, there are oscillations for all the species (Fig. 28), where it is assumed that the rate of predation for the species y (i.e.b) is greater than the rate of predation for the species x (i.e. a1 ). Here, as predation rate of the species y is greater than the species x, so number of species x is larger than the species y. Now in presence of predator, we assume that the logistic growth rate of prey species y is greater than the growth rate of prey species x (l2 > l1 ). As a result (when τi = 0, i = 1, 2, 3) it is seen in Fig. 29, that y survives better and remain larger than x. In Fig. 30, when τi =

Fig. 29 Time series plot with equal predation, equal carrying capacities between two preys and growth rate for prey species y is greater than the species x with τ1 = 0, τ2 = 0, τ3 = 0

Fig. 30 Time series plot with equal predation, equal carrying capacities between two preys and growth rate for prey species y is greater than the species x with τ1 = 0, τ2 = 0, τ3 = 0

0, i = 1, 2, 3, there are oscillations of all the species. From Fig. 31, it is seen that if τi = 0, i = 1, 2, 3 and l1 = l2 , a1 = b, but K 1 > K 2 , then x remain larger than y. Figure 32 shows that, for τi = 0, i = 1, 2, 3, all the species have oscillations and as time gets large x, y, z stabilise. The effect of the rate of cooperation (i.e. a2 ) is shown in Figs. 33 and 34. In Fig. 33 with a2 = 1.2, initially the species x is increasing and y, z are decreasing. After some time all the species come to a stable solution. In Fig. 34, a2 is greater than 1.2. Here though initially both the species y and z decrease, after some time they monotonically increase. The species x is seen to increase monotonically from the beginning. Thus, the number of species becomes unbounded.

123

S. Kundu, S. Maitra

Fig. 31 Time series plot of the system (6) with equal predation, equal growth rates but K 1 > K 2 when τ1 = 0, τ2 = 0, τ3 = 0

Fig. 34 Time series plot of the system (6) with equal predation, equal growth rates, equal carrying capacities between preys when τ1 = 0, τ2 = 0, τ3 = 0

7 Conclusion

Fig. 32 Time series plot of the system (6) with equal predation, equal growth rates but K 1 > K 2 when τ1 = 0, τ2 = 0, τ3 = 0

Fig. 33 Time series plot of the system (6) with equal predation, equal growth rates, equal carrying capacities between preys when τ1 = 0, τ2 = 0, τ3 = 0

123

In this paper we have investigated the dynamical behaviours of a three species (two preys, one predator) delayed predator–prey model. In Sect. 2, with this model we have derived sufficient condition to show the system has permanence. It is shown that the delays have a great impact on permanence, i.e. delay is an important factor to decide the permanence and stability of the system. Then in absence of delay, we have obtained the condition for local stability of the interior equilibrium point E(x ∗ , y ∗ , z ∗ ) (Sect. 3, Case 1). In presence of delays, we obtained a critical value of delays where the Hopf-bifurcation occurs. By taking τ1 = 0(τ2 = 0, τ3 = 0), τ2 = 0(τ1 = 0, τ3 = 0) and τ3 = 0(τ1 = 0, τ2 = 0), respectively (Cases 2, 3, 4), as bifurcation parameter, condition for existence of Hopfbifurcation has been obtained. By taking τ1 = 0, τ2 = 0, τ3 = 0, the local stability has been discussed analytically by taking a suitable Lyapunov function. In numerical section, taking particular values of model parameters and assuming delays as bifurcation parameters some numerical simulation have been done to understand our results better. Also, for τ1 = 0, τ2 = 0 and τ3 = 0 (Case 5) both analytical and numerical simulations are done to discuss the existence of Hopf-bifurcation with respect to the positive equilibrium E(x ∗ , y ∗ , z ∗ ). The effect of variable delay has been shown by numerical results. Finally in Sect. 6, the effect of model parameters has been shown in both presence and absence of delays.

Dynamical behaviour of a delayed three species predator–prey model Acknowledgements We thank the anonymous referee for valuable suggestions. The first author is thankful to DST, New Delhi, India, for its financial support under INSPIRE fellowship, without which this research would not have been possible.

Appendix

1 1 1 , M = ∗, N = ∗ x∗ y z 1 k1 = + (c2 L + l1 c2 τ1 2l1 K 1 L c1 c2 τ3 + 2c22 τ3 + 2c2 c3 τ3 a1 c2 τ1 + + 2 2 l2 c2 τ1 K 2 a2 L bc2 τ3 l1 c3 τ3 + + + + N 2K 1 2 2 a2 c2 τ3 a2 c2 τ3 + + 2L 2M  a2 c3 τ3  2a2 c2 τ3 + + (2l1 L (−a1 N N 2N a2 N a1 c2 τ3 a2 N + + − bN + L M 2 a1 c2 τ3 bc3 τ3 a2 c2 τ3 bc2 τ3 + + + + 2 2 2 2L a2 c2 τ3 a2 c3 τ3 a2 c3 τ3 + + + 2L 2M 2M  −2c1 N + c1 c2 τ3 + c1 c3 τ3 , + 2 1 k2 = + (c3 M + l2 c2 τ2 2l2 K 2 M c1 c3 τ3 + 2c32 τ3 + 2c2 c3 τ3 bc3 τ2 + + 2 2 l2 c3 τ2 K 2 a2 L bc3 τ3 + + + N 2K 1 2 a2 c2 τ3 a2 c3 τ3 l2 c3 τ3 + + + 2 2L  2M a2 c3 τ3  2a2 c3 τ3 + (2l1 L(−a1 N − bN + N 2N a2 N a1 c2 τ3 bc2 τ3 a2 N + + + + L M 2 2 bc3 τ3 a2 c2 τ3 a1 c2 τ3 + + + 2 2 2L a2 c2 τ3 a2 c3 τ3 a2 c3 τ3 + + + 2L 2M 2M  −2c1 N + c1 c2 τ3 + c1 c3 τ3 , + 2 bc2 τ3 a1 c3 τ1 a2 N k3 = c3 N + + + M 2 2 L=

bc3 τ2 a2 c3 τ1 a2 c3 τ2 + bc3 τ3 + + 2 2 2M a2 c3 τ3 a2 c3 τ3 + + 2L 2M 1 3a2 c3 τ3 + (2c2 c3 τ3 + 2L 2   a2 N 2 −a1 N − bN + + 2c3 τ3 L a2 N a1 c2 τ3 bc2 τ3 + + + M 2 2 bc3 τ3 a2 c2 τ3 a1 c2 τ3 + + + 2 2 2L a2 c2 τ3 a2 c3 τ3 a2 c3 τ3 + + + 2L 2M 2M  −2c1 N + c1 c2 τ3 + c1 c3 τ3 , + 2 k4 = k5 = k6 = 1/ (−a1 N − bN a2 N a2 N + + L M bc2 τ3 a1 c2 τ3 a1 c2 τ3 + + + 2 2 2 a2 c2 τ3 a2 c3 τ3 bc3 τ3 + + + 2 2L 2L a2 c3 τ3 a2 c2 τ3 + + 2M 2M  −2c1 N + c1 c2 τ3 + c1 c3 τ3 + , 2  2 l1 2l1 τ1 Λ1 = k1 − ∗+ K1 x K1  l2 a2 z ∗ τ2 + a1 τ1 − a2 y ∗ τ1 − a2 z ∗ τ1 − k2 K2   c 2 + k3 2c22 τ3 + 2c2 c3 τ3 − c1 c2 τ3 + k4 ∗ x a1 c2 τ1 a2 c2 y ∗ τ1 a2 c2 z ∗ τ1 − + + 2 2 2  c3 τ3 l1 − c2 τ1 + − K1 2  c2 τ1 c3 τ1 c1 τ1 − + + 2 2 2    ∗ a2 z τ1 l1l2 τ2 + + k5 + x∗ K1 K2 2 2  l1 ∗ ∗ (2a2 τ1 z − bτ1 + a2 x τ1 − 2K 1  c2 l2 τ2 a2 c3 3z ∗ τ2 + k6 a2 c2 z ∗ τ3 + − 2 2K 2 +

123

S. Kundu, S. Maitra

 bc2 τ3 l2 c2 τ3 a2 c2 x ∗ τ3 − − , 2 2 2K 2  −2 −l1 a2 z ∗ τ1 l2 − ∗ Λ2 = k1 + k2 K1 K2 y  2l2 τ2 + + bτ2 − a2 τ2 z ∗ − a2 2τ2 x ∗ K2   + k3 2c32 τ3 + 2c2 c3 τ3 − c1 c3 τ3 +

 +k4

a2 c2 z ∗ τ1 a3 c3 z ∗ τ3 c3l1 τ1 + − 2 2K 1 2

 l1 c3 τ3 a1 c3 τ3 a2 c3 y ∗ τ3 a2 c3 z ∗ τ3 − − + 2 2K 1 2 2    ∗ a2 z τ1 + τ2 l1 l2 + k5 + y∗ K1 K2 2    l2 − a2 τ2 z ∗ − a1 τ2 − a2 τ2 y ∗ − a1 τ2 z ∗ 2K 2  c3 a2 c3 τ2 z ∗ + k6 ∗ + y 2 ∗ bc3 τ2 a2 c3 τ2 x − + 2 2  l2 − (c2 τ3 + 4c3 τ2 + c2 τ2 − c1 τ2 ) , 2K 2 l1 l2 (a1 τ1 − a2 y ∗ τ1 ) + k2 (bτ2 − a2 x ∗ τ2 ) Λ3 = k1 K1 K2   2c1 + k3 − ∗ − c1 c2 τ3 − c1 c3 τ3 z  a2 y ∗ a1 c1 l1 τ1 a1 c2 τ1 + k4 − ∗ + − z∗ z 2K 1 2  ∗ ∗ a2 c3 y τ3 a1 c3 τ3 a2 c2 y τ1 + − + 2 2 2  l1 + k5 (bτ1 − a2 τ1 x ∗ ) 2K 1  l2 (a1 τ2 − a2 τ2 y ∗ ) + 2K 2  a2 x ∗ − b a2 c2 τ3 x ∗ + k6 + z∗ 2  ∗ a2 c3 τ2 x c1l2 τ2 bc3 τ2 bc2 τ3 + + . − − 2 2 2K 2 2 +

References

2.

3.

4.

5.

6.

7.

8.

9.

10. 11.

12. 13.

14.

15.

16.

17. 18. 19.

20. 21.

1. Alfred, J.: Lotka, elements of physical biology (baltimore: Williams and wilkins,); vito volterra. Variazioni e flut-

123

tuazioni del numero dindividui in specie animali conviventi (1925) Arino, J., Wang, L., Wolkowicz, G.S.: An alternative formulation for a delayed logistic equation. J. Theor. Biol. 241(1), 109–119 (2006) Begon, M., Harper, J.L., Townsend, C.R., et al.: Ecology. Individuals, Populations and Communities. Blackwell, New York (1986) Castro, R., Sierra, W., Stange, E.: Bifurcations in a predator– prey model with general logistic growth and exponential fading memory. Appl. Math. Model. 45, 134–147 (2017) Chen, F.: The permanence and global attractivity of Lotka– Volterra competition system with feedback controls. Nonlinear Anal. Real World Appl. 7(1), 133–143 (2006) Chen, F.: Permanence of a discrete n-species food-chain system with time delays. Appl. Math. Comput. 185(1), 719–726 (2007) Chen, F., You, M.: Permanence for an integrodifferential model of mutualism. Appl. Math. Comput. 186(1), 30–34 (2007) Choudhury, S.R.: On bifurcations and chaos in predator– prey models with delay. Chaos Solitons Fractals 2(4), 393– 409 (1992) Dhar, J., Jatav, K.S.: Mathematical analysis of a delayed stage-structured predator–prey model with impulsive diffusion between two predators territories. Ecol. Complex. 16, 59–67 (2013) Elettreby, M.: Two-prey one-predator model. Chaos Solitons Fractals 39(5), 2018–2027 (2009) Fan, Y.H., Li, W.T.: Permanence for a delayed discrete ratiodependent predator–prey system with holling type functional response. J. Math. Anal. Appl. 299(2), 357–374 (2004) Goh, B.: Stability in models of mutualism. Am. Nat. 113(2), 261–275 (1979) Greenhalgh, D., Khan, Q.J., Pettigrew, J.S.: An ecoepidemiological predator–prey model where predators distinguish between susceptible and infected prey. Math. Methods Appl. Sci. 40(1), 146–166 (2017) Hou, Z.: On permanence of Lotka–Volterra systems with delays and variable intrinsic growth rates. Nonlinear Anal. Real World Appl. 14(2), 960–975 (2013) Hu, H., Teng, Z., Jiang, H.: On the permanence in nonautonomous Lotka–Volterra competitive system with puredelays and feedback controls. Nonlinear Anal. Real World Appl. 10(3), 1803–1815 (2009) Hugie, D.M.: Applications of Evolutionary Game Theory to the Study of Predator–Prey Interactions. Simon Fraser University, Burnaby (1999) Hutchinson, G.E.: Circular causal systems in ecology. Ann. N. Y. Acad. Sci. 50(4), 221–246 (1948) Keshet, E.L.: Mathematical Models in Biology. McGrawHill, New York (1988) Kingsland, S.: The refractory model: the logistic curve and the history of population ecology. Q. Rev. Biol. 57(1), 29–52 (1982) Kot, M.: Elements of Mathematical Ecology. Cambridge University Press, Cambridge (2001) Kundu, S., Maitra, S.: Stability and delay in a three species predator–prey system. In: AIP Conference Proceedings, vol. 1751, p. 020004. AIP Publishing (2016)

Dynamical behaviour of a delayed three species predator–prey model 22. Kuniya, T., Nakata, Y.: Permanence and extinction for a nonautonomous seirs epidemic model. Appl. Math. Comput. 218(18), 9321–9331 (2012) 23. Li, C.H., Tsai, C.C., Yang, S.Y.: Analysis of the permanence of an sir epidemic model with logistic process and distributed time delay. Commun. Nonlinear Sci. Numer. Simul. 17(9), 3696–3707 (2012) 24. Liao, X., Zhou, S., Chen, Y.: Permanence and global stability in a discrete n-species competition system with feedback controls. Nonlinear Ana. Real World Appl. 9(4), 1661–1671 (2008) 25. Liu, M., Bai, C.: Optimal harvesting of a stochastic mutualism model with Lévy jumps. Appl. Math. Comput. 276, 301–309 (2016) 26. Liu, S., Chen, L.: Necessary-sufficient conditions for permanence and extinction in Lotka–Volterra system with distributed delays. Appl. Math. Lett. 16(6), 911–917 (2003) 27. Lynch, S.: Dynamical Systems with Applications Using MATLAB. Springer, Berlin (2004) 28. Ma, J., Zhang, Q., Gao, Q.: Stability of a three-species symbiosis model with delays. Nonlinear Dyn. 67(1), 567–572 (2012) 29. MacLean, M., Willard, A.: The logistic curve applied to Canada’s population. Can. J. Econ. Polit. Sci. 3(02), 241– 248(1937)

30. Muroya, Y.: Permanence and global stability in a Lotka– Volterra predator–prey system with delays. Appl. Math. Lett. 16(8), 1245–1250 (2003) 31. Pearl, R., Reed, L.J.: The logistic curve and the census count of i930. Science 72(1868), 399–401 (1930) 32. Saha, T., Bandyopadhyay, M.: Dynamical analysis of a delayed ratio-dependent prey–predator model within fluctuating environment. Appl. Math. Comput. 196(1), 458–478 (2008) 33. Teng, Z., Zhang, Y., Gao, S.: Permanence criteria for general delayed discrete nonautonomous n-species kolmogorov systems and its applications. Comput. Math. Appl. 59(2), 812–828 (2010) 34. Verhulst, P.F.: Notice sur la loi que la population suit dans son accroissement. Correspondance mathématique et physique publiée par a. Quetelet 10, 113–121 (1838) 35. Xu, C., Wu, Y.: Dynamics in a Lotka–Volterra predator–prey model with time-varying delays. In: McKibben, M. (ed.) Abstract and Applied Analysis, vol. 2013. Hindawi Publishing Corporation, Cairo, Egypt (2013). https://doi.org/10. 1155/2013/956703 36. Zhao, J., Jiang, J.: Average conditions for permanence and extinction in nonautonomous Lotka–Volterra system. J. Math. Anal. Appl. 299(2), 663–675 (2004)

123