A three species dynamical system involving prey–predation, competition and commensalism

Applied Mathematics and Computation 273 (2016) 54–67

Contents lists available at ScienceDirect

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

A three species dynamical system involving prey–predation, competition and commensalism Sunita Gakkhar, Komal Gupta∗ Department of Mathematics, Indian Institute of Technology, Roorkee, Uttarakhand 247667, India

a r t i c l e

i n f o

Keywords: Holling-type II response Commensalism Refuge Persistence Hopf bifurcation

a b s t r a c t In this paper, a three species dynamical system is explored. The system consisting of two logistically growing competing species and the third species acts as a predator as well as host. It is predating over second species with Holling type II functional response, while first species is benefited from the third species. In addition, the prey species move into a refuge to avoid high predation. The essential mathematical features of the proposed model are studied in terms of boundedness, persistence, local stability and bifurcation. The existence of transcritical bifurcations have been established about two axial points. It has been observed that survival of all three species may be possible due to commensalism. Numerical simulations have been performed to show the Hopf bifurcation about interior equilibrium point. The existence of period-2 solution is observed. Further, the bifurcations of codimension-2 have also been investigated. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The study of ecological models has created much interest among authors. Most of the two species models consider only one type of interaction at a time: prey–predation, competition or mutualism. The prey–predator models with Holling types I, II, III and IV are discussed in [1–6]. The Leslie–Gower, ratio-dependent and Beddington–DeAngelis functional responses are also extensively investigated in [7–12]. To avoid high risk of predation, prey species defend themselves by taking refuge. The term refuge can be generally defined to include any strategy that decreases predation risk. Many authors have studied refuge of prey species in [13–15]. Lot of work is available on two species competing systems but limited work is available on mutualism [16–19]. The stability, extinction, co-existence, bifurcations and limit cycles are the typical behaviors investigated for such nonlinear two dimensional systems. These models are extended for three species food chain and food-web systems [20–29]. Food-web systems have been studied with various combinations of prey–predation and competition in [30,31]. These three species systems exhibit rich dynamical behaviors as compared to their counter parts in two species systems. The existence of chaos, strange attractors and quasi-periodic solutions may be possible in food-webs and food chains. Few food-webs are also investigated considering mutualism in [32,33]. Commensalism is a form of mutualism between individuals of two species in which one is benefited and other is unaffected. The commensal may obtain food, shelter or support from the host species. Kumar and Pattabhiramacharyulu [34] have studied the three species system in combination with prey–predation and commensalism. One or more species taking refuge to avoid predation are considered by many authors in food-web models in [23,25,35]. ∗

Corresponding author. Tel.: +91 9012995115. E-mail addresses: [email protected] (S. Gakkhar), [email protected] (K. Gupta).

http://dx.doi.org/10.1016/j.amc.2015.09.036 0096-3003/© 2015 Elsevier Inc. All rights reserved.

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

55

Fig. 1. Schematic diagram of three species food-web system.

In this paper, a three species food-web system consisting of prey–predation, competition and commensalism is considered. The refuge is also incorporated in the prey species. 2. Model formulation Consider a three species food-web system, whose interactions are shown in schematic diagram (Fig. 1). Let X and Y are the two logistically growing competing species with intrinsic growth rates and carrying capacities ri and Ki (i = 1, 2), respectively. The species Z is predating over Y with Holling type II response. The species X is assumed to be benefited from the presence of Z, as such X is commensal of Z. The proportion p of species Y is refuged from predation. The coefficients αi j (i, j = 1, 2; i = j) and δ 13 are the inter-species competition coefficients and commensal coefficient, respectively. Accordingly, the dynamics of the system is governed by the following system of non-linear differential equations:



dX = r1 X 1 − dT  dY = r2 Y 1 − dT





X α12Y − + δ13 XZ K1 K1  Y a(1 − p)Y Z α21 X − − K2 K2 b + (1 − p)Y

ac(1 − p)Y Z dZ = Z −e + dT b + (1 − p)Y

(2.1)



The non-negative initial conditions are:

X (0) ≥ 0, Y (0) ≥ 0, Z (0) ≥ 0 The following non-dimensional variables and parameters are introduced:

X Y aZ b e , y= , z= , w1 = , w2 = , K1 K2 r1 K2 K2 r1 ac r2 α12 K2 α21 K1 δ13 K2 , r = , β12 = , β21 = , δ= w3 = r1 r1 K1 K2 a t = r1 T, x =

Accordingly, the non-dimensionalized equations and associated initial conditions are:

dx = x[1 − x − β12 y] + δ xz = x f1 (x, y, z) = F1 (x, y, z) dt dy (1 − p)yz = y f2 (x, y, z) = F2 (x, y, z) = ry[1 − y − β21 x] − dt w1 + (1 − p)y



w3 (1 − p)yz dz = z −w2 + dt w1 + (1 − p)y

(2.2)



= z f3 (x, y, z) = F3 (x, y, z)

x(0) ≥ 0, y(0) ≥ 0, z(0) ≥ 0. The system can be written in the matrix form as

X˙ = F (X ); X (0) = X0 ; X = (x, y, z)T ∈ R3 and X0 ∈ R3+

(2.3)

2.1. Preliminaries results In this section, some important results about existence, uniqueness and boundedness of the system (2.2) are presented. The interaction functions Fi (i = 1, 2, 3) are continuous and have continuous partial derivatives in the state space R3+ . Hence, the solution of the system with non-negative initial conditions exists and it is unique. The positiveness of the solution is guaranteed by Theorem (2.1.1). Theorem 2.1.1. The system (2.3) has non-negative solutions for non-negative initial conditions.

56

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

Proof. The result is a direct consequence of Nagumo’s theorem [36].  In ecology, boundedness of a system implies that the system is biologically well behaved. Theorem (2.1.2) ensures the boundedness of the system (2.2). The survival domain of species z is also obtained in Theorem (2.1.3). Theorem 2.1.2. The solution of the system (2.2) is uniformly bounded. Proof. From the second equation of (2.2)

dy ≤ ry(1 − y), y0 ≥ 0 dt The usual comparison theorem yields

y(t ) ≤

y0 , y0 + (1 − y0 )e−rt

∀t ≥0

⇒ 0 ≤ y(t ) ≤ 1, as t → ∞

(2.4)

Consider

W1 (t ) = y(t ) +

1 z(t ), W1 (0) = W10 ≥ 0 w3

(2.5)

then

w2 dW1 ≤ ry(1 − y) − z dt w3 After simplifying using max[0,1] y(1 − y) =

r

dW1 + w2W1 ≤ + w2 dt 4 Therefore,

W1 (t ) ≤

r 4

+ w2

⇒ 0 < W1 <

 1

r

w2

1 4,

we get

(1 − e−w2 t ) + W10 e−w2 t ∀ t ≥ 0

+ w2

4



 1

as t → ∞

w2

(2.6)

Again consider

W2 (t ) = x(t ) + y(t ) +

1 z(t ); W2 (0) = W20 ≥ 0 w3

(2.7)

then

w2 dW2 ≤ x(1 − x) + δ xz + ry(1 − y) − z dt w3 Using max[0,1] y(1 − y) =



1 and (2.6) 4



dW2 r + w2W2 ≤ + w2 + f (x) dt 4



where f (x) = x(1 − x) + w2 x + δ w3 x 1 + maximum value of f(x) is

W2 (t ) ≤ 1 +

(r + a2 ) 

.



r 4w2





4w2

⇒ 0 < W2 < 1 +

4w2

a2 where a is given as a = 1 + w2 + δ w3 a + 4

Accordingly, from (2.8)



(2.8)

 r

(1 − e−w2 t ) + W20 e−w2 t ∀ t ≥ 0

r + a2 4w2

 as t → ∞

Hence, the solution of the system is uniformly bounded. Further, the system is dissipative.  Theorem 2.1.3. The necessary condition for the survival of the species z is

0 <

w1 w2

(1 − p)(w3 − w2 )

< 1

(2.9)

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

57

Proof. From Theorem (2.1.2), y(t) ≤ 1. Then for the species z



dz w3 (1 − p) ≤z − w2 dt w1 + (1 − p)



< z[(1 − p)(w3 − w2 ) − w1 w2 ] Accordingly z(t ) < z0 e[(1−p)(w3 −w2 )−w1 w2 ]t If w3 ≤ w2 , then limt→∞ z(t ) = 0 And if w3 > w2 , with (1 − p)(w3 − w2 ) − w1 w2 < 0, then limt→∞ z(t ) = 0 w1 w2 Accordingly, the necessary condition for z to survive is 0 < < 1  (1 − p)(w3 − w2 ) 3. Analysis In the absence of species z, the system (2.2) reduces to well-studied Lotka–Volterra competition model. In the absence of y, species x will be logistically growing while z will go to extinction in absence of food. 3.1. Kolmogorov analysis In the absence of species x, which is commensal of z species, the system reduces to Kolmogorov’s system:



dy (1 − p)z = y r(1 − y) − dt w1 + (1 − p)y



w3 (1 − p)y dz = z −w2 + dt w1 + (1 − p)y





(3.1)

under the following condition:

w2 <

w3 (1 − p) [w1 + (1 − p)]

(3.2)

Local stability analysis for the Kolmogorov system (3.1) gives the following results: ¯ = (0, 0) always exists and it is a saddle point. (i) The trivial equilibrium point E21 ¯ = (1, 0) always exists and it is a stable point if (ii) The equilibrium point E22

(1 − p)w3 [w1 + (1 − p)]

w2 >

(3.3)

¯ will always be saddle. Otherwise, it will be a saddle point. It is clear that when system is Kolmogorov under (3.2), E22 ¯ = (y, ¯ z¯) is given by: (iii) The nontrivial positive equilibrium point E23

y¯ =

w1 w2

(1 − p)(w3 − w2 )

, z¯ =

r(1 − y¯ )(w1 + (1 − p)y¯ ) (1 − p)

The associated Jacobian is

⎡

y¯ z¯ (w1 + (1 − p)y¯)2 w1 w3 (1 − p)z¯ (w1 + (1 − p)y¯)2

−ry¯ + (1 − p)2

J =

⎢ ⎣

−

(3.4)

⎤ (1 − p)y¯ w1 + (1 − p)y¯ ⎥ ⎦ 0

¯ is The necessary and sufficient condition for linear stability of E23

w2 >

w3 [(1 − p) − w1 ] (1 − p) + w1

(3.5)

It is observed that the above condition will be trivially satisfied when p > 1 − w1 . Combining with the Kolmogorov condi¯ is stable equilibrium point if tion (3.2) E23

w3 [(1 − p) − w1 ] < w2 < (1 − p) + w1

(1 − p)w3 (1 − p) + w1

(3.6)

Eq. (3.6) can be written in terms of refuge parameter as follows:

1−

w1 (w3 + w2 ) w1 w2 < p < 1− w3 − w2 w3 − w2

(3.7)

58

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

Periodic solutions iny–z plane are possible when

w3 [(1 − p) − w1 ] (1 − p) + w1

w2 <

(3.8)

¯ occurs at The Hopf point about E23

w2 =

w3 [(1 − p) − w1 ] (1 − p) + w1

(3.9)

threshold parametric values of other parameters for Hopf bifurcation are:

w3 =

w2 [(1 − p) + w1 ] w1 (w3 + w2 ) , p=1− w3 − w2 (1 − p) − w1

(3.10)

¯ depends upon prey refuge p. The Hopf bifurcation point with respect to E23 3.2. Local stability analysis The local stability analysis of biologically feasible equilibrium points are discussed in this section. The Jacobian matrix for an equilibrium point E = (x, y, z) is given as:

⎡1 − 2x − β y + δ z 12 ⎢ ⎢ −rβ21 y JE = ⎢ ⎣

−β12 x

w3 (1 − p)z [w1 + (1 − p)y]2 w1 w3 (1 − p)z [w1 + (1 − p)y]2

r(1 − y − β21 x) − ry −

0

δz

⎤

⎥ −(1 − p)y ⎥ w1 + (1 − p)y ⎥ ⎦ w3 (1 − p)y − w2 w1 + (1 − p)y

(1) The trivial equilibrium point E0 = (0, 0, 0) always exists. The eigenvalues at E0 are 1, r and −w2 , hence origin is a saddle point having stable manifold in z-direction and unstable manifold in x–y plane. (2) The axial equilibrium point E1 = (1, 0, 0) always exists. The eigenvalues at E1 are −1, r(1 − β21 ) and −w2 . Accordingly, the point E1 is locally asymptotically stable if

β21 > 1

(3.11)

otherwise it will be a saddle point having unstable manifold in y-direction. Further, the system admits transcritical bifurcation corresponding to zero eigenvalue if

β21 = 1 (3) The second axial equilibrium point E2 = (0, 1, 0) always exists. The eigenvalues at E2 are (1 − β12 ), −r and −w2 + w3 (1 − p) . The point E2 is locally asymptotically stable provided w1 + (1 − p)

β12 > 1 and p > 1 −

w1 w2 w3 − w2

(3.12)

E2 becomes saddle point if one of the above conditions are violated. Further, the system undergoes a transcritical bifurcation if

β12 = 1 or p = 1 −

w1 w2 w3 − w2

(3.13)

ˆ y, ˆ 0) in x–y plane is given as: (4) The boundary equilibrium point E12 = (x,

xˆ =

1 − β21 1 − β12 , yˆ = 1 − β12 β21 1 − β12 β21

The equilibrium point E12 exists if one of the following conditions are satisfied:

β12 < 1, β21 < 1 β12 > 1, β21 > 1

(3.14a) (3.14b)

The eigenvalues of JE12 are obtained as:

  (1 − β12 ) + r(1 − β21 ) λˆ1 + λˆ2 =− (1 − β12 β21 ) λˆ1 λˆ2 =

r(1 − β12 )(1 − β21 ) (1 − β12 β21 )

(3.15a) (3.15b)

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

λˆ3 =

w3 (1 − p)yˆ − w2 w1 + (1 − p)yˆ

59

(3.15c)

Clearly, E12 is locally asymptotically stable when

β12 < 1, β21 < 1 and

w1 w2

(1 − p)(w3 − w2 )

>

(1 − β21 ) (1 − β12 β21 )

(3.16)

Note that, when E12 is stable then no axial point is stable. Also, when axial points E1 and E2 are stable, then system admits an unstable boundary equilibrium point E12 . ¯ z¯), where y¯ and z¯ are given by (3.4), exists in the y–z plane under Kol(5) The nonnegative equilibrium point E23 = (0, y, mogorov’s condition (3.2). The eigenvalues of JE23 are given as follows:

λ¯1 =1 − β12 y¯ + δ z¯

(3.17a)

(1 − p)2 y¯z¯ ¯2 [w1 + (1 − p)y] w w (1 − p)2 y¯ z¯ λ¯2 λ¯3 = 1 3 >0 ¯3 [w1 + (1 − p)y] λ¯2 + λ¯3 =−ry¯ +

(3.17b) (3.17c)

¯ . However, it is stable when λ¯1 is negative. Clearly, from the Eq. (3.17b) E23 has same stability behavior in y–z plane as E23 It has unstable manifold in the x-direction orthogonal to the y–z plane when λ¯1 is positive. (6) The interior equilibrium point E ∗ = (x∗ , y∗ , z∗ ) is obtained as:

x∗ =

w1 w2 (1 − p)(w3 − w2 )(1 − β12 y∗ ) + rδ w1 w3 (1 − y∗ ) , y∗ = (1 − p)(w3 − w2 ) + rδβ21 w1 w3 (1 − p)(w3 − w2 )

z∗ =

rw1 w3 [(1 − β21 ) − y∗ (1 − β12 β21 )] (1 − p)(w3 − w2 ) + rδβ21 w1 w3

(3.18)

E∗ is feasible under the following parametric restrictions:

w1 w2 , w3 − w2

p<1−

w1 w2

(1 − p)(w3 − w2 )

<

(1 − β21 ) (1 − β12 β21 )

(1 − p)(w3 − w2 )(1 − β12 y∗ ) + rδ w1 w3 (1 − y∗ ) > 0 ˆ y, ˆ 0) guarantees nonexistence of interior point E ∗ = (x∗ , y∗ , z∗ ). The Clearly, the stability of E2 = (0, 1, 0) and/or E12 = (x, Jacobian at E ∗ = (x∗ , y∗ , z∗ ) is



JE ∗ =

H11 H21 0

H12 H22 H32

H13 H23 0



where

H11 = −x∗ < 0 H12 = −β12 x∗ < 0 H13 = δ x∗ > 0 H21 = −rβ21 y∗ < 0

(1 − p)2 y∗ z∗ (w1 + (1 − p)y∗ )2 (1 − p)y∗ H23 = − <0 w1 + (1 − p)y∗ w1 w3 (1 − p)z∗ >0 H32 = (w1 + (1 − p)y∗ )2 H22 = −ry∗ +

The characteristic equation is λ3 + a0 λ2 + a1 λ + a2 = 0, with

a0 = −(H11 + H22 ) a1 = H11 H22 − H12 H21 − H23 H32 a2 = H11 H23 H32 − H13 H21 H32 According to Routh–Hurwitz criteria, E5 = (x∗ , y∗ , z∗ ) is locally asymptotically stable if

a0 > 0,

a2 > 0,

a0 a1 − a2 > 0

(3.19)

60

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

Further, the existence of Hopf bifurcation, generalized Hopf bifurcation, Bogdanov–Takens bifurcation etc. are possible provided

a0 a1 − a2 = 0

(3.20)

Since the expressions in (3.20) involve many parameters r, β12 , β21 , δ, p, w1 , w2 and w3 , its further analysis from bifurcation point of view is carried out in the Section 4 on numerical simulation. 3.3. Transcritical bifurcation In the following theorems existence of transcritical bifurcation has been established using Sotomayor’s theorem [37] about the two equilibrium points E1 = (1, 0, 0) and E2 = (0, 1, 0). Theorem 3.3.1. The dynamical system (2.2) undergoes a transcritical bifurcation at β21 = 1 about E1 = (1, 0, 0) when β 12 = 1. Proof. The Jacobian matrix about E1 is evaluated as



JE1 =

−β12 r(1 − β21 ) 0

−1 0 0



δ 0 −w2

Let V = ( − β12 , −1, 0)T and W = (0, 1, 0)T be the eigenvectors of JE1 and JET , respectively corresponding to zero eigenvalue. 1

Denoting F = [F1 F2 F3 ]T . Now the three conditions of Sotomayor’s theorem are computed below: ∗ )=0 (i) W T Fβ21 (E1 , β21 T ∗ )V ] = r = 0, where (ii) W [DFβ21 (E1 , β21



DFβ21 (E1 , β

∗ 21

) =

0 0 0

0 −r 0

3 (iii) Let us denote D2 F(V, V) as i=1  3j=1

0 0 0



∂ 2F ∗ )(V, V )] = 2r (β − 1) = 0 when β v v then W T [D2 F (E1 , β21 12 12 = 1 ∂ xi ∂ x j i j

∗ = 1. Accordingly, the system (2.2) admits transcritical bifurcation at β21 = β21



Theorem 3.3.2. The dynamical system (2.2) undergoes a transcritical bifurcation at β12 = 1 about E2 = (0, 1, 0) when β 21 = 1. Proof. The Jacobian matrix about E2 = (0, 1, 0) is computed as

JE2

⎡(1 − β ) 12 ⎢ ⎢ −rβ21 = ⎢ ⎣

0

⎥ 1− p ⎥ w1 + (1 − p) ⎥ ⎦ w3 (1 − p) − w2 w1 + (1 − p)

−r

0

⎤

0 −

0

Let V = ( − 1, −β21 , 0)T and W = (1, 0, 0)T be the eigenvectors of JE2 and JET , respectively corresponding to zero eigenvalue. 2

Denoting F = [F1 F2 F3 ]T . Now the three conditions of Sotomayor’s theorem are computed below: ∗ )=0 (i) W T Fβ12 (E2 , β12 T ∗ )V ] = −1 = 0, where (ii) W [DFβ12 (E2 , β12



DFβ12 (E2 , β

∗ 12

) =

−1 0 0

0 0 0



0 0 0

∗ )(V, V )] = 2(β − 1) = 0 when β (iii) W T [D2 F (E2 , β12 21 21 = 1 ∗ = 1.  Accordingly, the system (2.2) admits transcritical bifurcation at β12 = β12

It is noted that JE2 admits zero eigenvalue also at w2 = respect to w2 is stated in the following theorem:

w3 (1−p) . w1 +(1−p)

The existence of transcritical bifurcation about E2 with

Theorem 3.3.3. The dynamical system (2.2) undergoes a transcritical bifurcation at w2 = w∗2 about E2 = (0, 1, 0) where

w∗2 =

w3 (1 − p) w1 + (1 − p)

The proof is on similar lines and hence omitted from the text. Note 1:The expression (3.21) can be used to obtain the transcritical bifurcation with respect to w1 , w3 and p also. Note 2:The system (2.2) neither admits saddle- node nor pitchfork bifurcation about E1 and E2 .

(3.21)

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

61

3.4. Persistence of the system Biologically a system persists if all the populations of the system survive in future time. In order to establish the persistence of the system, the abstract theorem of Freedman and Waltman [38] is applied in the following: Theorem 3.4.1. The system (2.2) persists in the absence of nontrivial periodic solutions in the boundary planes x–y and y–z, provided:

ˆ y, ˆ 0) = λˆ3 = f3 (x,

w3 (1 − p)yˆ − w2 > 0 (w1 + (1 − p)yˆ)

(3.22a)

¯ z¯) = 1 − β12 y¯ + δ z¯ > 0 λ¯1 = f1 (0, y,

(3.22b)

Proof. In the following the hypotheses (1)–(4) are being established with respect to system (2.2):

∂ f2 ∂ f3 −1 w1 w3 = = < 0, >0 ∂z ∂y (w1 + y) (w1 + y)2 ∂ f2 (x, y, 0) 2. f2 (0, 0, 0) = r > 0, f2 (0, 1, 0) = 0 and = −r < 0 ∂y 1.

Accordingly, prey population y grows to its carrying capacity in the absence of predator z. Further, f3 (x, 0, z) = −w2 < 0, i.e. predator will not survive without prey. ˆ y, ˆ 0) 3. The x–y boundary plane admits an equilibrium point E12 = (x, 4. There is no equilibrium point in x–z boundary plane. ¯ z¯) such that f2 (0, y, ¯ z¯) = f3 (0, y, ¯ z¯) = 0. Further, However, the boundary plane y–z admits an equilibrium point E23 = (0, y, f3 (0, 1, 0) > 0, therefore E2 = (0, 1, 0) is unstable in z-direction. By the above assumptions, there is at most one equilibrium point interior to each positive coordinate plane. The conditions (3.22a) and (3.22b) ensure that planar equilibrium points are unstable in the direction orthogonal to the coordinate plane. Further, the system is bounded. It follows the proof.  The system (2.2) admits limit cycles in the y–z boundary plane provided the condition (3.8) is satisfied. In this case the following theorem for persistence is established. Theorem 3.4.2. Let the conditions (3.22a) and (3.22b) hold and there exists finite number of limit cycles (φ y (t), φ z (t)) in y–z plane with period T. Then the system (2.2) persists provided



T 0

(1 − β12 φy (t ) + δφz (t ))dt > 0

(3.23)

Proof. Let there exists a limit cycle in y–z plane. The variational matrix V(x(t), y(t), z(t)) about (0, φ y (t), φ z (t)) is given as:

⎡

1 − β12 φy (t ) + δφz (t )

⎢ ⎢

V (0, φy (t ), φz (t )) = ⎢ ⎢

−rβ21 φy (t )

⎣

0

(1 − p)2

0

0

φy (t )φz (t ) − rφy (t ) (w1 + (1 − p)φy (t ))2 (1 − p)w1 w3 φz (t ) (w1 + (1 − p)φy (t ))2

−

⎤

(1 − p)φy (t ) ⎥ ⎥ w1 + (1 − p)φy (t ) ⎥ ⎥ ⎦ 0

Consider solution of the system with non-negative initial conditions (α 1 , α 2 , α 3 ) sufficiently closed to the limit cycle. From the above matrix we get

dx = (1 − β12 φy (t ) + δφz (t ))x dt

(3.24)

and its solution is given as



x(t ) = exp

T 0

 (1 − β12 φy (t ) + δφz (t ))dt α1

(3.25)

partially differentiate (3.25) with respect to α 1

 T  ∂x = exp (1 − β12 φy (t ) + δφz (t ))dt ∂α1 0

By Taylor’s expansion theorem (about α1 = 0)

x(t, α1 , α2 , α3 ) − x(t, 0, α2 , α3 )



∼ = exp

0

T

 (1 − β12 φy (t ) + δφz (t ))dt α1 .

T Thus x increases or decreases as 0 (1 − β12 φy (t ) + δφz (t ))dt is positive or negative, respectively. Since the equilibrium point ¯ z¯) and these limit cycles are the only possible attractors in y–z plane with positive initial conditions, the trajectories E23 = (0, y, go away from y–z plane if (3.23) hold.

62

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

0.25

z

0.2 0.15 0.1 0.05 0.8

1 0.6

0.9 0.8

0.4 0.7

0.2 y

0.6 0

0.5

x ∗

Fig. 2. Phase plane diagram showing local stability of E : (0.84044, 0.41514, 0.16238).

Fig. 3. (a) Bifurcation diagram with respect to w1 (b) phase plane showing limit cycle for the data set (4.1) at w1 = 0.11.

Thus, the proof is complete.  It may be observed from (3.4.1) that survival of three species depend on prey–predation, competition (β 12 ), commensalism (δ ) and refuge (p). The species x may not survive in the absence of commensalism due to strong competition (when β 12 > 1) from y. However, its survival may be possible for suitably large value of δ . The impact of refuge on the survival may be more complex. 4. Numerical simulation In this section, some numerical results are given to discuss the dynamical aspects of system (2.2) and to verify analytical findings. Consider the following set of parameters:

w1 = 0.14, w2 = 0.32009, w3 = 0.5, p = 0.4, r = 0.5, δ = 0.04, β12 = 0.4 and β21 = 0.1 The positive equilibrium point

E∗

(4.1)

is found as (0.84044, 0.41514, 0.16238). The three eigenvalues −0.849670 and −0.019512 ±

(0.16238)i have negative real parts. Accordingly, the interior point E∗ is locally asymptotically stable. Further, the phase plane

diagram, drawn in Fig. 2 verifies the local stability of E∗ . The bifurcation diagram is drawn with respect to parameter w1 in Fig. 3(a), using XPPAUT [39]. The diagram has two different lines and two different circles. Thick lines show line of stable fixed points, stable periodics are solid circles, unstable fixed points are thin lines and unstable periodics are open circles. It is observed that Hopf bifurcation occurs at w1 ∼ = 0.12323 and the limit cycle is drawn in Fig. 3(b) for w1 = 0.11. Further, the Hopf bifurcation occurs at w2 ∼ = 0.30047, w3 ∼ = 0.53268 and p∼ =0.31821 and diagrams with respect to w2 , w3 and p are plotted in Fig. 4. In case of p = 0.31, the point E∗ will be unstable and oscillating behavior is shown in Fig. 5. When the amount of refuge is p = 0.33, then Fig. 6 shows the stability of the equilibrium state. Now consider the following set of data:

w1 = 0.2001, w2 = 0.32009, w3 = 0.5, p = 0.4, r = 0.5, δ = 0.04, β12 = 1.5 and β21 = 0.49

(4.2)

The Hopf bifurcation appears to occur at w1 ∼ = 0.1815, w2 ∼ = 0.30678, w3 ∼ = 0.5228 and p∼ =0.33856. The period-2 solution is obtained for w1 in the interval (0.1195, 0.13818) and p∼ The bifurcation diagrams with respect to parameters w1 , w2 , w3 =0.13119. and p are plotted in Fig. 7. Further, the time series and phase plane diagram shows period-2 solution in Figs. 8 and 9, respectively.

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

Fig. 4. (a)–(c): Bifurcation diagrams of y with respect to w2 , w3 and p for the set (4.1), showing maxima and minima of limit orbits.

Fig. 5. Oscillating behavior for the set (4.1) at p = 0.31.

Fig. 6. Time series showing asymptotically stability of E∗ for the set (4.1) at p = 0.33.

63

64

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

Fig. 7. Bifurcation diagrams of y with respect to parameters w1 , w2 , w3 and p showing maxima and minima of limit orbits.

Fig. 8. Time series for the data (4.2).

To see the effect of commensalism on coexistence of three species, consider the following choice of data:

w1 = 0.4, w2 = 0.32009, w3 = 0.5, r = 0.5, β12 = 1.5, β21 = 0.4 with p = 0 and δ = 0

(4.3)

It may be noted that the condition (3.22b) is not satisfied giving extinction of x. However, when δ is increased to 0.5, the condition (3.22b) is satisfied leading to coexistence of all the three species. The Fig. 10 clearly shows the coexistence of three species in the form of stability of E ∗ = (0.01137, 0.71167, 0.157737). It may be concluded that commensalism may be responsible for the coexistence. To further explore bifurcations in the system, if any, the software package MATCONT is used. Consider the following data set:

β12 = 1.5, β21 = 0.4, r = 0.5, δ = 0.5, p = 0, w1 = 0.2001, w2 = 0.32009, w3 = 0.5

(4.4)

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

65

0.25 0.2

z

0.15 0.1 0.05 0 1 0.5 y

0

0

0.2

0.4

0.6

0.8

1

x

Fig. 9. Phase plane diagram for w1 = 0.12.

Fig. 10. The figure shows coexistence of three species at δ = 0.5 and p = 0 for (4.3).

The Hopf bifurcation occurs at w1 = 0.28455632 with first Lyapunov coefficient as −1.019195e − 001. The generalized Hopf bifurcation occurs when the values of w1 and w3 are taken as 0.048203494 and 0.34427351, respectively with second Lyapunov coefficient as −1.234255e + 000. No Bogdanov–Takens bifurcation has been observed in positive domain for this data set, see Fig. 11(a). Now, consider another data set:

β12 = 1.5, β21 = 0.4, r = 0.5, δ = 0.6, p = 0.2, w1 = 0.2001, w2 = 0.32009, w3 = 0.5

(4.5)

For this, the Hopf bifurcation occurs at w1 = 0.22793998 with first Lyapunov coefficient as −9.159184e − 002. The generalized Hopf bifurcation occurs in (w1 , w3 ) plane at (0.04553632, 0.34886613) with second Lyapunov coefficient as −7.195912e − 001. Here, again no Bogdanov–Takens bifurcation has been observed, see Fig. 11(b). Now considering the data set (4.5) with no prey refuge, the Hopf bifurcation occurs when w1 = 0.28494025 with first Lyapunov coefficient as −9.159621e − 002. The generalized Hopf bifurcation is obtained in (w1 , w3 ) plane at (0.056918101, 0.34886489) with second Lyapunov coefficient as −7.196872e−001. At these values, the interior point will become E ∗ = (0.1165629, 0.6331532, 0.11048788) and corresponding eigenvalues are −0.286234, ± (0.0431737)i. The Bogdanov–Takens bifurcation occurs in (w1 , w3 ) plane for (0.00018736171, 0.3201749), then the interior point will be E∗ ∼ =(0.0024, 0.7063, 0.103383) and eigenvalues will be −0.20936, ± (0.00039)i. The normal form coefficients are (a,b)=(−1.246728e−007, 1.631720e−002). The bifurcation diagram in (w1 , w3 ) plane is drawn in Fig. 11(c). 5. Discussion In this paper, we have considered a three species system in which two competing species are logistically growing. One of the competitors is being predated by third species with Holling type II functional response while the other competitor is in commensal relationship with the third species acting as its host. The prey species take refuge to avoid predation. The stability analysis has been carried out for all possible feasible equilibrium points. Also, sufficient conditions for persistence of the system have been

66

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67

a

b 0.5

0.5 H

H 0.45 w

3

w3

0.45

0.4

0.4

GH

GH

0.35

0.35 BT

BT 0

0.05

0.1

0.15

0.2

0.25

0.3

0

0.05

0.1

0.15

0.2

w

w1

1

c 0.5 H

w

3

0.45

0.4 GH 0.35 BT 0

0.05

0.1

0.15 w

0.2

0.25

0.3

1

Fig. 11. (a–c) The codimension-2 bifurcation diagrams: (a) for the data set (4.4), (b) for the data set (4.5) with p = 0, (c) for the data set (4.5) with p = 0.

established. The existence of stable limit cycle and period-2 solution have been obtained. Also, the food-web system exhibits the generalized Hopf bifurcation but Bogdanov–Takens bifurcation point is detected only for a certain value of commensal coefficient. Through numerical simulation, commensalism has been found to be responsible for coexistence. The stable effect of the refuge parameter has also been shown. Acknowledgment The authors are thankful to reviewers for their valuable suggestions. The second author (K. Gupta) would like to thank “University Grants Commission (UGC)” for providing Junior Research Fellowship through grant no. 6405-11-044. References [1] F. Chen, Z. Ma, H. Zhang, Global asymptotical stability of the positive equilibrium of the Lotka–Volterra prey–predator model incorporating a constant number of prey refuges, Nonlinear Anal. Real World Appl. 13 (2012) 2790–2793. [2] G. Seo, D.L. DeAngelis, A predator–prey model with a Holling type i functional response including a predator mutual interference, J. Nonlinear Sci. 21 (2011) 811–833. [3] G. Peng, Y. Jiang, C. Li, Bifurcations of a Holling-type ii predator–prey system with constant rate harvesting, Int. J. Bifurc. Chaos 19 (2009) 2499–2514. [4] S. Yu, Global asymptotic stability of a predator–prey model with modified Leslie–Gower and Holling-type II schemes, Discret. Dyn. Nat. Soc. 2012 (2012) 8. [5] 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 (2006) 672–683. [6] R.K. Naji, R.N. Shalan, The dynamics of Holling type IV prey–predator model with intra-specific competitions, Iraqi J. Sci. 54 (2013) 386–396. [7] M. Aziz-Alaoui, M.D. Okiye, Boundedness and global stability for a predator–prey model with modified Leslie–Gower and Holling-type II schemes, Appl. Math. Lett. 16 (2003) 1069–1075. [8] S. Yu, Global stability of a modified Leslie–Gower model with Beddington–DeAngelis functional response, Adv. Differ. Equ. 2014 (2014) 1–14. [9] D. Xiao, W. Li, M. Han, Dynamics in a ratio-dependent predator–prey model with predator harvesting, J. Math. Anal. Appl. 324 (2006) 14–29. [10] M. Sen, M. Banerjee, A. Morozov, Bifurcation analysis of a ratio-dependent prey–predator model with the Allee effect, Ecol. Complex. 11 (2012) 12–27. [11] S.-B. Hsu, T.-W. Hwang, Y. Kuang, Global analysis of the Michaelis–Menten-type ratio-dependent predator–prey system, J. Math. Biol. 42 (2001) 489–506. [12] R.S. Cantrell, C. Cosner, On the dynamics of predator–prey models with the Beddington–DeAngelis functional response, J. Math. Anal. Appl. 257 (2001) 206–222. [13] F. Chen, L. Chen, X. Xie, On a Leslie–Gower predator–prey model incorporating a prey refuge, Nonlinear Anal. Real World Appl. 10 (2009) 2905–2908. [14] L. Chen, F. Chen, L. Chen, Qualitative analysis of a predator–prey model with Holling type II functional response incorporating a constant prey refuge, Nonlinear Anal. Real World Appl. 11 (2010) 246–252.

S. Gakkhar, K. Gupta / Applied Mathematics and Computation 273 (2016) 54–67 [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]

67

T.K. Kar, Stability analysis of a prey–predator model incorporating a prey refuge, Commun. Nonlinear Sci. Numer. Simul. 10 (2005) 681–691. J. Pastor, Mathematical Ecology of Populations and Ecosystems, John Wiley & Sons, 2011. J. Chattopadhyay, Effect of toxic substances on a two-species competitive system, Ecol. Model. 84 (1996) 287–289. O. Flaaten, Bioeconomics of sustainable harvest of competing species, J. Environ. Econ. Manag. 20 (1991) 163–180. T.K. Kar, K. Chaudhuri, On non-selective harvesting of two competing fish species in the presence of toxicity, Ecol. Model. 161 (2003) 125–137. M. Mamat, W.M. Sanjaya, Z. Salleh, M. Ahmad, Numerical simulation dynamical model of three-species food chain with Lotka–Volterra linear functional response, J. Sustain. Sci. Manag. 6 (2011) 44–50. E. Chauvet, J.E. Paullet, J.P. Previte, Z. Walls, A Lotka–Volterra three-species food chain, Math. Mag. 75 (2002) 243–255. A.B. Peet, P.A. Deutsch, E. Peacock-López, Complex dynamics in a three-level trophic system with intraspecies interaction, J. Theor. Biol. 232 (2005) 491–503. D. Mukherjee, The effect of prey refuges on a three species food chain model, Differ. Equ. Dyn. Syst. 22 (2014) 413–426. M. Sunaryo, Z. Salleh, M. Mamat, Mathematical model of three species food chain with Holling type-III functional response, Int. J. Pure Appl. Math. 89 (2013) 647–657. S. Sarwardi, P.K. Mandal, S. Ray, Analysis of a competitive prey–predator system with a prey refuge, Biosystems 110 (2012) 133–148. R.K. Upadhyay, S.N. Raw, Complex dynamics of a three species food-chain model with Holling type IV functional response, Nonlinear Anal. Model. Control 16 (2011) 353–374. M.A. Aziz-Alaoui, Study of a Leslie–Gower-type tritrophic population, Chaos Solitons Fractals 14 (2002) 1275–1293. A. Priyadarshi, S. Gakkhar, Dynamics of Leslie–Gower type generalist predator in a tri-trophic food web system, Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 3202–3218. S. Gakkhar, R.K. Naji, Order and chaos in a food web consisting of a predator and two independent preys, Commun. Nonlinear Sci. Numer. Simul. 10 (2005) 105–120. J. Alebraheem, Y. Abu-Hasan, Persistence of predators in a two predators-one prey model with non-periodic solution, Appl. Math. Sci. 6 (2012) 943–956. S. Gakkhar, B. Singh, R.K. Naji, Dynamical behavior of two predators competing over a single prey, Biosystems 90 (2007) 808–817. B. Goh, Stability in models of mutualism, Am. Nat. 113 (1979) 261–275. M. Dhakne, A. Munde, Stability analysis of mutualistic interactions among three species with limited resources for first species and unlimited resources for second and third species, Differ. Equ. Dyn. Syst. 20 (2012) 405–414. N.P. Kumar, N.C. Pattabhiramacharyulu, A three species ecosystem consisting of a prey, predator and a host commensal to the prey, Int. J. Open Probl. Comput. Math 3 (2010) 92–113. S. Gakkhar, A. Priyadarshi, S. Banerjee, Role of protection in a tri-trophic food chain dynamics, J. Biol. Syst. 20 (2012) 155–175. M. Nagumo, Über die lage der integralkurven gewöhnlicher differentialgleichungen, Nippon Sugaku-Buturigakkwai Kizi Dai 3 Ki 24 (1942) 551–559. L. Perko, Differential Equations and Dynamical Systems, 7, Springer Science & Business Media, 2013. H. Freedman, P. Waltman, Persistence in models of three interacting predator–prey populations, Math. Biosci. 68 (1984) 213–231. B. Ermentrout, Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students, 14, SIAM, 2002.