Study of a tri-trophic prey-dependent food chain model of interacting populations

Mathematical Biosciences 246 (2013) 55–71

Contents lists available at ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

Study of a tri-trophic prey-dependent food chain model of interacting populations Mainul Haque a,⇑, Nijamuddin Ali b, Santabrata Chakravarty c a

School of Clinical Sciences, Queen’s Medical Centre Campus, Nottingham University Hospital, Derby Road, Nottingham, UK Department of B.Ed., Katwa College, WB, India c Department of Mathematics, Visva-Bharati, Santiniketan, WB, India b

a r t i c l e

i n f o

Article history: Received 26 March 2013 Received in revised form 8 July 2013 Accepted 17 July 2013 Available online 14 August 2013 Keywords: Intra-specific competition Saddle-node Transcritical Hopf–Andronov Lyapunov functions

a b s t r a c t The current paper accounts for the influence of intra-specific competition among predators in a prey dependent tri-trophic food chain model of interacting populations. We offer a detailed mathematical analysis of the proposed food chain model to illustrate some of the significant results that has arisen from the interplay of deterministic ecological phenomena and processes. Biologically feasible equilibria of the system are observed and the behaviours of the system around each of them are described. In particular, persistence, stability (local and global) and bifurcation (saddle-node, transcritical, Hopf–Andronov) analysis of this model are obtained. Relevant results from previous well known food chain models are compared with the current findings. Global stability analysis is also carried out by constructing appropriate Lyapunov functions. Numerical simulations show that the present system is capable enough to produce chaotic dynamics when the rate of self-interaction is very low. On the other hand such chaotic behaviour disappears for a certain value of the rate of self interaction. In addition, numerical simulations with experimented parameters values confirm the analytical results and shows that intra-specific competitions bears a potential role in controlling the chaotic dynamics of the system; and thus the role of self interactions in food chain model is illustrated first time. Finally, a discussion of the ecological applications of the analytical and numerical findings concludes the paper. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Proposing several mathematical models (including predator– prey, food chain and other ecological models) have become a popular field of research in theoretical ecology. [1] showed that Rosenzweig–MacArthur food chain model admits a sequence of pairs of Bogdanov bifurcation. [2] shows that a three species food chain model exhibits stable steady state, periodic and chaotic solutions. Hopf bifurcation around positive interior equilibrium incorporating delay environment in food chain has been studied in [3]. Kooi and his colleagues illustrated the long-term behaviour of some popular food chain models through numerical bifurcation analysis [4,5] and identified some crucial parameters that control the dynamical behaviour of the model systems. [4] studied for the dynamics of food chains with a nutrient at the base and detected three important regions; a region where the predator goes to extinction, a region for stable equilibrium and a region where a stable limit cycle exists; the global picture for bifurcation diagram by numerical treatments. [6] studied a class of bioenergetic ecological ⇑ Corresponding author. Tel.: + 44 7768971898. E-mail address: [email protected] (M. Haque). 0025-5564/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2013.07.021

models to investigate the dynamics of food chains models with a nutrient at the base where the intake rate of the nutrient and the flux rate for all trophic levels are assumed to be constant. They illustrated that there are three important regions; a region where the predator goes to extinction, a region where there is a stable equilibrium and a region where a stable limit cycle exists. Bifurcation diagrams for tri-trophic food chains are shown which facilitates the study of population model and their dynamical behaviour of a food chain model. [7] studied the asymptotic behavior of a tri-trophic food chain model by finding both local and global bifurcations of equilibria and limit cycles of the proposed model. They have shown the transversal homoclinic orbits to a limit cycle. The appearance of both the homoclinic and heteroclinic orbits from an equilibrium to a limit cycle are computed which has important consequences on the domains of positive equilibrium points. Among all the papers which deal with the food chain models and cited above, the closest one to the current study was first introduced by [8] which describes a prey-dependent food chain model through the following set of differential equations

  dX X M 1 XY  ¼ rX 1  ; dT K A1 þ X

ð1:1aÞ

56

dY E1 M 1 XY M2 YZ ¼  D1 Y  ; dT A1 þ X A2 þ Y dZ E2 M 2 YZ ¼  D2 Z; dT A2 þ Y Xð0Þ > 0; Yð0Þ > 0; Zð0Þ > 0;

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

ð1:1bÞ ð1:1cÞ

where X; Y; Z stand for the population density of prey, predator and top predator respectively. The parameters Mi ; Ai ; Ei ; Di ; Hi (i ¼ 1; 2) represent the maximal predator growth rates, half saturation constants, predator efficiencies, predator’s death rates and predator’s self interaction rates respectively. This model is some times know as standard Rosenzweig–MacArthur model of tri-trophic food chains. This model is studied by several authors by nondimensionalising it according to x ¼ KX ; y ¼ KEY 1 ; z ¼ KEZ1 E2 ; t ¼ rT which lead the model to

dx a1 xy ¼ xð1  xÞ  ; dt 1 þ b1 x dy a1 xy a2 yz ¼  d1 y  ; dt 1 þ b1 x 1 þ b2 y dz a2 yz ¼  d2 z; dt 1 þ b2 y xð0Þ > 0; yð0Þ > 0; zð0Þ > 0;

ð1:2aÞ ð1:2bÞ ð1:2cÞ ð1:2dÞ

where the dimensionless parameters are given (2.3). Next we summarize a brief historical account of the food chain model (1.2). The complex behaviour of the solutions of three-level food chain model (1.2) has been studied in [8,9]. The model (1.2) was first formulated by [8]. Stability analysis were performed and Hopf bifurcation was determined on the same paper. [10] investigated the dynamical behavior of a three species food chain model where the consumer species has a Holling type II functional response. The authors examine how the densities of the top-level population change with the change of carrying capacity of the bottom level. Using a simple geometric approach, [11] detects the singular homoclinic orbits in the model (1.2), as well as parameter combinations for which these orbits exist. [12] found that complex behaviour of this food chain model depends on the functional form chosen to model the interaction between the two highest species in the food chain. The authors concluded that a slight changes in ecological structure may lead significant changes in the bifurcations analysis; although, most ecological data are inadequate to illustrate such phenomena. The model (1.2) revisited by [13] where the authors studied the dynamics of a three species food chain using bifurcation theory to demonstrate the existence of chaotic dynamics in the neighborhood of the equilibrium where the top species in the food chain is absent. The particular goal of this study was to demonstrate the presence of chaos in a class of ecological models, rather than just in a specific model. [14] studied a three-species food chain model of Holling type II functional responses and allometric relationships is analyzed mathematically. A local bifurcation analysis leads to several regions in parameter space that are characterized by rich dynamical behaviours including multiple domains of attraction, quasi-periodicity, chaos, homoclinic events, and transient chaos. [15] focused on two-parameter bifurcation diagram of the model (1.2) which contains Hopf, fold, and transcritical bifurcation curves. The appearance of chaos in this model is proved to be related to a Hopf bifurcation. [16] have shown the global stability of the model around its positive interior equilibrium using a general Lyapunov functions. 2. Mathematical model formulation [17] modified the famous predator–prey model of Lotka-Volterra by incorporating an intra-specific competition in the predator species. Bazykin and his colleagues [18] studied a predator–prey

model with Holling type II predation response function with an intra-specific competition in predator. [19] studied a ratio-dependent predator–prey models of interacting populations which contains an intra-specific competition term in the predator population. This paper showed that predator’s self interaction term has stabilization effect. The author also illustrated how an intra-specific competition among predator could be beneficial for a predator–prey system. [20] studied the original Beddington– DeAngelis predator–prey model which accounts the intra-specific competition among predators. This study also shows that competition among the predator population is beneficial for a number of predator–prey models by keeping them stable around its positive interior steady states when both populations co-exist. The authors showed that predator intra-species competition benefits the predator–prey system under both deterministic and stochastic environments. Observing the importance of self-interaction which has a great role in predator–prey system, we were influenced to modify the famous Rosenzweig–MacArthur model of tri-trophic food chains model Eqs. (1.1) by incorporating intra-specific competitions among predator. Therefore with the above modification the model (1.1) takes the form

  dX X M1 XY  ¼ rX 1  ; dT K A1 þ X dY E1 M1 XY M2 YZ ¼  D1 Y   H1 Y 2 ; dT A1 þ X A2 þ Y dZ E2 M2 YZ ¼  D2 Z  H2 Z 2 ; dT A2 þ Y Xð0Þ > 0; Yð0Þ > 0; Zð0Þ > 0;

ð2:1aÞ ð2:1bÞ ð2:1cÞ

where all the variables and parameters are defined earlier except Hi ’s which stand for the predator’s intra-specific competitions. The system of the Eq. (2.1) has been nondimensionalised according to x ¼ KX ; y ¼ KEY 1 ; z ¼ KEZ1 E2 ; t ¼ rT which lead them to

dx a1 xy ¼ xð1  xÞ  ; dt 1 þ b1 x dy a1 xy a2 yz ¼  d1 y   h1 y2 ; dt 1 þ b1 x 1 þ b2 y dz a2 yz ¼  d2 z  h2 z2 ; dt 1 þ b2 y xð0Þ > 0; yð0Þ > 0; zð0Þ > 0;

ð2:2aÞ ð2:2bÞ ð2:2cÞ ð2:2dÞ

where the dimensionless parameters are given by

M 1 KE1 K D1 H1 KE1 ; b1 ¼ ; d 1 ¼ ; h1 ¼ ; A1 rA1 r r M 2 KE1 E2 KE1 D2 H2 KE1 E2 ; b2 ¼ ; d2 ¼ ; h2 ¼ : a2 ¼ rA2 A2 r r

a1 ¼

ð2:3aÞ ð2:3bÞ

In the current paper we focus the local and global stability of equilibria, persistence of three species, bifurcations and chaos of the system (2.2). Note that the model (1.2) and (2.2) have same number of variables and parameters except hi ’s which represent the intra-specific competition among the predator population in the model (2.2). The rest of the paper is organized as follows. Some preliminary results such as boundedness and dissipativeness is presented in Section 3. Section 4 presents the behaviours of the system (2.2) around all positive equilibrium points including the local stability, global stability and bifurcations. Section 5 illustrates the analytical findings through numerical simulations. The concluding Section 6 contains a final discussion and interpretation of our analytical and numerical results.

57

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

hyperbolic equilibrium point E0 ð0; 0; 0Þ is an unstable saddle point with the unstable manifold in the direction orthogonal to the y–z co-ordinate plane.

3. Preliminary results 3.1. Equilibria

4.2. Behaviour of the system around E1 ð1; 0; 0Þ

The system (2.2) has the following positive equilibria: (i) (ii) (iii) (iv)

Trivial equilibrium: E0 ð0; 0; 0Þ, Axial equilibrium: E1 ð1; 0; 0Þ, Boundary equilibrium: E2 ðx2 ; y2 ; 0Þ, Positive interior equilibrium: E3 ðx3 ; y3 ; z3 Þ, where 1 x2 Þ y2 ¼ ð1x2 Þð1þb and x2 satisfies the following cubic equation a1

3 X Ai xi ¼ 0;

ð3:1Þ

Lemma 1. ðiÞ The equilibrium point E1 ð1; 0; 0Þ of the system is locally asymptotically stable if and only if a1 < d1 ð1 þ b1 Þ. ðiiÞ The system experiences transcritical bifurcation around the axial equilibrium ½tc  ½tc  point E1 at a1 ¼ a1 1 , where a1 1 ¼ d1 ð1 þ b1 Þ. ðiiiÞ E1 ð1; 0; 0Þ is 1 globally asymptotically stable if d1 > 2a ; d2 > ab22 . b1 Proof of Lemma 1: Proof is given in appendix A.6.1.

i¼0

where the coefficients Ai ’s are given

A3 ¼

2 h1 b1 ;

A2 ¼ h1 b1 ð2  b1 Þ;

4.3. Behaviour of the system around E2 ðx2 ; y2 ; 0Þ

A1 ¼ h1  2h1 b1  a1 b1 d1 þ

a21 ;

A0 ¼ ða1 d1 þ h1 Þ: This cubic equation has in general three roots. We find sufficient conditions for the existence of at least one real positive root in Appendix A.2. For positive interior equilibrium E3 ðx3 ; y3 ; z3 Þ, we 1 x3 Þ 2 d2 Þy3 d2 have y3 ¼ ð1x3 Þð1þb ; z3 ¼ ða2hb and x3 satisfies the followa1 2 ð1þb2 y3 Þ ing equation 6 X

Lemma 2. (i) The equilibrium point E2 is locally asymptotically stable d2 1 1Þ if the conditions x2 < ðb2b and y2 < a2 b are satisfied. (ii) The local 1 2 d2

stability of the system ensures its global stability around the boundary equilibrium point E2 provided a1 b1 y2 6 ð1 þ b1 x2 Þ. (iii) The system ½sn2 

experiences a saddle-node bifurcation around E2 at a1 ¼ a1 ½sn  a1 2

Bi xi ¼ 0;

ð3:2Þ

i¼0

is a root of the equation h1 x2 y2 

þ

a21 x2 y2 ð1þb1 x2 Þ

3

, where

¼ 0 and the ½tc 

system experiences a transcritical bifurcation around E2 at a1 ¼ a1 2 , ½tc2 

where are the coefficients Bi ’s are given in Appendix A.1. The Eq. (3.2) has in general six roots. The sufficient conditions for the existence of a positive real root of the Eq. (3.2) are given in Appendix A.3.

a1 b1 h1 x2 y22 2 ð1þb1 x2 Þ

where a1

2 y2 Þ ¼ d2 ð1þb . (iv) The system enters into Hopf-bifurcation y 2

around E2 a1 ¼

at

½hb2 

a1 ¼ a1 2

ðx2 þh1 y2 Þð1þb1 x2 Þ b1 x2 y2

,

where

½hb2 

a1

satisfies the equality

.

Proof of Lemma 2: Proof is given in appendix A.6.2. 3.2. Boundedness 4.4. Behaviour of the system around E3 ðx3 ; y3 ; z3 Þ

Proposition 1. The solution sðtÞ ¼ ðxðtÞ; yðtÞ; zðtÞÞ of system (2.2) which initiate in R3þ (that is, s0 2 R3þ ) are uniformly bounded.

Lemma 3. ðiÞ The equilibrium point is locally asymptotically stable if 2

2

a1 b1 < ð1þby1 x3 Þ and a2 b2 < ð1þbz23 y3 Þ . ðiiÞ The local stability of the system 3

Proof. The proof is provided in the Appendix A.4.

h

ensures its global stability around the interior equilibrium point E3 , provided a1 y3 6 ð1 þ b1 x3 Þ and a2 z3 6 h1 ð1 þ b2 y3 Þ. ðiiiÞ The system ½sn3 

experiences a saddle-node bifurcation around E3 at a1 ¼ a1

3.3. Dissipativeness

, where

½sn  a1 3

satisfies the equation m11 m22 m33  m12 m21 m33  m11 m23 m32 ¼ 0, where mij ’s are given by Proposition 2. If ai > bi di , i ¼ 1; 2, then the system (2.2) is dissipative.

m21 Proof. The proof is provided in the Appendix A.5.

h

m23 4. Stability and bifurcation analysis In this section, we deal with the local stability, global stability and the bifurcation analysis of the system (2.2). The proofs of the stability and bifurcations results are described in appendix A.6. Here we report the results only.

a1 b1 x3 y3 a1 x3 ; m12 ¼  ; m13 ¼ 0; 2 1 þ b1 x3 ð1 þ b1 x3 Þ a1 y3 a2 b2 y3 z3 ¼ ; m33 ¼ h2 z3 ; m22 ¼  h1 y3 ; 2 2 ð1 þ b1 x3 Þ ð1 þ b2 y3 Þ a2 y3 a2 z3 ¼ ; m32 ¼ ; m31 ¼ 0; m32 ¼ 0: 2 1 þ b2 y3 ð1 þ b2 y3 Þ

m11 ¼ x3 þ

½hb 

(iv) The system exhibits a Hopf-bifurcation at E3 for a1 ¼ a1 3 , where ½hb  ½hb  ½hb  ½hb  a1 3 satisfies the equation ðm11  a1 3 Þðm22  a1 3 Þðm33  a1 3 Þ ½hb  ½hb  m12 m21 ðm33  a1 3 Þ  m23 m32 ðm11  a1 3 Þ ¼ 0. Proof of Lemma 3: Proof is given in A.6.3. 4.5. Persistence

4.1. Behaviour of the system around the origin E0 ð0; 0; 0Þ A straightforward calculation shows that the eigenvalues of the system (2.2) around E0 ð0; 0; 0Þ are 1; d1 and d2 . Therefore the

We have already shown that system (2.2) is uniformly bounded. In order to prove the persistence of the system, we first show that all the boundary equilibria are repellers.

58

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

Proposition 3. If thee following conditions hold ðiÞ d2 1 1Þ a1 > d1 ð1 þ b1 Þ; ðiiÞ x2 < ðb2b and y > then the system is 2 a b d 1 2 2 2 persistent.

Proof. From Lemma 1, we observe that if a1 > d1 ð1 þ b1 Þ then E1 is unstable. From Lemma 2, we notice that E2 is unstable for d2 1 1Þ x2 < ðb2b and y2 > a2 b . Hence combination of all these condi1 2 d2 tions yield the result. h

5. Numerical simulations One of the aims of this section to support the analytical findings with the help of experimented parameter values from published literature. For this purpose we report, in the Table 1, the parameter values of the eight parameters of the model (2.2) used by several authors to simulate the behavior of the nondimentionalised food chains model (1.2). We tried to take the parameter values within the corresponding ranges shown in The Table 1. Here we report the parameter values that have been used in the current study, that is, in the modified model (2.2); and they are a1 ¼ 0:5; 3:39; 5; 6:109; 6:39; a2 ¼ 0:01; 2:3; 2:8; 3; b1 ¼ 2; 2:05; 2:2; 2:3; 3; b2 ¼ 0:5; 2; d1 ¼ 0:5; 0:7; 1; d2 ¼ 0:4; 0:5; 1:055. The parameters values used in the current study are within the corresponding ranges given in the Table 1 except b1 ¼ 2 and d2 ¼ 0:4; 0:5; 1:055. [19–21] took h1 ¼ 0:01 per year. We used h1 ¼ 0:03; 0:05; 0:06; 0:075; 0:09; 0:4; 0:56; 0:75 and h2 ¼ 0:00024; 0:000255; 0:00026; 0:0003; 0:000315; 0:0004; 0:05, that are closely related to the values hi ’s that have been calculate from the corresponding nondimentionalised parameters hi ’s of the Eq. (2.3) using the

Table 1 The table reports the parameter values of the six parameters of the model (2.2) used by several authors to simulate the behavior of the model (1.2). The values indicated by ⁄ are calculated as the functional form of predator response function are slightly different from that of (1.2). The values indicated by   are the center values of the range used in the numerical simulations. a1 y

1:75 5:0 1:81 0:22 2:4 8:0

a2

b1

0:3

2:615

0:1 0:181 0:011 0:022 2:88

4:1 4:5 5:0 2:98 6:66

b2 y

1:749 2:0 0:46 0:25 0:44 2:4

d1

d2

Referred Papers

0:1

0:04

Abrams and Roth [10]

0:4 0:16 0:02 0:4 0:87

0:01 0:08 0:02 0:01 0:25

Hastings and Powell [8] Hogeweg and Hesper [28] Rai and Sreenivasan [29] Mccann and Yodzis [14] Scheffer [30]

parameter values in Table 1. One of the integral purposes of this section is to verify the analytical conditions with the help of reasonably close experimental data reported in the published literature. Besides verification of our analytical findings, these numerical simulations shows important features of the system from practical point of view. Next we present numerical simulations of the proposed model (2.2). Fig. 1 exhibits the global stability of the model (2.2) around E1 with the parameter values (given in the corresponding figures) for which the conditions for local stability and global stability of Lemma 1 are satisfied (see Fig. 2). Fig. 3 demonstrates the transcritical bifurcation of the system (2.2) around the predator-free axial equilibrium E1 at ½tc  a1 ¼ a1 1 ¼ 1:5 which also matches with the figure drawn for other dynamical systems at a transcritical point (cf. [22]). Fig. 7(a) illustrates the global stability of the system (2.2) around E2 which indicates that without a top-predator, predator species coexist with prey species in some situations in a stable form. The parameter values satisfy the condition mentioned in the Lemma 2. In some ecological situation the top predator represents the pest and thus the extinction of top-predator is desirable; and Lemma 2 gives the conditions for which we observe this situation. Lemma 2 describes that the system (2.2) experiences a trans½tc  critical bifurcation around E2 for a1 ¼ a1 2 ; the Fig. 3 demonstrates ½tc2  this fact at a1 ¼ a1 ¼ 1:57. System (2.2) attains a saddle-node ½sn  bifurcation around E2 at a1 ¼ a1 2 ¼ 6:9603; see Fig. 4 which is similar to the figure described for other predator–prey systems at saddle-node point (cf. [22]). The solution of the system (2.2) settle down to a stable steady state when the parameters of the system are chosen as h1 ¼ 0:56; h2 ¼ 0 and b1 2 ð2:3; 6:2 which has been well illustrated in Fig. 5(a). A little further increase of h2 from 0:000315 to 0:0004 without varying other parameters yields the system to posses a steady stable state as noted in Fig. 5(b). System (2.2) passes through a saddle-node bifurcation around E3 at ½sn3 a1 ¼ a1 ¼ 7:094; see Fig. 6 which is similar to the figure described for other predator–prey systems at the saddle-node point (cf. [22]). The Fig. 7(a) confirms the extinction of the top-predator. These extinction results suggest that biological control is possible in such situations. However, in practice, we often want to reduce the pest y to an acceptable level in a finite time. In order to accomplish that, we need to give an estimate as how much z can be introduced to control y. This confirms the need of early planning of a biological control. The occurrence of Hopf-bifurcation in absence of top-predator is shown in Fig. 8. Another biologically interesting question remains open for the model (2.2) is under what conditions, all three species coexist. This

1.6 Prey Predator Top−predator

1.4 1.2

1.5

1

1

0.8

Species

Top−predator

2

0.5

0.6 0.4

0 2

0.2 1.5

2 1.5

1

Predator

0

1

0.5

0.5 0

0

−0.2 Prey

0

200

400

600

800

1000

Time

Fig. 1. Global stability of the system around E1 for parameter values a1 ¼ 0:5; a2 ¼ 0:1; b1 ¼ 3:0; b2 ¼ 2:0; d1 ¼ 0:7; d2 ¼ 0:5; h1 ¼ 0:05; h2 ¼ 0:05. Note that for this set of 1 parameter values we have a1 ¼ 0:7 < d1 ð1 þ b1 Þ ¼ 2:8; d1 ¼ 0:7 > 2a ¼ 0:33 and d2 ¼ 0:5 > ba22 ¼ 0:05. b1

59

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71 0.15 1.08 0.1 Transcritical tcBifurcation at a1=a1

1.06

Predator (y)

Prey (x)

0.05

Transcritical Bifurcation tc at a1=a1

1.04 1.02 1

0

−0.05

0.98

−0.1

0.96 1

1.2

1.4

1.6

1.8

2

1

1.2

1.4

1.6

a

1.8

2

a

1

1

½tc 

Fig. 2. Transcritical bifurcation around E1 at a1 ¼ a1 1 ¼ d1 ð1 þ b1 Þ ¼ 1:5. The red line represents the loci of stable equilibria, the black line the unstable ones. The other parameter values are a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. We do not report of top predator here but in the final figure where all bifurcation are combined in Figs. 10–12.

1

0.2

0.98

0.19

0.96

0.18

0.94

0.17

0.92

0.16

0.04 0.03

0.9 0.88

0.14 0.13

at a =atc2 1

2

0.84

0.15

Transcritical Bifurcation

Transcritical Bifurcation at a =atc2 1 1 around E

0.86

Top predator (z)

Predator (y)

Prey (x)

0.02

1.2

1.4 a

Transcitical Bifurcation

−0.02

1

at a =atc2

1.6

1.8

0.1

2

1

around E2

0.11 1

0 −0.01

0.12

0.82 0.8 0.8

0.01

1

−0.03

1.2

1.4

1.6

1.8

−0.04 1.5

2

1.52

1.54

a

1

1

around E2 1.56

1.58

1.6

a

1

1

½tc 

Fig. 3. Transcritical bifurcation around E2 at a1 ¼ a1 2 ¼ 1:57. The red line represents the loci of stable equilibria, the black line the unstable ones. The other parameter values are a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05.

0.4

0.16

0.3

0.159

0.35

0.25

0.158 0.2

0.3

Predator (y)

Prey (x)

0.25 0.2 0.15

Top predator (z)

0.157 0.156 0.155 0.154

0.15 0.1 0.05

0.153 0

0.1

0.152

0.05

−0.05

0.151 5

5.5

6

6.5

7

7.5

8

0.15

½sn2 

Fig. 4. Saddle-node bifurcation around E2 at a1 ¼ a1

5

5.5

6

6.5

7

7.5

8

−0.1

6

6.5

7 a1

7.5

8

¼ 6:9603. The red line represents the loci of stable equilibria, the black line the unstable ones.

scenario is shown in the global stability of the system (2.2) around E3 (cf. Fig. 7(b)). By changing the parameter values we have shown that the system (2.2) enters into a Hopf-bifurcation around E3 . This happens, for instance, if we choose a1 ¼ 6:109; a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05 and h2 ¼ 0:05. This situation has been illustrated in Figs. 9. Now we are in a position to combine all the bifurcations points in one figure to observe carefully how the system (2.2) experiences these bifurcations.Figs. 10–12 illustrates a detailed bifurcation behaviour of the system (2.2) while the parameter a1 varies from 1 and 8. Let us focus on the classifications of the Hopf-bifurcation of the system (2.2) around E2 and E3 which are observed, respectively, at

½hb 

½hb 

a1 ¼ a1 2 ¼ 3:353 and a1 ¼ a1 3 ¼ 6:839 . Using the method outlined in [23] one can easily determine the first Lyapunov coefficients to determine the stability of the limit cycles. Figs. 13–15 illustrate whether the Hopf-bifurcations around E2 and E3 are subcritical or supercritical. The branch of periodic orbits that comes out of the Hopf point E2 is unstable (blue open circles) and then it turns around at a limit point of oscillations at a1 ¼ 6:73 so that a stable (blue solid circles) and unstable limit cycle coalesce here. Similarly, the branch of periodic orbits that comes out of the Hopf point E3 is stable (green solid circles). We illustrate the greater view of this stable limit cycles in Fig. 9. The periodic solution loses stability and that is, there is a period-doubling instability. Now we

60

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

Fig. 5. Steady state solution of the system (2.2) a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01 and b1 ¼ 3:0.

0.45

for

(a)

h1 ¼ 0:56; h2 ¼ 0;

and

for

(b)

h1 ¼ 0; h2 ¼ 0:0004.

other

parameter

values

are

0.35

0.161

0.4

The

0.16

0.3

0.35 0.159 0.25

0.25 0.2

sn at a =asn3 1

1

0.158 sn at

Top predator (z)

Predator (y)

Prey (x)

0.3 sn a1=a1 3

0.157 0.156

0.15 0.1

0.2 sn 3 1

sn at a =a 1

0.15

0.155

0.1

0.05 0.154 6.5

7 a

7.5

0.05

6.95

7

7.05 a

1

7.1

7.15

7

7.02

7.04

7.06

7.08

7.1

a

1

1

½sn 

Fig. 6. Saddle-node bifurcation around E3 at a1 ¼ a1 3 ¼ 7:094. The red line represents the loci of stable equilibria, the black line the unstable ones. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

0.8 Prey Predator Top−predator

0.7

0.4

0.6

0.3 Top−predator

Species

0.5 0.4 0.3

0.2 0.1

0.2 0 0.4

0.1

0.3 0 −0.1

1 0.2

0.8 0.1

0

200

400

600

800

1000

Time

E2

Predator

0.6 0

0.4

Prey

E3

Fig. 7. (a) Global stability of the system around E2 for parameter values a1 ¼ 5:0; a2 ¼ 0:1; b1 ¼ 3:0; b2 ¼ 2:0; d1 ¼ 1:0; d2 ¼ 1:055. (b) Global stability of the system around E3 for parameter values a1 ¼ 3:39; a2 ¼ 3:0; b1 ¼ 2:9; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4. The other parameter values are h1 ¼ 0:05; h2 ¼ 0:05.

observe the exchange of different states of stability, namely the limit cycle, period doubling and finally the chaotic behaviour of the system (2.2). Subsequently, it has been observed how this chaotic dynamics gets transferred into a stable dynamics for different values of h1 and h2 . For h1 ¼ 0; h2 ¼ 0 and b1 2 ð0; 2:1Þ, the solutions of the system (2.2) remain stable steady state solution which

is depicted in Fig. 16(a). Chaotic dynamics of the system (2.2) occurs corresponding to the parameter value b1 2 ð2:3; 6:2 and h1 ¼ 0; h2 ¼ 0 as pointed out in Fig. 16(b). For h1 ¼ 0; h2 ¼ 0 and b1 2 ð2:1; 2:3Þ the trajectories enter into limit cycle oscillations from steady state position. For instance b1 ¼ 2:2, we observe the Fig. 17(a). Further enter into a situation of period doubling for

61

0.8

0.8

0.7

0.7

0.7

0.6

0.6

0.6

0.5

0.5

0.5

0.4

Top−predator

0.8

Predator

Prey

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

0.4

0.4

0.3

0.3

0.2

0.2

0.2

0.1

0.1

0.1

0

0

200

400

600

800

1000

0

0

200

400

Time

600

800

1000

0.3

0

0

200

400

600

Time

800

1000

Time

Fig. 8. Hopf-bifurcation around E2 for parameter values a1 ¼ 6:39; a2 ¼ 2:3; b1 ¼ 3:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:5; h1 ¼ 0:05; h2 ¼ 0:05.

0.7

0.6

0.6

0.5

0.5

0.4

0.4

Prey

Predator

0.7

0.3

0.3

0.2

0.2

0.1

0.1

0

0

200

400

600

800

0

1000

0

200

400

600

800

1000

Time

Time

0.18 0.16 0.2

0.14

0.15 Top predator

Top−predator

0.12 0.1 0.08

0.1 0.05

0.06 0 0.8

0.04

0.6 0.02 0

0.8 0.6

0.4

0.4

0.2 0

200

400

600

800

1000

Time

Predator

0.2 0

0

Prey

Fig. 9. Hopf-bifurcation around E3 where solution trajectories are functions of time (a,b,c) a1 ¼ 6:109; a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05 and h2 ¼ 0:05. (d) Phase portart of Hopf-bifurcation.

b1 ¼ 2:3 (cf. Fig. 17(b)). Moreover, for h1 ¼ 0:075; h2 ¼ 0 and b1 2 ð2:3; 6:2 the trajectories of the dynamical system (2.2) enter into period doubling situation from chaotic dynamics as evident from Fig. 18(a). The Fig. 18(b) describes the behaviour of the trajectories of the system (2.2) as how they enter into the different phases of stability depending upon the non-zero parametric variations of both h1 and h2 for a constant parameter values of a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01 and b1 ¼ 3:0. The limit cycle behaviour of the system (2.2) for h1 ¼ 0:4; h2 ¼ 0 and b1 2 ð2:3; 6:2 is illustrated in Fig. 19(a), where the system enters

into the state of limit cycle oscillations from period doubling phase. Fig. 19(b) depicts a similar behaviour for h1 ¼ 0:06; h2 ¼ 0:00026, other parameter values are same as that of 19(a). Furthermore, the results of Fig. 20(a) exhibit the transfer of trajectories of the system (2.2) into the period doubling phase from chaotic dynamics when the parameters h1 ; h2 ; b1 are selected to have their respective values as 0; 0:000255 and 3:0 while they are found to enter into the state of limit cycle oscillations from period doubling phase for a slight change of the parameter h2 value from 0:000255 to 0:000315 by keeping other parameters unperturbed as observed

62

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71 0.5

1

0.45

tc at E

hb at E2

0.4

1

sn at E

3

0.35

tc at E

2

Predator (y)

0.8 hb at E

Prey (x)

3

0.6

0.3 0.25

sn at E3

0.1

sn at E

2

0.2

tc at E1

sn at E2

0.2

0.05

2

hb at E3

0.15

0.4 hb at E

tc at E2

0 1

2

3

4

5

6

7

8

a1

0

1

2

3

4

5

6

7

a

1

0.17 hb at E3 0.165

0.35

0.16

0.3

hb at E

sn at E

3

Prey (x)

Predator (y)

0.4

3

sn at E

3

0.155

0.25

0.15

0.2

0.145

sn at E

2

0.15

0.14

sn at E

2

6.75

0.1 0.05 6.5

6.8

6.85

6.9

6.95

7

7.05

7.1

a1

6.6

6.7

6.8

6.9

7 a

7.1

7.2

7.3

7.4

7.5

1

Fig. 10. All bifurcation situation of the system (2.2) for prey (x) around all possible positive steady states E0 ; E1 ; E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. Here tc, sn and hb stand for transcritical bifurcation, saddle-node bifurcation, and Hopf-bifurcation.

in Fig. 20(b). It is important to note that the solutions of the system (2.2) enter the period doubling phase from chaotic situation for particular values of h1 ¼ 0:03 and h2 ¼ 0:00024 which subsequently alter to enter into the state of limit cycle oscillations from the phase of period doubling based upon the slightly increased values of h1 ¼ 0:06 and h2 ¼ 0:00026 and eventually they settle down to a steady state stable condition for a final choice of the significant parameters h1 ¼ 0:09 and h2 ¼ 0:0003. The other parameter b1 has been kept unperturbed with its constant value 3:0 over the entire results of these figures.

6. Discussion Classical prey dependent food chain models (2.2) were studied by many authors since the pioneering work of [24,25] to get a better insight of the behaviour of the solution trajectories. In the present investigation, an attempt has been made to discuss the phenomena of the trajectories of a more suitable natural situation by using a modified classical prey-dependent food chain model where we incorporate the predator and top predator inter-specific competitions to the existing models ([24,25]). To summarise the analytical findings and the evidence from numerical simulations, the proposed model exhibited four biologically feasible equilibria. The first equilibrium point is point where all populations are extinct which is always unstable. The second one is where the prey is at its highest population while the predator and the top predator

Fig. 11. All bifurcation situation of the system (2.2) for predator (y) around all possible positive steady states E0 ; E1 ; E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. Here tc, sn and hb stand for transcritical bifurcation, saddle-node bifurcation, and Hopf-bifurcation.

are extinct. The third one is where the prey and the predator are persist while the top predator is extinct. Finally, the fourth one is the equilibrium where all three species exist. In Table 2, we give a comparison between the results of the system (1.2) and the results of the system (2.2). The results of the system (1.2) is obtained in Appendix A.7. The first two columns represent the results of the system (1.2) in terms of our system parameters. The third and fourth columns represent the summary of the system (2.2) around E2 . The first two rows represent the results about local and global stability of the axial equilibrium E1 . This steady sate is globally asymptotically stable in the parametric region a1 < d1 ð1 þ b1 Þ; 1 ; d2 > ba22 . Third row of Table 2 shows that local stability of d1 > 2a b1 the system around E1 influences by self-interaction among predators. Fourth row represents the global stability results around E2 . The third is the boundary equilibrium point where the top-predator is extinct but the predator and prey persist. When h1 ¼ 0 and h2 ¼ 0, [16] showed (Sec. 3 page-377) that the boundary equilibrium is locally asymptotically stable for asymptotically stable for

KA1 2

M 2 y A2 þy

6 D2 and globally

< x < K which in terms of our sys-

d2 and tem parameter reduces to y2 6 a2 b 2 d2

1b1 2

< x2 < 1. However

in our model (where h1 – 0 and h2 – 0), the boundary equilibrium d2 1 1Þ is globally asymptotically stable for x2 < ðb2b ; y2 > a2 b and 1 2 d2

a1 b1 y2 6 ð1 þ b1 x2 Þ. We conclude that stability of E2 needs more conditions when predators intra-species competition is considered. Again, for the parameter values

b1 1 2b1

2

1 < 1b , i.e. for, b1 < 1, 2

our system is unstable whereas their system is still globally stable, which implies that extinction of predator population is possible for

63

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71 Y

0.4

0.7

0.35 hb at E3

0.3 Top predator (z)

0.6

0.25

sn at E

3

0.2

0.5

0.15

tc at E

1

0.1

sn at E2

0.4

hb at E2

0.3

tc at E

2

0.05 0 −0.05

1

2

3

4

5

6

0.2

7

a1 0.1

0.4 0

0.35

1

hb at E

Top predator (z)

0.3

3

4

5

6

7

8

a1

Fig. 14. The Hopf-bifurcation situation of the predator of the system (2.2) around E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. The branch of periodic orbits that comes out of the Hopf point E2 is unstable (blue open circles) and then it turns around at a limit point of oscillations at a1 ¼ 6:73 so that a stable (blue solid circles) and unstable limit cycle coalesce here. Similarly, the branch of periodic orbits that comes out of the Hopf point E3 is stable (green solid circles).

0.25 0.2

2

3

sn at E3

0.15 sn at E2

0.1 0.05 0 −0.05 6.5

6.6

6.7

6.8

6.9

7 a1

7.1

7.2

7.3

7.4

7.5

Z 0.45 0.4

Fig. 12. All bifurcation situation of the system (2.2) for top predator (z) around all possible positive steady states E0 ; E1 ; E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. Here tc, sn and hb stand for transcritical bifurcation, saddle-node bifurcation, and Hopf-bifurcation.

0.35 0.3 0.25 0.2

X

0.15 0.1

1.4

0.05 1.2 0 1

-0.05 0

1

2

0.8

3

4 a1

5

6

7

8

Fig. 15. The Hopf-bifurcation situation of the top predator of the system (2.2) around E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. The branch of periodic orbits that comes out of the Hopf point E2 is unstable (blue open circles) and then it turns around at a limit point of oscillations at a1 ¼ 6:73 so that a stable (blue solid circles) and unstable limit cycle coalesce here. Similarly, the branch of periodic orbits that comes out of the Hopf point E3 is stable (green solid circles). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

0.6

0.4

0.2

0 0

2

4

6

8

10

a1

Fig. 13. The Hopf-bifurcation situation of the prey of the system (2.2) around E2 and E3 for the parameter values a2 ¼ 2:8; b1 ¼ 2:0; b2 ¼ 0:5; d1 ¼ 0:5; d2 ¼ 0:4; h1 ¼ 0:05; h2 ¼ 0:05. The branch of periodic orbits that comes out of the Hopf point E2 is unstable (blue open circles) and then it turns around at a limit point of oscillations at a1 ¼ 6:73 so that a stable (blue solid circles) and unstable limit cycle coalesce here. Similarly, the branch of periodic orbits that comes out of the Hopf point E3 is stable (green solid circles).

three predator–prey models where inter-specific competitions have been incorporated ([18–20]) The system (1.2) has already been studied in ([16–26]) but unfortunately no one showed the global stability of the interior equilibrium point E3 where all three individuals were present. We showed that E3 of the system (2.2) incorporating the self-interaction term is globally asymptotically stable

their model but is not possible for our model. Therefore, one may remark that self-interaction among the predator population (with h1 – 0 and h2 – 0) might be beneficial for the predator species under certain conditions. Similar phenomena have been observed in

in

the

region

2

2

a1 b1 < ð1þby13x3 Þ ; a2 b2 < ð1þbz23y3 Þ ; a1 y3 6

ð1 þ b1 x3 Þ and a2 z3 6 h1 ð1 þ b2 y3 Þ. In addition, the conditions for persistence of the system (2.2) are also been determined successfully. It has been observed that the axial equilibrium E1 exhibits ½tc 

½tc1 

a transcritical bifurcation at a1 ¼ a1 1 , where a1

¼ d1 ð1 þ b1 Þ;

64

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

Fig. 16. Chaotic behaviour of the system (2.2) for (a) b1 ¼ 2:05 and (b) b1 ¼ 3:0. The other parameter values are a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01; h1 ¼ 0 and h2 ¼ 0.

Fig. 17. (a) Limit cycle behaviour of the system (2.2) for parameter values b1 ¼ 2:2; (b) Period doubling behaviour of the system (2.2) for b1 ¼ 2:3. The other parameter values are a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01; h1 ¼ 0 and h2 ¼ 0.

Fig. 18. Period doubling behaviour of the system (2.2) for (a) h1 ¼ 0:075; h2 ¼ 0; and for (b) h1 ¼ 0:03; h2 ¼ 0:00024. The other parameter values are a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01 and b1 ¼ 3:0.

65

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

Fig. 19. Limit cycle behaviour of the system (2.2) a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01 and b1 ¼ 3:0.

for

(a)

h1 ¼ 0:4; h2 ¼ 0;

and

for

(b)

h1 ¼ 0:06; h2 ¼ 0:00026.

The

other

parameter

values

are

Fig. 20. (a) Period doubling behaviour of the system (2.2) for h2 ¼ 0:000255. (b) Limit cycle behaviour of the system (2.2) for h2 ¼ 0:000315. The other parameter values are a1 ¼ 5:0; a2 ¼ 0:1; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01; b1 ¼ 3:0 and h1 ¼ 0.

Table 2 The table describes the results of the system (1.2) and the system (2.2). Conditions of the system (1.2)

Results of the system (1.2)

Conditions of the system (2.2)

Results of the system (2.2)

a1 < d1 ð1 þ b1 Þ

E1 is LAS E1 is GAS

a1 < d1 ð1 þ b1 Þ

E1 is LAS E1 is GAS

d2 y2 6 a2 b 2 d2

E2 is LAS

d2 1 1Þ x2 < ðb2b ; y2 < a2 b 1 2 d2

E2 is LAS

1b1 2

E2 is GAS

a1 b1 y2 6 ð1 þ b1 x2 Þ

E2 is GAS

E3 is LAS

a1 b1 < ð1þby1 x3 Þ ; a2 b2 < ð1þbz23 y3 Þ

1 d1 > 2a ; d2 > ab22 b1

< x2 < 1 2

2

a1 b1 < ð1þby1 x3 Þ ; ½x3 ð1 þ b1 x3 Þ þ a1 b1 x3 y3 a2 b2 y3 z3 3

a21 x3 y3 ð1

2

1 d1 > 2a ; d2 > ba22 b1

2

2

E3 is LAS

3

2

þ b2 y3 Þ  a2 y3 z3 ð1 þ b1 x3 Þ < 0 a1 y3 6 ð1 þ b1 x3 Þ; a2 z3 6 h1 ð1 þ b2 y3 Þ

the boundary equilibrium E2 exhibits saddle-node and transcritical bifurcation for the suitable value of a1 (see Lemma 2ðiiiÞ), while the interior equilibrium E3 exhibits saddle-node bifurcation for specific values of a1 (see Lemma 3ðiiiÞ). The conditions for Hopf bifurcations to arise around E2 and E3 are obtained. Fifth row says the results about local stability of the interior equilibrium E3 of the system (2.2). Sixth row represents the results about the global stability of the interior equilibrium E3 of the considered system (2.2). By putting h1 ¼ 0 and h2 ¼ 0 in the model (2.2), one can reduce it to the same model studied by ([8–13]). They showed that the model

E3 is GAS

exhibits chaotic dynamics in long-term behaviour when biologically reasonable parameter values (a1 ¼ 5:0; a2 ¼ 0:1; b1 ¼ 3:0; b2 ¼ 2:0; d1 ¼ 0:4; d2 ¼ 0:01) are chosen. Our simulation shows that when h1 – 0 or h2 – 0, the system bears the potential ability to control chaotic dynamics for the same parameter values used by [8]. That is, self interactions among predators are largely responsible to transform the system dynamics from chaos to limit cycle, stable steady state situation. Finally, the present analytical results were confirmed by appropriate numerical simulation with realistic parameter values. Therefore, the observations that have

66

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

been made in the present investigation will be of some help to the experimental ecologists who are engaged to perform their work and develop to understand some ecological phenomena of nature one step ahead of the real situation in near future. Under the conditions mentioned in the Table 2, we observe that the axial equilibrium E1 is stable (both locally and globally). Ecologically this means that the both the predator and top-predator are unable to reproduce fast enough to compensate their death rates di (see the second column of Table 2 which shows that the consumption rate ½tc 

a1 should be smaller than a certain threshold, such as a1 1 ). Whatever the amount of prey consumed (that is, a1 ), this food does not allow them to reproduce fast enough; such predators and top-predators are, therefore, called ‘‘maladapted to the prey’’. We show the global stability of both the systems (1.2) and (2.2) around its axial equilibrium point by suitable Lyapunov function. Table 2 also reflects that the qualitative behaviours of the systems (1.2) and (2.2) are same around E1 , more clearly, all these systems attain transcritical bifurcation, local asymptotically stable and global asymptotically stable around its axial equilibrium point. However, it significantly changes the quantitative dynamics which can also be seen from third column of Table 2. This column also illustrates that the critical value of prey consumption rate a1 of (1.2) and (2.2) for which these systems experience transcritical bifurcation are same. The ecologically it reflects that there is no effect of predator intra-species competition rate hi in the qualitative dynamics of the system as well as in critical values of a1 around their axial steady state E1 . Lemma 2 illustrates that if the consumption rate is med½tc 

½hb 

equilibrium. Lemma 2(iii) gives the value of a1 2 for which the system has oscillatory behaviour around E2 ; these are reasonable behaviours of a three species food-chain models that are supported by experimental evidence, see [27] and references therein. These are quite natural behaviors of a food-chain model and have also been observed in several food chain models studied earlier ([24– 26]). Once the values of a1 is increased then system (2.2) has trans½tc2 

and it becomes stable around E3

½hb3 

< a1 < a1

critical bifurcation at a1 ¼ a1 ½tc2 

½hb3 

< a1 < a1

cation at a1 ¼

½hb  a1 3 .

and a1

½sn 

½sn2 

½sn  a1 2

and keeps stable until ½sn 

it a1 < a1 3 . It has saddle-node bifurcation at a1 ¼ a1 3 and the system (2.2) becomes again stable around E3 once a1 becomes greater ½sn3 

A.1. The coefficients Bi ’s The coefficients Bi ’s are given by 3 2

B6 ¼ b1 b2 h1 h2 ;

ðA:1aÞ

2 2

2 2

B5 ¼ b1 b2 h1 h2 ðb1  1Þ þ b1 b2 a1 h2 ða1  b1 d1 Þ;

ðA:1bÞ

2

B4 ¼ h1 h2 b1 b2 ½b1 b2  3b2 ðb1  1Þ  2b1 ða1 þ b2 Þ þ a1 b1 b2 ½b1 b2 d1  2b2 ðb1  1Þða1  b1 d1 Þ; 2

2

ðA:1cÞ

3

B3 ¼ h1 h2 ½2b1 b2 ðb1  1Þ þ b2 ðb1  1Þ  4b1 b2 ðb1  1Þ 2

2

 ða1 þ b2 Þ þ a1 h2 ½b2 ða1  b1 d1 Þðb1  1Þ  2b1 b2 ða1 2

3

 b1 d1 Þða1 þ b2 Þ þ 2b1 b2 d1 ðb1  1Þ ; 2

2

ðA:1dÞ 2

B2 ¼  h1 h2 ½b2 ðb1  1Þ  2b1 b2 ða1 þ b2 Þ þ 2b2 ðb1  1Þ ða1 þ b2 Þ 2

 b1 ða1 þ b2 Þ  þ a21 a2 b1 ða2  b2 d2 Þ 2

2

þ a1 h2 ½2b2 ðb1  1Þða1 þ b2 Þ þ 2b1 b2 d1 ða1 þ b2 Þ  b2 d1 ðb1  1Þ ; ðA:1eÞ B1 ¼  h1 h2 ðb1  1Þða1 þ b2 Þ½2b2 þ ða1 þ b2 Þ  a21 a2 ða2  b2 d2 Þðb1  1Þ þ a1 h2 ða1 þ b2 Þ½ða1  b1 d1 Þða1 þ b2 Þ  2b2 d1 ðb1  1Þ; ðA:1fÞ 2

2

B0 ¼ h1 h2 ða1 þ b2 Þ  a21 a2 ða2  b2 d2  a1 d2 Þ  d1 ða1 þ b2 Þ : ðA:1gÞ

and has a Hopf bifur-

The system (2.2) becomes stable around E2

once a1 passes through the values a1 ¼

than a1

Appendix A

½tc 

ium (a1 1 < a1 < a1 2 ), then both the predator (but not the toppredator) and prey populations coexist in the form of boundary

while a1

Assistance Programme (SAP-II) sponsored by the University Grants Commission, India.

. Therefore the system keeps stable around E3 for all val½sn 

ues of a1 > a1 3 . In this article, we obtained both analytical and the numerically conditions for which the system experiences all kinds of stable and bifurcating situations. Note that the parameter values for which the system (2.2) experiences all kind of bifurcations lies within parameter ranges reported in existing published literature (see Table 1). The ecological significance of these results is that when actual ecological dynamics in the field are observed and when physical parameters are varied qualitatively with time, such as outbreaks of disease, then these systems can show either kind of stable or bifurcating situations. In other words, modification of the model (1.2) by introducing an extra term for predator and toppredator competition enables them to produce richer dynamics that may arise from several natural ecological calamities for which parameter values vary greatly.

A.2. Condition of at least one positive root of Eq. (3.1) Here we find sufficient conditions for the existence of a positive real root of the Eq. (3.1) and leave other possibilities for further investigation. Assuming that there shall be a complex root c and also its conjugate c so that this pair of complex roots build the quadratic with negative discriminant x2 þ bx þ d ¼ ðx  cÞðx  c Þ ¼ x2  2ReðcÞx þ jcj2 , where b and d are real constants with d ¼ jcj2 > 0. Eq. (3.1) will then have the following factorisation 3 X Ai xi ¼ A3 ðx  aÞðx2 þ bx þ dÞ i¼0

¼ A3 ½x3 þ ðb þ aÞx2 þ ðba þ dÞx þ da;

where a is to be determined. By equating coefficients of like powers of x on the left and the right, we find

bþa¼

The authors Nijamuddin Ali and Santabrata Chakravarty gratefully acknowledge the financial support in part from Special

A2 A1 A0 ; d þ ba ¼ ; da ¼ ; A3 A3 A3

ðA:3aÞ

from which we have

a¼ Acknowledgement

ðA:2Þ

A2 A0 ¼ A3 b A3 jcj2

ðA:4aÞ

The other root of Eq. (A.2) is thus found as x2 ¼ a. By imposing conditions a < 0, we obtain x2 > 0, which ensures the uniqueness and feasibility of the boundary equilibrium.

67

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

A.3. The sufficient conditions for a positive real root of the Eq. (3.2) We find sufficient conditions for the existence of a positive real root of the Eq. (2.3) and leave other possibilities for further investigation. Assuming that there shall be two complex roots d and g along with their conjugates d and g respectively so that d; d build the quadratic with negative discriminant x2 þ mx þ n ¼ ðx  dÞðx  d Þ ¼ x2  2ReðdÞx þ jdj2 , where m and n are real constants with n ¼ jdj2 > 0 and g; g build the quadratic with negative discriminant x2 þ rx þ s ¼ ðx  gÞðx  g Þ ¼ x2  2ReðgÞx þ jgj2 , where r and s are real constants with s ¼ jgj2 > 0. Eq. (3.2) will then have the following factorisation

Bi xi ¼B6 ðx2 þ mx þ nÞðx2 þ rx þ sÞðx2 þ px þ qÞ ¼ B6 ½x6 þ ðm þ r þ pÞx5 þ ðq þ pm þ pr þ n þ mr þ sÞx4

dWðtÞ 6 /  gWðtÞ ¼ Cðt; WðtÞÞ 8 t 2 ð0; tb Þ: dt  nÞ ¼ /  g  n and  nð0Þ ¼ Wð0Þ ¼ W 0 . The Let us assume ddtn ¼ Cðt;  solution of this differential equation is

nðtÞ ¼ / ð1  egt Þ þ W 0 egt :

/ WðtÞ 6 nðtÞ ¼ ð1  egt Þ þ W 0 egt

þ ðqm þ qr þ ms þ nr þ pn þ pmr þ psÞx3 þ ðqn þ qmr þ qs þ pms þ pnrÞx2 þ ðqms þ qnr þ pnsÞx þ nsq ðA:5Þ where p and q are to be determined. By equating coefficients of like powers of x on the left and the right, we find

mþrþp¼

B5 ; B6 B4 ; B6

qðn þ mr þ sÞ þ ns þ pðms þ nrÞ ¼ qðms þ nrÞ þ pns ¼

Suppose t b < 1, then Wðtb Þ 6  nðt b Þ < 1, but then the solution exists uniquely for some interval ð0; tb Þ by the Picard–Lindelof Theorem. This contradicts the fact t b < 1. Thus WðtÞ must be bounded for all non-negative t and therefore sðtÞ is uniformly bounded in R3þ .

lim sup xðtÞ < 1:

B3 ; B6

t!þ1

Similarly, from the Eq. (2.2b), we observe that

B2 ; B6

    a1 xy a2 yz a1 6y y_ ¼ d1 y  h1 y2 þ   d1  h1 y : 1 þ b1 x 1 þ b2 y b1

B1 ; B6

Therefore, by using similar arguments, we have

B0 qns ¼ : B6

lim sup yðtÞ < t!þ1

From which we have

B5 B5  ðm þ rÞ ¼ þ 2ReðdÞ þ 2ReðgÞ; B6 B6 B0

q¼

ðA:7aÞ

ð2Þ x3

¼

½p 

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2  4q : 2

    a2 yz a2 d2 z  h2 z2 þ 6z  d2  h2 z : 1 þ b2 y b2

lim sup zðtÞ < ðA:8Þ

ð1Þ

By imposing conditions q < 0; p > 0, we obtain x3 > 0, which ensures the uniqueness and feasibility of the interior equilibrium. A.4. The proof of boundedness of the system (2.2) We define a positive definite function

WðtÞ ¼ x þ y þ z: From definition, WðtÞ is differentiable in some maximal interval ð0; t b Þ. For each g > 0 we have

dWðtÞ þ gWðtÞ ¼xðg þ ð1  xÞÞ þ yðg  d1  h1 yÞ dt 2 ðg þ 1Þ2 ðg  d1 Þ þ zðg  d2  h2 zÞ 6 þ 4 4h1

z_ ¼

Therefore, we have

The other two roots of Eq. (A.5) are thus found as

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2  4q ; 2

t!þ1

ðg  d2 Þ : 4h2

a2  b2 d2 ¼ ~z; b 2 h2

where ~z denotes an upper bound of zðtÞ which is positive if a2 > b2 d2 . A.6. Proof of stability and bifurcation analysis Before we show the proof we adopt the following notations. This system can be written in the form X_ ¼ FðXÞ, where X ¼ ðx; y; zÞT . The Jacobian matrix J  DFðXÞ of the system (2.2) at any arbitrary point ðx; y; zÞT is given by J ¼ ðmij Þ33 , where mij ’s are given by

a1 y a1 b1 xy þ ; 1 þ b1 x ð1 þ b1 xÞ2 a1 x a1 y ¼ ; m21 ¼ ; 2 1 þ b1 x ð1 þ b1 xÞ

m11 ¼ 1  2x  m12

a1 x a2 z a2 b2 yz  d1  þ  2h1 y; 1 þ b1 x 1 þ b2 y ð1 þ b2 yÞ2 a2 y a2 z a2 y ¼ ; m32 ¼ ; m33 ¼  d2  2h2 z: 2 1 þ b2 y 1 þ b2 y ð1 þ b2 yÞ

m22 ¼

2

þ

a1  b1 d1 ~; ¼y b 1 h1

~ denotes an upper bound of yðtÞ which is positive provided where y a1 > b1 d1 . Again, from the Eq. (2.2c), we observe that

B0 B6 ns

B6 jdj2 jgj2

¼

8t 2 ð0; t b Þ:

By comparison theorem ([31]), we obtain from Eq. (2.2a) that

qðm þ rÞ þ ðms þ nrÞ þ pðn þ mr þ sÞ ¼

½p þ

g

A.5. The proof of dissipativeness of the system (2.2)

q þ pðm þ rÞ þ mr þ s ¼

ð1Þ x3

Assume Cðt; yÞ ¼ /  gs. Then Cðt; yÞ satisfies Lipschtiz condition in R3þ . Clearly,

 nðtÞ is bounded in ð0; t b Þ. By comparison theorem (Birkhoff and Rota [31])

i¼0

¼

dWðtÞ þ gWðtÞ 6 / for each t 2 ð0; t b Þ: dt

g

6 X

p¼

Hence we can find / > 0 such that

m23

68

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

R3 ¼ ½ðx; y; zÞ 2 R3þ ; x > 0; y > 0; z P 0 and consider the scalar function L2 : R3 ! R defined by

A.6.1. Proof of Lemma 1 Proof of Lemma 1:

L2 ¼ (i) The eigenvalues of the Jacobian matrix J 1 of the system (2.2) a1 at the axial equilibrium point E1 are 1; 1þb  d1 and d2 . 1 Clearly, the equilibrium point E1 is locally asymptotically stable if a1 < d1 ð1 þ b1 Þ. (ii) J 1 has a zero eigenvalue iff det J 1 ¼ 0 which provides ½tc  a1 ¼ a1 1 . The other two eigenvalues are 1 and d2 . Let c1 and /1 are the eigenvectors corresponding to the eigenvalue 0 of the matrices J 1 and J T1 (the transpose of J 1 ), respecT T ½1 ½1 ½1 tively. Then we obtain that c1 ¼ ðh1 ; h2 ; 0Þ ; /1 ¼ ð0; d2 ; 0Þ , ½1 ½1 ½1 ½1 a1 where h1 ¼  1þb1 h2 , and h2 ; d2 are any two nonzero real ½tc  numbers. Since /T1 ½F a1 ðE1 ; a1 1 Þ ¼ 0, the system does not experience any saddle-node bifurcation. Again since ½tc  ½tc  /T1 ½DF a1 ðE1 ; a1 1 Þc1  – 0 and /T1 ½D2 FðE1 ; a1 1 Þðc1 ; c1 Þ – 0, then the system possesses a transcritical bifurcation ([32]), where D2 F is defined in Appendix A.8. (iii) Let R3 ¼ ½ðx; y; zÞ 2 R3þ ; x > 0; y P 0; z P 0 and consider the scalar function L1 : R3 ! R defined by

L1 ¼

1 1 ðx  1Þ2 þ y2 þ z2 þ y: 2 2

ðA:9Þ

The derivative of the Eq. (A.9) along the solution of the system (2.2) is given by

dL1 dx dy dz dy a1 x2 y ¼ðx  1Þ þ y þ z þ ¼ xð1  xÞ2  dt dt dt dt 1 þ b1 x dt 2a1 xy a1 xy2 a2 yz2 a2 yz2 2 3 þ  d1 y   h1 y þ þ 1 þ b1 x 1 þ b1 x 1 þ b2 y 1 þ b2 y a2 yz  d2 z2  h2 z3  d1 y   h1 y2 6 xð1  xÞ2 1 þ b2 y a1 x2 y 2a1 a1 a2 y2 z a2  þ y þ y2  d1 y2   h1 y3 þ z2 1 þ b1 x b1 1 þ b2 y b1 b2 a2 yz 2 2 3 2  d2 z  h2 z  d1 y   h1 y 6 xð1  xÞ 1 þ b2 y     2 a1 x y 2a1 a1  þ  d1 y þ  d1  h1 y2 1 þ b1 x b1 b1   a2 y2 z a a2 yz 2   h1 y3 þ  d2 z2  h2 z2  6 0; 1 þ b2 y 1 þ b2 y b2   2a1 a1 2a1 a2 ;  h1 ¼ ; d2 > ; ðA:10Þ if d1 > max b1 b1 b1 b2 also dLdt1 ¼ 0 when ðx; y; zÞ ¼ ð1; 0; 0Þ. Hence by Lyapunov–Lasalle’s invariance principle [33], E1 is globally asymptotically stable.

    x y þ k1 y  y2  y2 ln : x  x2  x2 ln x2 y2

The derivative of the Eq. (A.11) along the solution of the system (2.2) is given by

  dL2  x2  dx y dy ¼ 1 þ k1 1  2 dt x dt y dt   a1 y ¼ ðx  x2 Þ 1  x  þ k1 ðy 1 þ b1 x   a1 x  d1  h1 y :  y2 Þ 1 þ b1 x

c11 ¼ x2 þ c21 c31

a1 b1 x2 y2

;

c12 ¼ 

a1 x2 ; 1 þ b 1 x2

c13 ¼ 0; ð1 þ b1 x2 Þ a1 y2 a2 y2 ¼ ; c22 ¼ h1 y2 ; c23 ¼  ; 2 1 þ b2 y2 ð1 þ b1 x2 Þ a2 y2 ¼ 0; c32 ¼ 0; c33 ¼  d2 : 1 þ b2 y2 2

Now, tr J 2 ¼ c11 þ c22 þ c33 and det J2 ¼ ðc11 c22  c12 c21 Þc33 , where we notice that c22 < 0; c12 < 0; c21 > 0. Using Routh–Hurwitz criterion we obtain that the system is locally asymptotically stable around the equilibrium point E2 if tr J2 < 0 and det J 2 < 0, that is, d2 1 1Þ if c11 < 0 and c33 < 0, that is, if x2 < ðb2b and y2 < a2 b .ðiiÞ Let 1 2 d2

ðA:12Þ

Now at the equilibrium point E2 , we have

1 ¼ x2 þ

a1 y2 ; 1 þ b 1 x2

d1 ¼

a1 x2  h1 y 2 : 1 þ b1 x2

ðA:13Þ

The Eqs. (A.13) reduce the Eq. (A.12) to

  dL2 a1 y2 a1 y þ k1 ðy ¼ ðx  x2 Þ ðx  x2 Þ þ  dt 1 þ b1 x2 1 þ b1 x   a1 x a1 x2  y2 Þ   h1 ðy  y2 Þ 1 þ b 1 x 1 þ b 1 x2 ¼ ðx  x2 Þ2  k1 h1 ðy  y2 Þ2 þ

a1 b1 y2 ðx  x2 Þ2 ð1 þ b1 xÞð1 þ b1 x2 Þ



a1 b1 x2 ðx  x2 Þðy  y2 Þ a1 ðx  x2 Þðy  y2 Þ  ð1 þ b1 xÞð1 þ b1 x2 Þ ð1 þ b1 xÞð1 þ b1 x2 Þ

þ

k1 a1 ðx  x2 Þðy  y2 Þ ð1 þ b1 xÞð1 þ b1 x2 Þ

¼ ðx  x2 Þ2  k1 h1 ðy  y2 Þ2 þ þ

a1 b1 y2 ðx  x2 Þ2 ð1 þ b1 xÞð1 þ b1 x2 Þ

a1 ðx  x2 Þðy  y2 Þ ðk1  1  b1 x2 Þ: ð1 þ b1 xÞð1 þ b1 x2 Þ

Let us chose k1 ¼ ð1 þ b1 x2 Þ, then we obtain

dL2 a1 b1 y2 ðx  x2 Þ2 ¼ ðx  x2 Þ2  k1 h1 ðy  y2 Þ2 þ dt ð1 þ b1 xÞð1 þ b1 x2 Þ   a1 b1 y2 ðx  x2 Þ2 ¼ 1 þ ð1 þ b1 xÞð1 þ b1 x2 Þ  k1 h1 ðy  y2 Þ2   a 1 b1 y2 ðx  x2 Þ2  k1 h1 ðy  y2 Þ2 6 1 þ ð1 þ b1 x2 Þ 6 0 if a1 b1 y2 6 ð1 þ b1 x2 Þ; and

A.6.2. Proof of Lemma 2 Proof of Lemma 2: ðiÞ The Jacobian matrix J 2 of the system (2.2) around the boundary equilibrium point E2 is given by J 2 ¼ ðcij Þ33 , where cij are given by

ðA:11Þ

dL2 dt

¼ 0 when ðx; y; zÞ ¼ ðx2 ; y2 ; 0Þ:

ðA:14Þ

The proof follows from the Eq. (A.14) and Lyapunov–Lasalle’s invariance principle mentioned in [33]. ½sn  ðiiiÞ J 2 has a zero eigenvalue iff det J 2 ¼ 0 which gives a1 ¼ a1 2 , ½sn2  where a1 ¼ a1 . The other eigenvalues of J 2 is evaluated at ½sn  a1 ¼ a1 2 and one of them should be negative in order to get a saddle-node bifurcation. Let c2 and /2 are the eigenvectors corresponding to the eigenvalue 0 of the matrix J 2 and its transpose, ½2

½2

T

½2

½2

½2 T

½2

½2

respectively. We obtain that c2 ¼ ðh1 ; h2 ; 0Þ ; /2 ¼ ðd1 ; d2 ; d3 Þ , ½2

 cc23 33

½2 d2

½2

½2

½2

½2

h1 ¼  cc12 h ¼  cc22 h ; d1 ¼  cc21 d ¼  cc22 d ; d3 ¼ 11 2 21 2 11 2 12 2

where and

½2 ½2 h2 ; d2

½sn  /T2 ½F a1 ðE2 ; a1 2 Þ

¼

½2 ðd2

are 

any

½2 x2 y 2 d1 Þ 1þb 1 x2

two –

nonzero real

numbers.

½sn  0; /T2 ½D2 FðE2 ; a1 2 Þð 2 ;

c c2 Þ – 0.

Therefore, the system experiences a saddle-node bifurcation ½sn2 

around E2 at a1 ¼ a1

, [32]. To prove the transcritical bifurcation

69

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

let us assume c3 and /3 are the eigenvectors corresponding to the eigenvalue 0 of the matrix J 2 and its transpose respectively. We obtain ½3 h1

¼

½3 h2

½3 ½3 ½3 T ½3 ½3 ½3 T 3 ¼ ðh1 ; h2 ; h3 Þ ; /3 ¼ ðd1 ; d2 ; d3 Þ ½3 ½3 ½3 ½3 c22 ½3  c21 h2 ; h3 ¼ 0 and d1 ¼ 0 ¼ d2 ; h2 and

c

that  cc12 11

¼

where ½3

d2 are

any two real numbers. Now by the result of [32] we characterize the equilibrium point E2 for the parameter a2 . Since ½tc 

½tc 

/T3 ½DFðE2 ; a2 2 Þc3  – 0 and /T3 ½D2 FðE2 ; a2 2 Þðc3 ; c3 Þ – 0, the system ½tc 

attains a transcritical bifurcation around E2 at a2 ¼ a2 2 . ðiv Þ The characteristic equation of J 2 is given by

At the equilibrium point E3 of the system (2.2), we have

ðc33  kÞ½k2  ðc11 þ c22 Þk þ ðc11 c22  c12 c21 Þ ¼ 0

a1 y3 a1 x3 a2 z3 ; d1 ¼   h1 y 2 ; 1 þ b1 b3 1 þ b1 x3 1 þ b2 y3 a 2 y3 ¼  h2 z3 : 1 þ b2 y3

1 ¼ x3 þ

The eigenvalues are

 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ðc11 þ c22 Þ  ðc11 þ c22 Þ2  4ðc11 c22  c12 c21 Þ ; 2 a2 y2 ¼  d2 : 1 þ b2 y2

k1;2 ¼

k3

Using (A.20) in

Clearly, k3 is a real root, k1 and k2 are purely imaginary if and only if 2 ½hb  ½hb  1 x2 Þ there is an a1 ¼ a1 2 such that a1 2 ¼ ðx2 þh1by12xÞð1þb . Next we 2 y2 calculate

  dki b1 x2 y2 Re ¼ – 0; da1 a1 ¼a½hb2  ð1 þ b1 x2 Þ2 1

   dL3  x3  dx y dy z3  dz ¼ 1 þ k2 1  3 þ k3 1  dt x dt y dt z dt   a1 y þ k2 ðy ¼ ðx  x3 Þ 1  x  1 þ b1 x   a1 x a2 z  y3 Þ  d1   h1 y þ k3 ðz 1 þ b1 x 1 þ b2 y   a2 y  d2  h2 z  z3 Þ 1 þ b2 y

i ¼ 1; 2:

Therefore, the system enters into a Hopf-bifurcation around E2 at 2 ½hb  ½hb  1 x2 Þ a1 ¼ a1 2 , where a1 2 satisfies the equality a1 ¼ ðx2 þh1by12xÞð1þb . 2 y2

we obtain

a1 b1 y3 ðx  x3 Þ2 ð1 þ b1 xÞð1 þ b1 x3 Þ ðx  x3 Þðy  y3 Þ ða1  a1 b1 x3 þ a1 k2 Þ þ k2 þ ð1 þ b1 xÞð1 þ b1 x3 Þ

¼ ðx  x3 Þ2 þ

a2 b2 z3 ðy  y3 Þ2  k2 h1 ðy  y3 Þ2 ð1 þ b2 yÞð1 þ b2 y3 Þ a2 ðy  y3 Þðz  z3 Þ þ ðk2  k2 b2 y3 þ k3 Þ  k3 h2 ðz  z3 Þ2 : ð1 þ b2 yÞð1 þ b2 y3 Þ 

ðA:15Þ

Choose k2 ¼ ð1 þ b1 x3 Þ and k3 ¼ ð1 þ b1 x3 Þð1 þ b2 y3 Þ, we have

where Ai ’s are defined as

A1 ¼ m11  m22  m33 ;

ðA:20Þ

  dL3 a1 y a1 y3 þ k2 ðy ¼ ðx  x3 Þ ðx  x3 Þ  þ 1 þ b 1 x 1 þ b 1 x3 dt   a1 x a1 x3 a2 z a2 z3  y3 Þ   þ  h1 ðy  y3 Þ 1 þ b1 x 1 þ b1 x3 1 þ b2 y 1 þ b2 y3   a2 y a2 y3 þ k3 ðz  z3 Þ   h2 ðz  z3 Þ 1 þ b2 y 1 þ b2 y3

A.6.3. Proof of Lemma 3 Proof of Lemma 3: ðiÞ The Jacobian matrix J 3 of the system (2.2) around the interior equilibrium point E3 is given by J 3 ¼ ðmij Þ33 . Let ki ; i ¼ 1; 2; 3 be the roots of the characteristic equation of J 3 which is given by

k3 þ A1 k2 þ A2 k þ A3 ¼ 0;

dL3 , dt

d2

A3

¼ m12 m21 m33 þ m11 m23 m32  m11 m22 m33 ;

ðA:16Þ

A2 ¼ m22 m33 þ m11 m22 þ m11 m33  m12 m21  m23 m32 :

ðA:17Þ

  dL3 a1 b1 y3 ðx  x3 Þ2 ¼ 1 þ dt ð1 þ b1 xÞð1 þ b1 x3 Þ   a2 b2 z3 ðy  y3 Þ2 þ k2 h1 þ ð1 þ b2 yÞð1 þ b2 y3 Þ

 ðm22 Þ2 m33  m11 ðm22 Þ2  2m11 m22 m33

 k3 h2 ðz  z3 Þ2 ;   a1 b1 y3 ðx  x3 Þ2 6 1 þ ð1 þ b1 x3 Þ   a2 b2 z3 ðy  y3 Þ2  k3 h2 ðz  z3 Þ2 ; þ k2 h1 þ ð1 þ b2 y3 Þ

þ m22 m12 m21 þ m23 m32 m22  m22 ðm33 Þ2

6 0 if a1 b1 y3 6 ð1 þ b1 x3 Þ and a2 b2 z3

From Routh–Hurwitz criterion, E3 is locally asymptotically stable if and only if A1 > 0; A3 > 0 and A1 A2 > A3 . Now

A1 A2  A3 ¼ ðm11 Þ2 m22  ðm11 Þ2 m33 þ m11 m12 m21

2

 m11 ðm33 Þ þ m23 m32 m33 ;

ðA:18Þ

as m12 < 0; m21 > 0; m23 < 0; m32 > 0; m33 < 0, it is easy to check that A1 A2  A3 > 0 if m11 < 0,

2

m22 < 0, i.e if a1 b1 < ð1þby1 x3 Þ and 3

6 h1 ð1 þ b2 y3 Þ; and

dL3 ¼ 0 when ðx; y; zÞ dt

¼ ðx3 ; y3 ; z3 Þ:

ðA:21Þ

cally stable. ðiiÞ Let R3 ¼ ½ðx; y; zÞ 2 R3þ ; x > 0; y > 0; z > 0 and con-

The proof follows from the Eq. (A.21) and the Lyapunov–Lasalle’s invariance principle ([33]). ðiiiÞ J 3 has a zero eigenvalue if and only if det ðJ 3 Þ ¼ 0 which

sider the scalar function L3 : R3 ! R defined by

gives a1 ¼ a1

2

a2 b2 < ð1þbz23y3 Þ . Hence the equilibrium point E3 is locally asymptoti-

    x y þ k2 y  y3  y3 ln L3 ¼ x  x3  x3 ln x3 y3   z þ k3 z  z3  z3 ln z3

½sn3 

. The other eigenvalues of J 3 are evaluated at

½sn  a1 3

and one of them must be negative in order to get a sada1 ¼ dle-node bifurcation. Let c4 and /4 are the eigenvectors corresponding to the eigenvalue 0 of the matrix J 3 and its transpose,

ðA:19Þ

where k2 and k3 are two positive constants to be determined later. The derivative of the above Eq. (A.19) along the solution of the system (2.2) is given by

½4

½4

½4 T

½4

½4

½4 T

respectively. We obtain that c4 ¼ ðh1 ; h2 ; h3 Þ ; /4 ¼ ðd1 ; d2 ; d3 Þ where ½4

½4

½4

½4

½4

12 23 h1 ¼  m h ; h3 ¼  m h , m11 2 m33 2

½4

and

½4

½4

½4

21 d1 ¼  m d ; d3 ¼ m11 2

½4

23 m d , in which h2 and d2 are any two nonzero real numbers. m33 2

Since

½sn  /T4 ½F a1 ðE3 ; a1 3 Þ

½4

½4

x3 y 3 ¼ ðd2  d1 Þ 1þb – 1 x3

½sn3 

0; /T4 ½D2 FðE3 ; a1

Þ

70

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

ðc4 ; c4 Þ – 0, then the system experiences a saddle-node bifurca½sn3 

tion around E3 at a1 ¼ a1

([32]). Again, J 3 has a zero eigenvalue ½sn 

if and only if det ðJ 3 Þ ¼ 0 which gives a2 ¼ a2 3 . With a similar fashion we check that the system also experiences a saddle-node bifur½sn  a2 3 .

cation around E3 at a2 ¼ ðiv Þ The characteristic equation of J 3 is given by (A.15). By Routh– Hurwitz criterion, a set of necessary and sufficient conditions for all the roots of characteristic equation to have negative real parts is A1 > 0; A3 > 0 and A1 A2 > A3 . From the expression given in the Eqs. (A.16)–(A.18), we obtain that A1 > 0; A2 > 0; A3 > 0 if m11 < 0 and m22 < 0. Clearly, the characteristic Eq. (A.15) has two purely imagi½hb  nary roots if and only if A1 A2 ¼ A3 which provides a1 ¼ a1 3 . Since ½hb3  A2 > 0 at a1 ¼ a1 , there is an interval ða1  ; a1 þ Þ containing a1 for some  > 0 for which a1   > 0 such that A2 > 0 for a1 2 ða1  ; a1 þ Þ. Thus, for a1 2 ða1  ; a1 þ Þ, the characteristic equa½hb  tion can not have real positive roots. Now at a1 ¼ a1 3 , we have pffiffiffiffiffiffi 2 ðk þ A2 Þðk þ A1 Þ ¼ 0 which has three roots k1 ¼ i A2 ; k2 ¼ pffiffiffiffiffiffi i A2 ; k3 ¼ A1 . For a1 2 ða1  ; a1 þ Þ, the roots are in general of the form k1 ¼ aða1 Þ þ bða1 Þ; k2 ¼ aða1 Þ  bða1 Þ; k3 ¼ A1 ða1 Þ. Now we verify the transversality condition

  dki Re – 0; da1 a1 ¼a½hb3 

positive constants that stand for prey intrinsic growth rate and carrying capacity respectively. The system (E.1) has been nondimensionalised according to x ¼ KX ; y ¼ KY ; z ¼ KZ ; t ¼ cT as

dx a1 xy ¼ xð1  xÞ  ; dt 1 þ b1 x

ðA:24aÞ

dy a1 xy a2 yz ¼  d1 y  ; dt 1 þ b1 x 1 þ b2 y

ðA:24bÞ

dz a2 yz ¼  d2 z; dt 1 þ b2 y

ðA:24cÞ

xð0Þ > 0;

(i) (ii) (iii) (iv)

1

Substituting kj ¼ aða1 Þ þ bða1 Þ; j ¼ 1; 2 into the characteristic equation and calculating the derivative, we obtain

_ 1 Þ þ gða1 Þ ¼ 0; xða1 Þa_ ða1 Þ  /ða1 Þbða _ 1 Þ þ lða1 Þ ¼ 0; /ða1 Þa_ ða1 Þ þ xða1 Þbða

E3 ðx3 ; y3 ; z3 Þ, where qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

xða1 Þ ¼ 3a2 ða1 Þ þ 2A1 ða1 Þaða1 Þ þ A2 ða1 Þ  3b2 ða1 Þ; /ða1 Þ ¼ 6aða1 Þbða1 Þ þ A1 ða1 Þbða1 Þ; gða1 Þ ¼ a2 ða1 ÞA_ 1 ða1 Þ þ A_ 2 ða1 Þaða1 Þ þ A_ 3 ða1 Þ  A_ 1 ða1 Þb2 ða1 Þ; lða1 Þ ¼ 2aða1 Þbða1 ÞA_ 1 ða1 Þ þ A_ 2 bða1 Þ:

1 1

2b1

;

Þx3 d1 ð1þb2 y3 Þ 1 x3 Þ ; z3 ¼ ½ða1 b1 da12 ð1þb . y3 ¼ ð1x3 Þð1þb a1 1 x3 Þ

Lemma 4.

Proof. (i) The eigenvalues of the Jacobian matrix J 1 of the system (E.2) at the axial equilibrium point E1 is given by a1 1; 1þb  d1 ; d2 . Clearly, the equilibrium point E1 is locally 1 asymptotically stable if a1 < d1 ð1 þ b1 Þ. (ii) Let R3 ¼ ½ðx; y; zÞ 2 R3þ ; x > 0; y P 0; z P 0 and consider the scalar function L1 : R3 ! R defined by

Since /ða1 Þlða1 Þ þ xða1 Þgða1 Þ – 0, we have

j ¼ 1; 2;

a1 d2 1Þ 2 b2 d2

ðb1 1Þ2 4b1 ða

(i) The equilibrium point E1 ð1; 0; 0Þ of the system (E.2) is locally asymptotically stable iff a1 < d1 ð1 þ b1 Þ. 1 (ii) E1 ð1; 0; 0Þ is globally asymptotically stable if d1 > 2a ; d2 > ba22 . b1

where

k3 ða1 Þ

1

L1 ¼

¼ A1 ða1 Þ – 0: Hence, if m11 < 0 and m22 < 0, then the system exhibits a Hopf-bifurcation in the first octant leading to a family of solutions that bifurcates ½hb  from E3 for suitable values of a1 in a neighborhood of a1 3 . A.7. Study of the system in [16] The system studied by [16] is given by the following differential equations

  dX X M1 XY  ¼ cX 1  ; dT K A1 þ X

ðA:23aÞ

dY M1 XY M 2 YZ ¼  D1 Y  ; dT A1 þ X A2 þ Y

ðA:23bÞ

dZ M2 YZ ¼  D2 Z; dT A2 þ Y

ðA:23cÞ

Yð0Þ > 0;

Trivial equilibrium: E0 ð0; 0; 0Þ, Axial equilibrium: E1 ð1; 0; 0Þ, Boundary equilibrium: E2 ðx2 ; y2 ; 0Þ, Positive interior equilibrium:

d1 1 1 ; y2 ¼ aða1 bb1 d1ddÞ21 ; x3 ¼ b2b þ x2 ¼ a1 b 1 d1 1

1

Xð0Þ > 0;

zð0Þ > 0;

where the dimensionless parameters are given by a1 ¼ McA11K ; b1 ¼ AK1 ; d1 ¼ Dc1 ; a2 ¼ McA22K ; b2 ¼ AK2 ; d2 ¼ Dc2 . The system (E.2) has the following positive equilibria:

i ¼ 1; 2:

  dkj /l þ xg Re ¼ 2 – 0; da1 a1 ¼a½hb3  / þ x2

yð0Þ > 0;

Zð0Þ > 0;

where XðTÞ; YðTÞ; ZðTÞ represent the population density of prey, predator and top-predator at time T respectively while the parameters Ai ; Mi ; Di (i ¼ 1; 2) designate half-saturation constant, maximal predation rate and predator’s death rate, respectively and c; K are

1 1 ðx  1Þ2 þ y2 þ z2 þ y 2 2

ðA:25Þ

The derivative of the above Eq. (A.25) along with the solution of the system (E.2) is given by

dL1 ¼ ðx  1Þx_ þ yy_ þ zz_ þ y_ dt a1 x2 y 2a1 xy a1 xy2 ¼ xð1  xÞ2  þ þ  d1 y2 1 þ b1 x 1 þ b1 x 1 þ b1 x a2 yz2 a2 yz2 a2 yz  þ  d2 z2  d1 y  1 þ b2 y 1 þ b2 y 1 þ b2 y a1 x2 y 2a1 a1 þ y þ y2  d1 y2 1 þ b1 x b1 b1 a2 y2 z a2 2 a2 yz  þ z  d2 z2  d1 y  1 þ b2 y b2 1 þ b2 y     a1 x2 y 2a1 a1 þ  d1 y þ  d1 y2 6 xð1  xÞ2  1 þ b1 x b1 b1   a2 y2 z a2 a yz 2  þ  d2 z2  1 þ b2 y 1 þ b2 y b2   2a1 a1 2a1 a2 ¼ ; ; d2 > 6 0; if d1 > max b1 b1 b1 b2 6 xð1  xÞ2 

ðA:26Þ

Also dLdt1 ¼ 0 when ðx; y; zÞ ¼ ð1; 0; 0Þ. Hence by Lyapunov–Lasalle’s invariance principle [33], E1 is globally asymptotically stable.

M. Haque et al. / Mathematical Biosciences 246 (2013) 55–71

Lemma 5. The equilibrium point E3 is locally asymptotically stable if 2

2

a1 b1 < ð1þby1 x3 Þ and ½x3 ð1 þ b1 x3 Þ þ a1 b1 x3 y3 a2 b2 y3 z3  a21 x3 y3 ð1þ 3

2

b2 y3 Þ2  a2 y3 z3 ð1 þ b1 x3 Þ < 0. Proof. The characteristic equation of the Jacobian matrix J 3 ¼ ½mij ; i; j ¼ 1; 2; 3 of the system (E.2) at the equilibrium point E3 with m13 ¼ 0; m31 ¼ 0; m33 ¼ 0 is given by

ðm11  kÞðm22  kÞðkÞ  m12 m21 ðkÞ þ m23 m32 ðm11  kÞ ¼ 0 a1 b1 x3 y3 a1 y3 a1 x3 where m11 ¼ x3 þ ð1þb m12 ¼  1þb ; m21 ¼ ð1þb 2, 2, 1 x3 1 x3 Þ a y 1 x3 Þ a2 b2 y3 z3 a2 z3 2 3 m22 ¼ ð1þb y Þ2  h1 y3 , m23 ¼  1þb2 y3 ; m32 ¼ ð1þb y Þ2 . 2 3

2 3

Then the roots ki ; i ¼ 1; 2; 3 satisfy the equation k3 þ A1 k2 þ A2 kþ A3 ¼ 0 where

A1 ¼ m11  m22 ;

ðA:27Þ

A2 ¼ m11 m22  m12 m21  m23 m32 ;

ðA:28Þ

A3 ¼ m11 m23 m32 :

ðA:29Þ

From Routh–Herwitz criterion, the equilibrium point E3 is locally asymptotically stable if and only if A1 > 0; A3 > 0 and A1 A2 > A3 . A1 A2  A3 ¼ ðm11 Þ2 m22 þ m11 m23 m32 þ m11 m12 m21 

Now, 2

m11 ðm22 Þ þ m22 m12 m21 þ m23 m32 m22 > 0,

if

m11 < 0,

m12 m21  m11 m22 < 0, that is, if m11 < 0 ) a1 b1 <

m23 m32 þ

ð1þb1 x3 Þ2 y3

and

2

m23 m32 þ m12 m21  m11 m22 < 0 ) ½x3 ð1 þ b1 x3 Þ þ a1 b1 x3 y3 a2 b2 y3 2

2

z3  a21 x3 y3 ð1 þ b2 y3 Þ  a2 y3 z3 ð1 þ b1 x3 Þ < 0. Hence, the equilibrium point E3 is locally asymptotically stable. A.8. Notations In the current paper, we use the same notations that used in the Book of [34], D2 Fðx; kÞ is defined by 2 @2 F

3 @2 F 1 @2 F 1 @2 F 1 2 @2 F1 @2 F1 @2 F1 2 1 2 2 u1 þ @x@y u1 u2 þ @y@x u2 u1 þ @y2 u2 þ @y@z u2 u3 þ @z@y u3 u2 þ @z2 u3 6 @x 7 2 2 2 2 2 2 2 6 7 F2 F2 F2 F2 D2 Fðx;kÞðU;UÞ ¼ 6 @@xF22 u21 þ @@x@y u1 u2 þ @@y@x u2 u1 þ @@yF22 u22 þ @@y@z u2 u3 þ @@z@y u3 u2 þ @@zF22 u23 7 4 5 @2 F3 2 @2 F 3 @2 F 3 @2 F 3 2 @2 F3 @2 F3 @2 F3 2 u þ @x@y u1 u2 þ @y@x u2 u1 þ @y2 u2 þ @y@z u2 u3 þ @z@y u3 u2 þ @z2 u3 @x2 1

where U ¼ ðu1 ; u2 ; u3 ÞT is the eign vector corresponding to a zero eignvalue of the matrix, J (defined in Appendix A.6.1), and k is the bifurcation parameter. References [1] Y.A. Kuznetsov, O. De Feo, S. Rinaldi, Belyakov homoclinic bifurcations in a tritrophic food chain model, SIAM Journal on Applied Mathematics (2001) 462. [2] S. Lv, M. Zhao, The dynamic complexity of a three species food chain model, Chaos, Solitons & Fractals 37 (5) (2008) 1469.

71

[3] Y. Chen, J. Yu, C. Sun, Stability and hopf bifurcation analysis in a threelevel food chain system with delay, Chaos, Solitons & Fractals 31 (3) (2007) 683. [4] B.W. Kooi, M.P. Boer, S. Kooijman, Consequences of population models for the dynamics of food chains, Mathematical Biosciences 153 (2) (1998) 99. [5] B.W. Kooi, L.D.J. Kuijper, M.P. Boer, S. Kooijman, Numerical bifurcation analysis of a tri-trophic food web with omnivory, Mathematical Biosciences 177 (2002) 201. [6] M.P. Boer, B.W. Kooi, S. Kooijman, Multiple attractors and boundary crises in a tri-trophic food chain, Mathematical Biosciences 169 (2) (2001) 109. [7] M.P. Boer, B.W. Kooi, S. Kooijman, Homoclinic and heteroclinic orbits to a cycle in a tri-trophic food chain, Journal of Mathematical Biology 39 (1) (1999) 19. [8] A. Hastings, T. Powell, Chaos in a three-species food chain, Ecology (1991) 896. [9] G.R. Huxel, K. McCann, Food web stability: the influence of trophic flows across habitats, The American Naturalist 152 (3) (1998) 460. [10] P.A. Abrams, J.D. Roth, The effects of enrichment of three-species food chains with nonlinear functional responses, Ecology (1994) 1118. [11] O. De Feo, S. Rinaldi, Singular homoclinic bifurcations in tritrophic food chains, Mathematical Biosciences 148 (1) (1998) 7. [12] J.N. Eisenberg, D.R. Maszle, The structural stability of a three-species food chain model, Journal of Theoretical Biology 176 (4) (1995) 501. [13] A. Klebanoff, A. Hastings, Chaos in three species food chains, Journal of Mathematical Biology 32 (5) (1994) 427. [14] K. Mccann, P. Yodzis, Bifurcation structure of a three-species food-chain model, Theoretical Population Biology 48 (2) (1995) 93. [15] Y.A. Kuznetsov, S. Rinaldi, Remarks on food chain dynamics, Mathematical Biosciences 134 (1) (1996) 1. [16] C.H. Chiu, S.B. Hsu, Extinction of top-predator in a three-level food-chain model, Journal of Mathematical Biology 37 (4) (1998) 372. [17] E.C. Pielou, An Introduction to Mathematical Ecology, John Wiley and Sons, 1969. [18] A.D. Bazykin, A.I. Khibnik, B. Krauskopf, Nonlinear Dynamics of Interacting Populations, World Scientific Pub Co Inc, 1998. [19] M. Haque, Ratio-dependent predator-prey models of interacting populations, Bulletin of Mathematical Biology 71 (2) (2009) 430. [20] M. Haque, A detailed study of the Beddington–DeAngelis predator–prey model, Mathematical Bioscience 234 (1) (2011) 1. [21] M. Haque, Existence of complex patterns in the Beddington–DeAngelis predator–prey model, Mathematical Bioscience 239 (2012) 179. [22] Y.A. Kuznetsov, Elements of Applied Bifurcation Theory, Springer, 1998. [23] B.D. Hassard, N.D. Kazarinoff, Y.H. Wan, Theory and Applications of Hopf Bifurcation, Cambridge University Press, 1981. [24] H.I. Freedman, Deterministic Mathematical Models in Population Ecology, Marcel Dekker Inc, 1980. [25] R.M. May, Stability and Complexity in Model Ecosystems, Princeton Univ Pr, 2001. [26] P.A. Abrams, J. Roth, The responses of unstable food chains to enrichment, Evolutionary Ecology 8 (2) (1994) 150. [27] G.F. Fussmann, S.P. Ellner, K.W. Shertzer, N.G. Hairston Jr., Crossing the Hopf bifurcation in a live predator–prey, system, 2000. [28] P. Hogeweg, B. Hesper, Interactive instruction on population interactions, Computers in Biology and Medicine 8 (4) (1978) 319. [29] V. Rai, R. Sreenivasan, Period-doubling bifurcations leading to chaos in a model food chain, Ecological Modelling 69 (1) (1993) 63. [30] M. Scheffer, Should we expect strange attractors behind plankton dynamics– and if so, should we bother?, Journal of Plankton Research 13 (6) (1991) 1291 [31] G. Birkhoff, G.C. Rota, Ordinary Differential Equations, Ginn, Boston, 1978. [32] J. Sotomayor, Generic Bifurcations of Dynamical systems. Dynamical Systems, Academic Press, San Diego, 1973. pp. 549–560. [33] J.K. Hale, Ordinary Differential Equations, Wiley-Interscience, New York, 1969. [34] L. Perko, Differential Equations and Dynamical Systems, third ed., Springer, 2001.