- Email: [email protected]
Dynamical behavior of three species predator–prey system with mutual support between non refuge prey
Accepted Manuscript Dynamical behavior of three species predator–prey system with mutual support between non refuge prey D. Pal, P. Santra, G.S. Mahapatra PII:
S2405-9854(17)30002-2
DOI:
10.1016/j.egg.2017.05.001
Reference:
EGG 9
To appear in:
Ecological Genetics and Genomics
Received Date: 26 January 2017 Revised Date:
14 April 2017
Accepted Date: 15 May 2017
Please cite this article as: D. Pal, P. Santra, G.S. Mahapatra, Dynamical behavior of three species predator–prey system with mutual support between non refuge prey, Ecological Genetics and Genomics (2017), doi: 10.1016/j.egg.2017.05.001. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Dynamical Behavior of Three Species Predator–Prey System with Mutual Support Between Non Refuge Prey D. Pal1∗ , P. Santra2 and G. S. Mahapatra3 2
Chandrahati Dilip Kumar High School(H.S.), Chandrahati 712504, West Bengal, India.
RI PT
1
Department of Mathematics, Abada New Set Up Upper Primary School, Howrah, West Bengal, India 3
Department of Mathematics, National Institute of Technology Puducherry, Karaikal-609609, India. Tel. +919433135327, Fax: 04368231665
E-mail: [email protected]
SC
∗ Correspondence
1
M AN U
Abstract: In this paper, we study the effect of refuge and harvesting in a three species model of two prey and one predator, where the two prey refuge groups are helping each other and shielding from the predator group. The stability theory of nonlinear differential equations is used to analyze the model. Our study exhibits that, both the refuge and harvesting play a significant role on the stability of the system. Along with the local stability of the system, the global stability is also studied. Finally, numerical simulations are given to support the analytical conclusion. Keywords: Prey-Predator; Harvesting; Refuge; Stability; Mutualism.
Introduction
AC C
EP
TE D
In population dynamics, the relationship between prey-predator plays an importance for universal existence. Lotka [1] and Volterra [2] first studied the relationship between prey and predator. During the last few decades, several research articles have been published on the prey-predator system incorporating refuge concept and mutualism. The study of stability and harvesting of species including refuge is growing in importance. Due to the mutual support between the species, the study of harvesting of the prey-predator system is enhanced based on the study of refuge and mutualness of prey-predator system, as per the diverse nature of our ecological system. Recently, several researchers (Ding et. al. [3], Pal et al. ([4]-[6]), Pal and Mahapatra [7], Gupta and Chandra [8], Murphy and Smith [9], Palma and Olivares [10], Luo et al. [11]) concentrated their research work in the field of harvesting management of renewable resources. We have studied the refuge and harvesting of two prey and one predator model. Rai et al. [12] analyzed the model of mutualism in predator-prey and competitive systems in three species environment. Ruxton [13] studied the stability of prey-predator model using short term refuge. Huang et. al. [14] presented a prey–predator model with Holling type III response function incorporating a prey refuge. Frisvold and Reeves [15] studied the costs and benefits of refuge requirements. Ma et al. [16] discussed the effects of prey refuges on a prey-predator model. Ross and J´ozsef [17] studied a predator–prey refuge system. Chen et al. [18] studied qualitative analysis of a predator-prey model with Holling type II functional response incorporating a constant prey refuge. Ji and Wu [19] presented qualitative analysis of a prey-predator model with constant rate of prey harvesting incorporating a constant prey refuge. Tao et al.[20] discussed the effects of prey refuge on a harvested predator–prey model with generalized functional response. Wu et al. [21] derived a sufficient condition for the existence of at least one positive periodic solution of a generalized prey–predator model with harvesting term. Mostafa and Mahmoud [22] presented periodic solutions of a prey-predator system with Beddington–DeAngelis functional response [23] on time scales. Jia et al. [24] studied the existence of positive solutions for a prey-predator model with refuge and diffusion. Leard and Rebaza [25] presented a prey-predator model with continuous threshold harvesting and formulated an alternative form of threshold harvesting. Chen et al. [26] studied the stability property for the predator-free equilibrium point of prey-predator systems with a class of functional response [27] and prey refuges. In this paper, we consider a model which consists of two prey groups and one predator group. Prey groups help each other against the predator and few prey of each prey groups are the refuge. Here we consider an interesting situation, i.e., when a predator attacks the 1st prey group then some prey of the 1st prey group will be refuge, but the total population of 2nd prey group help 1st prey group against the predator and when a predator attacks the
1
ACCEPTED MANUSCRIPT
2
RI PT
2nd prey group then some prey of the 2nd prey group will be refuge but the total population of 1st prey groups help the 2nd prey groups against the predator. In this model, we also consider that the predator species is harvested. We have constructed the model, motivated by the work of Elettreby and El-Metwally [28] in which they have studied the local and global stability of multi-team prey-predator model. The rest of the paper is organized as follows: In section 2, we present the mathematical formulation of three species prey-predator system. Positivity of the model is presented in Section 3. Section 4 present the possible steady states of the proposed predator-pray model. In section 5, we study the local stability as well as global stability analysis of the equilibrium points. The Numerical simulation is presented in section 6. Finally, Section 7 gives concluding remark of our work.
Mathematical model of proposed predator–prey system
M AN U
SC
In our proposed model, two teams of prey with densities x(t) and y(t) respectively, interact with one team of predator with densities z(t). The assumptions of the prey-predator model are as follows: (i) Each team of preys grows logistically in absence of any predator, (ii) Due to effect of predation, the prey growth rate is reduced by a quantity proportional to prey and predator population, (iii) The non refuge teams of prey help each other against the predator, (iv) In absence of any prey for sustenance, the predator’s death rate results is in inverse decay, i.e., −cz 2 , (v) Prey’s contribution to predator growth rate is proportional to available prey and size of predator population, (vi) Predator species is harvested with the harvesting effort E, (vii) Some preys of each team are refuge, i.e., αx and βy, (viii) When the team of 2nd prey helps first prey against predator then not only the total populations of 2nd team involved in this interaction but also non refuge population of 1st team involve in this interaction. Same situation occurs for the first prey support to the second prey. The rate of primary resources available for the prey to feed plays an important role in this kind of problems ([29]-[33]). Using the above assumptions, the following predator-prey model is proposed:
TE D
dx = ax(1 − x) − (1 − α)xz + (1 − α)xyz dt dy = by(1 − y) − (1 − β)yz + (1 − β)xyz dt dz = −cz 2 + d(1 − α)xz + e(1 − β)yz − Ez dt
(1)
EP
where the coefficients a, b, c, d and e are positive constants. With initial densities and x(0) > 0, y(0) > 0 and z(0) > 0
(2)
3
AC C
It is clear that the two teams of preys help each other, e.g., in foraging and in early warning against predation. Note that the mutual support occurs only in the presence of the predator.
Positivity of the Proposed System
Theorem 1 Every solution of the system (1) with initial conditions (2) exists in the interval [0, ∞) and x(0) > 0, y(0) > 0 and z(0) > 0 for all t ≥ 0. Proof. Since the right-hand side of system (1) is completely continuous and locally Lipschitzian on C, the solution (x(t), y(t), z(t)) of (1) with initial conditions (2) exist and is unique on [0, ξ), where 0 < ξ < ∞ [34]. From system
2
ACCEPTED MANUSCRIPT (1) with initial conditions (2), we have t Z x (t) = x (0) exp {a(1 − x (θ)) − (1 − α)z (θ) + (1 − α)y (θ) z (θ)} dθ > 0 0t Z y (t) = y (0) exp {b(1 − y (θ)) − (1 − β)z (θ) + (1 − β)x (θ) z (θ)} dθ > 0
RI PT
0t Z z (t) = z (0) exp {−cz (θ) + d(1 − α)x (θ) + e(1 − β)y (θ) − E} dθ > 0 0
which completes the proof.
4
Steady States Analysis of the Proposed Predator-Prey System
M AN U
SC
The possible steady states of the predator-prey dynamical system (1) are as follows: G0 (0, 0, − Ec ) (this equilibrium point does not biologically relevant as − Ec < 0), G1 (0, 0, 0), G2 (1, 0, 0), G3 (0, 1, 0), G4 (1, 1, 0) and the other equilibrium points are be(1−β)−Eb cb+E(1−β) G5 (0, y5 , z5 ): where y5 = cb+e(1−β) 2 > 0 and z5 = cb+e(1−β)2 . The equilibrium point G5 exists if e(1 − β) > E. ac+E(1−α) ad(1−α)−aE ac+d(1−α)2 > 0 and z6 = ac+d(1−α)2 . The equilibrium point G6 exists if d(1 − α) > E. . The equilibrium point G6 exists if d(1 − α) + e(1 − β) > E. G7 (1, 1, z7 ): where z7 = d(1−α)+e(1−β)−E c q q 2 ab ab bc+e(1−β) ac+d(1−α)2 }+a{E−d(1−α)} }+b{E−e(1−β)} (1−α)(1−β) { (1−α)(1−β) { q q G8 (x8 , y8 , z8 ): where x8 = , y = and 8 ab ab 2 2 e(1−β) d(1−α) +bd(1−α) +ae(1−β) (1−α)(1−β) (1−α)(1−β)
G6 (x6 , 0, z6 ): where x6 =
z8 =
ab (1−α)(1−β)
> 0. The equilibrium point G8 exists if
s ab ab 2 bc + e(1 − β) +b {E − e(1 − β)} > 0 and ac + d(1 − α)2 +a {E − d(1 − α)} > 0 (1 − α)(1 − β) (1 − α)(1 − β) −
q
TE D
s
q
ab (1−α)(1−β)
{bc+e(1−β)2 }+b{E−e(1−β)}
−
q
ab (1−α)(1−β)
q G9 (x9 , y9 , z9 ): where x9 = , y9 = ab −e(1−β)2 (1−α)(1−β) +bd(1−α) q ab < 0.The equilibrium point G9 does not exist since z9 < 0. z9 = − (1−α)(1−β)
−d(1−α)2
q
ab +ae(1−β) (1−α)(1−β)
and
Stability Analysis
EP
5
{ac+d(1−α)2 }+a{E−d(1−α)}
5.1
AC C
In this section, we present local as well as global stability of the equilibrium points of the proposed predator-prey System.
Local Stability analysis
Here, we study the local stability analysis of our proposed model. It is clear that the equilibrium points G1 (0, 0, 0), G2 (1, 0, 0), G3 (0, 1, 0) are unstable. Now we analyze the local stability at the other equilibrium points of the system based on the following Theorems. The variational matrix M of the system (1) corresponding to different equilibrium points is given by M11 (1 − α)xz (1 − α) (y − 1) x M22 (1 − β) (x − 1) y M = (1 − β)yz d(1 − α)z e(1 − β)z M33
Here,
M11 = a − 2ax − (1 − α)z + (1 − α)yz, M22 = b − 2by − (1 − β)z + (1 − β)xz, M33 = −2cz + d(1 − α)x + e(1 − β)y − E,
3
ACCEPTED MANUSCRIPT Theorem 2 The equilibrium point G4 (1, 1, 0) is stable if d(1 − α) + e(1 − β) < E. Proof. The variational matrix M corresponding to the equilibrium point G4 (1, 1, 0) has eigenvalues −a, −b and d(1 − α) + e(1 − β) − E respectively. Since −a < 0 and −b < 0, the equilibrium point G4 (1, 1, 0) is stable if d(1 − α) + e(1 − β) − E < 0 i.e., d(1 − α) + e(1 − β) < E. a Theorem 3 The equilibrium point G5 (0, y5 , z5 ) is locally asymptotically stable if (y5 − 1)z5 < − (1−α) .
SC
The characteristic equation of the above matrix MG5 (0,y5 ,z5 ) is given by det MG5 (0,y5 ,z5 ) − λI = 0
RI PT
Proof. The variational matrix at the equilibrium point G5 (0, y5 , z5 ) is given by a + (1 − α)(y5 − 1)z5 0 0 (1 − β)y5 z5 −by5 −(1 − β)y5 MG5 (0,y5 ,z5 ) = d(1 − α)z5 e(1 − β)z5 −cz5
(3)
One root of the characteristic equation (3) i.e. one eigenvalue is given by a + (1 − α)(y5 − 1)z5 . This eigenvalue is a . The other two eigenvalues are the root of the quadratic equation negative if (y5 − 1)z5 < − (1−α) (4)
M AN U
λ2 + λ(cz5 + by5 ) + y5 z5 {cb + e(1 − β)2 } = 0
Now, sum of the roots of the equation (4) is −(cz5 + by5 ), which is negative. The product of the roots of the equation (4) is y5 z5 {cb + e(1 − β)2 }, which is always positive. Hence, the roots of the quadratic equation (4) are real and negative or complex conjugates with negative real parts. Therefore, G5 (0, y5 , z5 ) is locally asymptotically stable if a . (y5 − 1)z5 < − (1−α) b . Theorem 4 The equilibrium point G6 (x6 , 0, z6 ) is locally asymptotically stable if (x6 − 1)z6 < − (1−β)
TE D
Proof. The variational matrix at the equilibrium point G6 (x6 , 0, z6 ) is given by −ax6 (1 − α)x6 z6 −(1 − α)x6 MG6 (x6 ,0,z6 ) = 0 b + (1 − β)(x6 − 1)z6 0 d(1 − α)z6 e(1 − β)z6 −cz6
EP
The characteristic equation of the above matrix MG6 (x6 ,0,z6 ) is given by det MG6 (x6 ,0,z6 ) − λI = 0
(5)
One root of the characteristic equation (5) is given by b + (1 − β)(x6 − 1)z6 . This eigenvalue is negative if (x6 − 1)z6 < b .The other two eigenvalues are the root of the quadratic equation − (1−β)
AC C
λ2 + λ(cz6 + ax6 ) + x6 z6 {ca + d(1 − α)2 } = 0
(6)
The sum of the roots of the equation (6) is−(cz6 + ax6 ) < 0. The product of the roots of the equation (6) is x6 z6 {ca + d(1 − α)2 }, which is always positive. Hence, the roots of the above quadratic equation (6) are real and negative or complex conjugates with negative real parts. Therefore, G6 (x6 , 0, z6 ) is locally asymptotically stable if b (x6 − 1)z6 < − (1−β) . Theorem 5 The equilibrium point G7 (1, 1, z7 ) is locally asymptotically stable if ab > (1 − α)(1 − β)z72 . Proof. The variational matrix at the equilibrium point G7 (1, 1, z7) is given by −a (1 − α)z7 0 −b 0 MG7 (1,1,z7 ) = (1 − β)z7 d(1 − α)z7 e(1 − β)z7 −cz7
The characteristic equation of the matrix MG6 (x6 ,0,z6 ) is given by
det MG7 (1,1,z7 ) − λI = 0 4
(7)
ACCEPTED MANUSCRIPT One root of the characteristic equation (7) is −cz7 , which is negative. The other two eigenvalues are the root of the quadratic equation λ2 + λ(a + b) + {ab − (1 − α)(1 − β)z72 } = 0 (8) The sum of the roots of the equation (8) is −(a + b) < 0. The product of the roots of the equation (8) is{ab − (1 − α)(1 − β)z72 }. The product of the roots is positive if ab > (1 − α)(1 − β)z72 . Therefore, if ab > (1 − α)(1 − β)z72 then the roots of the above quadratic equation (8) are real and negative or complex conjugates with negative real parts. Hence, G7 (1, 1, z7 ) is locally asymptotically stable if ab > (1 − α)(1 − β)z72 .
RI PT
Theorem 6 The equilibrium point G8 (x8 , y8 , z8 ) is locally asymptotically stable if {ac + d(1 − α)2 }x8 + {bc + e(1 − β)2 }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2(y8 − 1)} + (1 − α)(1 − β)z8{d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0.
SC
Proof. The variational matrix at the equilibrium point G8 (x8 , y8 , z8 ) is given by −ax8 (1 − α)x8 z8 (1 − α)x8 (y8 − 1) −by8 (1 − β)y8 (x8 − 1) MG8 (x8 ,y8 ,z8 ) = (1 − β)y8 z8 d(1 − α)z8 e(1 − β)z8 −cz8 The characteristic equation of the the above matrix MG8 (x8 ,y8 ,z8 ) is given by
(9)
M AN U
f (λ) = 0
Global stability analysis of predator–prey system
EP
5.2
TE D
where f (λ) = λ3 + λ2 ω1 + λω2 + ω3 , ω1 = (ax8 + by8 + cz8 ), ω2 = [{ac + d(1 − α)2 }x8 z8 + {bc + e(1 − β)2 }y8 z8 − {e(1 − β)2 + d(1 − α)2 }x8 y8 z8 ], ω3 = −[x8 y8 z8 {ae(1 − β)2 (x8 − 1) + bd(1 − α)2 (y8 − 1)} + (1 − α)(1 − β)x8 y8 z82 {d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}]. Now ω1 > 0, ω2 > 0 if {ac + d(1 − α)2 }x8 + {bc + e(1 − β)2 }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and ω3 > 0 if [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2 (y8 − 1)} + (1 − α)(1 − β)z8 {d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0. The roots λi (i = 1, 2, 3), of the polynomial equation (9)satisfies the following conditions (i) λ1 + λ2 + λ3 = −ω1 < 0 (ii)λ1 λ2 + λ2 λ3 + λ3 λ1 = ω2 > 0 (iii)λ1 λ2 λ3 = −ω3 < 0. Since ω1 , ω2 and ω3 are all positive, so there is no change of sign in the polynomial f (λ) for λ > 0. Therefore, f (λ) = 0 does not possess any positive root. Again in f (−λ), there are three changes of sign. Therefore, the nature of the roots f (λ) = 0 has three possibilities such as (i) no positive root (ii) one negative root and two complex roots with negative real parts or (iii) three negative roots. Therefore, the equilibrium point G8 (x8 , y8 , z8 ) is locally asymptotically stable if {ac + d(1 − α)2 }x8 + {bc + e(1 − 2 β) }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2(y8 − 1)} + (1 − α)(1 − β)z8{d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0.
AC C
Here, we study the global stability of the interior equilibrium points G7 (1, 1, z7 ) and G8 (x8 , y8 , z8 ) by considering suitable Lyapunov function. Theorem 7 Equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable. Proof. In order to test whether the equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable or not, we construct a suitable Lyapunov function as follows z V = [(x − 1) − log x] + l1 [(y − 1) − log y] + l2 (z − z7 ) − z7 log (10) z7 It is obvious that V is positive definite. Differentiating both sides of (10) with respect to time t along the solutions of (1) we get, dV dt
⇒
= dV dt
dy z−z7 dz + l1 y−1 y dt + l2 z dt 2 = −a(x − 1) + (1 − α)(x − 1)(yz − z7 ) − l1 b(y − 1)2 + l1 (1 − β)(y − 1)(xz − z7 ) − l2 c(z − z7 )2
x−1 dx x dt
The above equation can be written in real quadratic form dV = −X T P X dt 5
ACCEPTED MANUSCRIPT where X T = [(x − 1) , (y − 1) , (z − z7 )] and 7) a −(1 − α) (yz−z x−1 P = 0 bl1 0 0
0 7) −l1 (1 − β) (xz−z y−1 cl2
Theorem 8 Equilibrium point G8 (x8 , y8 , z8 ) is globally asymptotically stable.
RI PT
Now dV dt < 0, if the matrix P is positive definite. The matrix P is positive definite as all the principal minors are positive and hence dV dt < 0. Therefore, the interior equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable.
Proof. we construct a suitable Lyapunov function as follows: x y z W = (x − x8 ) − x8 log + l1 (y − y8 ) − y8 log + l2 (z − z8 ) − z8 log x8 y8 z8
(11)
dW dt
dW dt
z−z8 dz 8 dy + l1 y−y y dt + l2 z dt = −a(x − x8 )2 + (1 − α)(x − x8 )(yz − y8 z8 ) − l1 b(y − y8 )2 + l1 (1 − β)(y − y8 )(xz − x8 z8 ) − l2 c(z − z8 )2
x−x8 dx x dt
M AN U
⇒
=
SC
It is obvious that W is positive definite. Differentiating both sides of (11) with respect to time t along the solutions of (1) we get,
The above equation can be written in real quadratic form
dW = − − Y T QY dt
TE D
where Y T = [(x − x8 ) , (y − y8 ) , (z − z8 )] and 8 z8 ) a −(1 − α) (yz−y x−x8 Q= 0 bl1 0 0
0 8 z8 ) −l1 (1 − β) (xz−x y−y8 cl2
Now dV dt < 0, if the matrix Q is positive definite. The matrix Q is positive definite as all the principal minors are positive and hence dW dt < 0. Therefore, the equilibrium point G8 (x8 , y8 , z8 ) is globally asymptotically stable.
Numerical Simulation
EP
6
AC C
Analytical studies can never be completed without numerical verification of the derived results. In this section we present computer simulation of some solutions of the systems (1). Beside verification of our analytical findings, these numerical solutions are very important from practical point of view. For the purpose of simulation experiments we mainly use the software MATLAB R2008a and Wolfram Mathematica 8. As the problem is not a case study for a particular species, here some hypothetical data are taken for the sole purpose of illustrating the analytical results of the previous sections. Moreover, it may be noted that as the parameters of the model are not based on real world observations, the main features described by the simulations presented in this section should be considered from a qualitative, rather than a quantitative point of view. We consider the case when E = 0.2 with the biological parameters values given in Table 1 and and the initial conditions (x(0), y(0), z(0)) = (0.2, 0.4, 0.6). Table 1: Parameter values for the proposed Prey-Predator system Parameter Value
0.2 0.3 0.2 0.5 0.5 0.5 0.5
a b c d e α β 6
ACCEPTED MANUSCRIPT Considering, biological parameter values given in Table 1, the equilibrium points, corresponding eigenvalues and the nature of the equilibrium points and the stability of the proposed prey-predator model are presented in Table 2.
RI PT
Table 2: Equilibrium points and their stability of the proposed prey-predator model Equilibrium point Eigenvalues Stable/Unstable (0.85, 0.0, 0.60) 0.2550, −0.3337 + 0.2290i, −0.3337 + 0.2290i Unstable (1, 1, 0) −0.3, −0.2, 0.3 Unstable (0, 0.86, 0.81) 0.1433, −0.4650 + 0.2505i, −0.4650 − 0.2505i Unstable (1, 1, 1.5) −0.3, −1.0017, 0.5017 Unstable (0.56, 0.64, 0.49) −0.2997, −0.0502 + 0.1619i, −0.0502 − 0.1619i Asymptotically stable
EP
TE D
M AN U
SC
From Table 1 it is observed that the different steady state points behave in a different way. All the equilibrium points are unstable except the nontrivial equilibrium point (0.56, 0.64, 0.49). The equilibrium point (0.56, 0.64, 0.49) is locally asymptotically stable as the eigenvalues corresponding to (0.56, 0.64, 0.49) are negative and complex with negative real parts (see Theorem 6). The variations of the three species with respect to time for different α (0.1, 0.5, 0.9) and β (0.1, 0.5, 0.9) are depicted in the Figure 1 starting with initial condition (x(0), y(0), z(0)) = (0.2, 0.4, 0.6).
Figure 1: Variation of the prey and predator species with time for different α and β
AC C
From Figure 1 we observe a natural phenomenon that when refuge population of both the prey groups increases, the stability value of both prey groups also increases. But an interesting case arises when the refuge population of both prey groups increases, then the stability value of predator also increases up to a certain value of α and β (0.1, 0.5). Again, when α and β are equal to 0.9, the predator is going to extinct. Again, the variations of the three species with respect to time for the same values of all the parameters given in Table 1, together with E = 0.1, E = 0.3 and E = 0.5 beginning with initial condition (x(0), y(0), z(0)) = (0.2, 0.4, 0.6) is presented in Figure 2.
7
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Figure 2: Variation of the prey and predator species with time for different E
AC C
EP
TE D
Figure 2. shows that when harvesting effort E of the predator species increases then the stability value of both the prey species also increase. On the contary, when the harvesting effort E of the predator increases then stability value of the predator species gradually decreases. Also, when E = 0.5 then the predator species goes to extinct. The phase space trajectories of the proposed prey–predator model corresponding to the artificially chosen data for different α and β (0.1, 0.5 and 0.9) are presented in Figure 3. It is clear from Figure 3, using different initial population sizes of the species the trajectories for all the species indicate that the nontrivial equilibrium point are spiral for α and β equal to 0.1 and 0.5. But when α and β equal to 0.9 the trajectories for all the species converge at a point. Spiral figure concludes that the stability comes after the oscillation of population size of the three species. It is an interesting fact that when refuge is low then oscillation occurs but when refuge level is high then there does not exist any oscillation.
Figure 3: Phase space trajectories corresponding to different α and β Next, the phase space trajectories of the proposed prey–predator model corresponding to different E (0.1, 0.3 and 0.5) are presented in Figure 4.
8
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 4: Phase portrait for different E and initial points
7
M AN U
Figure 4 shows that using different initial population sizes of the species the trajectories for all the species indicate that the nontrivial equilibrium point are spiral for E equals to 0.1 and 0.3. But when E = 0.5 the trajectories for all the species converge at a point. In this case also spiral figure concludes that the stability comes after the oscillation of population size of three species. It is also an interesting fact that when harvesting E effort is low (0.1 and 0.3) then oscillation occurs but when harvesting effort level is high (E = 0.5) then there does not exist any oscillation.
Conclusion
TE D
From our studies, we see that harvesting and refuge affect the stability of the model. Refuge help prey to increase their population but predator faces a problem of food for this region. When refuge crosses a certain level then predator will be extinct. Similar case arises during harvesting. Here when harvesting increases, predator decreases and prey increases. We also see from the numerical simulation that when refuge and harvesting are low then oscillation in population size occur but when refuge and harvesting are high then there exists no oscillation.
References
EP
Acknowledgments We are grateful to the anonymous referees and the Editor for their careful reading, valuable comments and helpful suggestions which have helped us to improve the presentation of this work significantly.
AC C
[1] A. J. Lotka, Elements of Physical Biology (Williams and Wilkins, Baltimore, 1925). [2] V. Volterra, Variations and fluctuations of a number of individuals in animal species living together, Translation by R. N. Chapman in Animal Ecology (McGraw Hill, NewYork, 1931), pp. 409-448.
[3] W. Ding, H. Finotti, S. Lenhart, Y. Louc, Q. Yed, Optimal control of growth coefficient on a steady-state population model, Nonlinear Anal: Real World Appl. 11 (2010) 688-704.
[4] D. Pal, G. S. Mahapatra, G. P. Samanta, Stability and bionomic analysis of fuzzy prey–predator harvesting model in presence of toxicity: a dynamic approach, Bull Math Biol 78 (2016) 1493–1519.
[5] D. Pal, G. S. Mahapatra, G. P. Samanta, Stability and bionomic analysis of fuzzy parameter based prey predator harvesting model using UFM, Nonlinear Dyn 79 (2014)1939–1955.
[6] D. Pal, G. S. Mahapatra, G. P. Samanta, Optimal harvesting of prey–predator system with interval biological parameters: a bioeconomic model, Math. Biosci. 241 (2013) 181–187.
[7] D. Pal, G. S. Mahapatra, A bioeconomic modeling of two-prey and one-predator fishery model with optimal harvesting policy through hybridization approach, Appl. Math. Comput. 242 (2014) 748–763.
9
ACCEPTED MANUSCRIPT [8] R. P. Gupta, P. Chandra, Bifurcation analysis of modified Leslie–Gower predator–prey model with Michaelis–Menten type prey harvesting, J. Math. Anal. Appl. 398 (2013) 278–295.
[9] L. F. Murphy, S. J. Smith, Optimal harvesting of an age-structured population, J. Math. Biol. 29 (1990) 77–90. [10] A. R. Palma, E. G. Olivares, Optimal harvesting in a predator-prey model with Allee effect and sigmoid functional response, Appl. Math. Model. 5 (2012) 1864–1874.
[11] Z. Luo, X. Fan, Optimal control for an age-dependent competitive species model in a polluted environment, Appl. Math.
RI PT
Comput. 228 (2014) 91–101.
[12] B. Rai, H.L. Freedman and J.F. Addicott, Analysis of three species models of mutualism in predator-prey and competitive systems, Math. Biosci. 65(1983) 13-50.
[13] G. D. Ruxton, Short term refuge use and stability of predator–prey models, Theor. Popul. Biol. 47 (1995) 1–17. [14] H. Huang, F. Chen, L. Zhong, Stability analysis of a prey–predator model with Holling type III response function
SC
incorporating a prey refuge, Appl. Math. Comput. 182 (2006) 672–683.
[15] G. B. Frisvold, Jeanne M. Reeves, The costs and benefits of refuge requirements: The case of Bt cotton, Ecological Economics, 65(1) (2008) 87-97.
[16] Z. Ma, W. Li, Y. Zhao, W. Wang, H. Zhang, Z. Li, Effects of prey refuges on a predator–prey model with a class of
M AN U
function responses: The role of refuges, Math. Biosci. 218 (2009) 73–79.
[17] C. Ross, G. J´ozsef, A predator–prey refuge system: Evolutionary stability in ecological systems, Theor. Popul. Biol. 76(4) (2009) 248–257.
[18] 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.
[19] L. Ji, C. Wu, Qualitative analysis of a predator–prey model with constant-rate prey harvesting incorporating a constant prey refuge, Nonlinear Anal. Real World Appl. 11 (2010) 2285–2295.
TE D
[20] Y. Tao, X. Wang and X. Song, Effects of prey refuge on a harvested predator–prey model with generalized functional response, Commun. Nonlinear Sci. Numer. Simul.16 (2011) 1052–1059.
[21] X. M. Wu, J. W. Li and Z. C. Wang, Existence of positive solutions for a generalized predator–prey model with harvesting term, Comput. Math. Appl. 55 (2008) 1895–1905.
[22] F. Mostafa, H. Mahmoud, Periodic solutions for predator–prey systems with Beddington–DeAngelis functional response
EP
on time scals, Nonlinear Anal. Real World Appl. 9 (2008) 1244–1235.
[23] M. Haque, A detailed study of the Beddington–DeAngelis predator–prey model, Math. Biosci. 234 (2011) 1–16.
AC C
[24] Y. Jia, H. Xu, R. P. Agarwal, Existence of positive solutions for a prey-predator model with refuge and diffusion, Appl. Math. Comput. 217(21) (2011) 8264-8276.
[25] B. Leard, J. Rebaza, Analysis of predator–prey models with continuous threshold harvesting, Appl. Math. Comput. 217 (2011) 5265–5278.
[26] F. D. Chen, Y. M. Wu, Z. Z. Ma, Stability property for the predator-free equilibrium point of predator–prey systems with a class of functional response and prey refuges, Discrete Dynamics in Nature and Society, 2012 (2012) 5 (Article ID 148942).
[27] C. S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Canada 45 (1965) 3-60.
[28] M. F. Elettreby and H. El-Metwally, Multi-team prey-predator model, Int. J. Mod. Phys. C. 18(10) (2007) 1609-1617. [29] T. Park, M. Lloyd, Natural Selection and the Outcome of Competition, The American Naturalist. 89 (1955) 235-240. [30] J. A. Prevedello, C. R. Dickman, M. V. Vieira, E. M. Vieira, Population responses of small mammals to food supply and predators: a global meta-analysis. J. Animal Ecology. 82 (2013) 927-936.
10
ACCEPTED MANUSCRIPT [31] E. Sciubba, F. Zullo, Is sustainability a thermodynamic concept?, Int. J. Exergy. 8 (2011) 68-85. [32] E. Sciubba, F. Zullo, An exergy-based analysis of the co-evolution of different species sharing common resource, Ecol. Model. 273 (2014) 277–283.
[33] Y. M. Svirezhev, Nonlinearities in mathematical ecology: Phenomena and models. Would we live in Volterra’s world?, Ecol. Model. 216 (2008) 89–101.
AC C
EP
TE D
M AN U
SC
RI PT
[34] J. K. Hale, Theory of functional differential equations. Springer, New York, 1977.
11
S2405-9854(17)30002-2
DOI:
10.1016/j.egg.2017.05.001
Reference:
EGG 9
To appear in:
Ecological Genetics and Genomics
Received Date: 26 January 2017 Revised Date:
14 April 2017
Accepted Date: 15 May 2017
Please cite this article as: D. Pal, P. Santra, G.S. Mahapatra, Dynamical behavior of three species predator–prey system with mutual support between non refuge prey, Ecological Genetics and Genomics (2017), doi: 10.1016/j.egg.2017.05.001. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
Dynamical Behavior of Three Species Predator–Prey System with Mutual Support Between Non Refuge Prey D. Pal1∗ , P. Santra2 and G. S. Mahapatra3 2
Chandrahati Dilip Kumar High School(H.S.), Chandrahati 712504, West Bengal, India.
RI PT
1
Department of Mathematics, Abada New Set Up Upper Primary School, Howrah, West Bengal, India 3
Department of Mathematics, National Institute of Technology Puducherry, Karaikal-609609, India. Tel. +919433135327, Fax: 04368231665
E-mail: [email protected]
SC
∗ Correspondence
1
M AN U
Abstract: In this paper, we study the effect of refuge and harvesting in a three species model of two prey and one predator, where the two prey refuge groups are helping each other and shielding from the predator group. The stability theory of nonlinear differential equations is used to analyze the model. Our study exhibits that, both the refuge and harvesting play a significant role on the stability of the system. Along with the local stability of the system, the global stability is also studied. Finally, numerical simulations are given to support the analytical conclusion. Keywords: Prey-Predator; Harvesting; Refuge; Stability; Mutualism.
Introduction
AC C
EP
TE D
In population dynamics, the relationship between prey-predator plays an importance for universal existence. Lotka [1] and Volterra [2] first studied the relationship between prey and predator. During the last few decades, several research articles have been published on the prey-predator system incorporating refuge concept and mutualism. The study of stability and harvesting of species including refuge is growing in importance. Due to the mutual support between the species, the study of harvesting of the prey-predator system is enhanced based on the study of refuge and mutualness of prey-predator system, as per the diverse nature of our ecological system. Recently, several researchers (Ding et. al. [3], Pal et al. ([4]-[6]), Pal and Mahapatra [7], Gupta and Chandra [8], Murphy and Smith [9], Palma and Olivares [10], Luo et al. [11]) concentrated their research work in the field of harvesting management of renewable resources. We have studied the refuge and harvesting of two prey and one predator model. Rai et al. [12] analyzed the model of mutualism in predator-prey and competitive systems in three species environment. Ruxton [13] studied the stability of prey-predator model using short term refuge. Huang et. al. [14] presented a prey–predator model with Holling type III response function incorporating a prey refuge. Frisvold and Reeves [15] studied the costs and benefits of refuge requirements. Ma et al. [16] discussed the effects of prey refuges on a prey-predator model. Ross and J´ozsef [17] studied a predator–prey refuge system. Chen et al. [18] studied qualitative analysis of a predator-prey model with Holling type II functional response incorporating a constant prey refuge. Ji and Wu [19] presented qualitative analysis of a prey-predator model with constant rate of prey harvesting incorporating a constant prey refuge. Tao et al.[20] discussed the effects of prey refuge on a harvested predator–prey model with generalized functional response. Wu et al. [21] derived a sufficient condition for the existence of at least one positive periodic solution of a generalized prey–predator model with harvesting term. Mostafa and Mahmoud [22] presented periodic solutions of a prey-predator system with Beddington–DeAngelis functional response [23] on time scales. Jia et al. [24] studied the existence of positive solutions for a prey-predator model with refuge and diffusion. Leard and Rebaza [25] presented a prey-predator model with continuous threshold harvesting and formulated an alternative form of threshold harvesting. Chen et al. [26] studied the stability property for the predator-free equilibrium point of prey-predator systems with a class of functional response [27] and prey refuges. In this paper, we consider a model which consists of two prey groups and one predator group. Prey groups help each other against the predator and few prey of each prey groups are the refuge. Here we consider an interesting situation, i.e., when a predator attacks the 1st prey group then some prey of the 1st prey group will be refuge, but the total population of 2nd prey group help 1st prey group against the predator and when a predator attacks the
1
ACCEPTED MANUSCRIPT
2
RI PT
2nd prey group then some prey of the 2nd prey group will be refuge but the total population of 1st prey groups help the 2nd prey groups against the predator. In this model, we also consider that the predator species is harvested. We have constructed the model, motivated by the work of Elettreby and El-Metwally [28] in which they have studied the local and global stability of multi-team prey-predator model. The rest of the paper is organized as follows: In section 2, we present the mathematical formulation of three species prey-predator system. Positivity of the model is presented in Section 3. Section 4 present the possible steady states of the proposed predator-pray model. In section 5, we study the local stability as well as global stability analysis of the equilibrium points. The Numerical simulation is presented in section 6. Finally, Section 7 gives concluding remark of our work.
Mathematical model of proposed predator–prey system
M AN U
SC
In our proposed model, two teams of prey with densities x(t) and y(t) respectively, interact with one team of predator with densities z(t). The assumptions of the prey-predator model are as follows: (i) Each team of preys grows logistically in absence of any predator, (ii) Due to effect of predation, the prey growth rate is reduced by a quantity proportional to prey and predator population, (iii) The non refuge teams of prey help each other against the predator, (iv) In absence of any prey for sustenance, the predator’s death rate results is in inverse decay, i.e., −cz 2 , (v) Prey’s contribution to predator growth rate is proportional to available prey and size of predator population, (vi) Predator species is harvested with the harvesting effort E, (vii) Some preys of each team are refuge, i.e., αx and βy, (viii) When the team of 2nd prey helps first prey against predator then not only the total populations of 2nd team involved in this interaction but also non refuge population of 1st team involve in this interaction. Same situation occurs for the first prey support to the second prey. The rate of primary resources available for the prey to feed plays an important role in this kind of problems ([29]-[33]). Using the above assumptions, the following predator-prey model is proposed:
TE D
dx = ax(1 − x) − (1 − α)xz + (1 − α)xyz dt dy = by(1 − y) − (1 − β)yz + (1 − β)xyz dt dz = −cz 2 + d(1 − α)xz + e(1 − β)yz − Ez dt
(1)
EP
where the coefficients a, b, c, d and e are positive constants. With initial densities and x(0) > 0, y(0) > 0 and z(0) > 0
(2)
3
AC C
It is clear that the two teams of preys help each other, e.g., in foraging and in early warning against predation. Note that the mutual support occurs only in the presence of the predator.
Positivity of the Proposed System
Theorem 1 Every solution of the system (1) with initial conditions (2) exists in the interval [0, ∞) and x(0) > 0, y(0) > 0 and z(0) > 0 for all t ≥ 0. Proof. Since the right-hand side of system (1) is completely continuous and locally Lipschitzian on C, the solution (x(t), y(t), z(t)) of (1) with initial conditions (2) exist and is unique on [0, ξ), where 0 < ξ < ∞ [34]. From system
2
ACCEPTED MANUSCRIPT (1) with initial conditions (2), we have t Z x (t) = x (0) exp {a(1 − x (θ)) − (1 − α)z (θ) + (1 − α)y (θ) z (θ)} dθ > 0 0t Z y (t) = y (0) exp {b(1 − y (θ)) − (1 − β)z (θ) + (1 − β)x (θ) z (θ)} dθ > 0
RI PT
0t Z z (t) = z (0) exp {−cz (θ) + d(1 − α)x (θ) + e(1 − β)y (θ) − E} dθ > 0 0
which completes the proof.
4
Steady States Analysis of the Proposed Predator-Prey System
M AN U
SC
The possible steady states of the predator-prey dynamical system (1) are as follows: G0 (0, 0, − Ec ) (this equilibrium point does not biologically relevant as − Ec < 0), G1 (0, 0, 0), G2 (1, 0, 0), G3 (0, 1, 0), G4 (1, 1, 0) and the other equilibrium points are be(1−β)−Eb cb+E(1−β) G5 (0, y5 , z5 ): where y5 = cb+e(1−β) 2 > 0 and z5 = cb+e(1−β)2 . The equilibrium point G5 exists if e(1 − β) > E. ac+E(1−α) ad(1−α)−aE ac+d(1−α)2 > 0 and z6 = ac+d(1−α)2 . The equilibrium point G6 exists if d(1 − α) > E. . The equilibrium point G6 exists if d(1 − α) + e(1 − β) > E. G7 (1, 1, z7 ): where z7 = d(1−α)+e(1−β)−E c q q 2 ab ab bc+e(1−β) ac+d(1−α)2 }+a{E−d(1−α)} }+b{E−e(1−β)} (1−α)(1−β) { (1−α)(1−β) { q q G8 (x8 , y8 , z8 ): where x8 = , y = and 8 ab ab 2 2 e(1−β) d(1−α) +bd(1−α) +ae(1−β) (1−α)(1−β) (1−α)(1−β)
G6 (x6 , 0, z6 ): where x6 =
z8 =
ab (1−α)(1−β)
> 0. The equilibrium point G8 exists if
s ab ab 2 bc + e(1 − β) +b {E − e(1 − β)} > 0 and ac + d(1 − α)2 +a {E − d(1 − α)} > 0 (1 − α)(1 − β) (1 − α)(1 − β) −
q
TE D
s
q
ab (1−α)(1−β)
{bc+e(1−β)2 }+b{E−e(1−β)}
−
q
ab (1−α)(1−β)
q G9 (x9 , y9 , z9 ): where x9 = , y9 = ab −e(1−β)2 (1−α)(1−β) +bd(1−α) q ab < 0.The equilibrium point G9 does not exist since z9 < 0. z9 = − (1−α)(1−β)
−d(1−α)2
q
ab +ae(1−β) (1−α)(1−β)
and
Stability Analysis
EP
5
{ac+d(1−α)2 }+a{E−d(1−α)}
5.1
AC C
In this section, we present local as well as global stability of the equilibrium points of the proposed predator-prey System.
Local Stability analysis
Here, we study the local stability analysis of our proposed model. It is clear that the equilibrium points G1 (0, 0, 0), G2 (1, 0, 0), G3 (0, 1, 0) are unstable. Now we analyze the local stability at the other equilibrium points of the system based on the following Theorems. The variational matrix M of the system (1) corresponding to different equilibrium points is given by M11 (1 − α)xz (1 − α) (y − 1) x M22 (1 − β) (x − 1) y M = (1 − β)yz d(1 − α)z e(1 − β)z M33
Here,
M11 = a − 2ax − (1 − α)z + (1 − α)yz, M22 = b − 2by − (1 − β)z + (1 − β)xz, M33 = −2cz + d(1 − α)x + e(1 − β)y − E,
3
ACCEPTED MANUSCRIPT Theorem 2 The equilibrium point G4 (1, 1, 0) is stable if d(1 − α) + e(1 − β) < E. Proof. The variational matrix M corresponding to the equilibrium point G4 (1, 1, 0) has eigenvalues −a, −b and d(1 − α) + e(1 − β) − E respectively. Since −a < 0 and −b < 0, the equilibrium point G4 (1, 1, 0) is stable if d(1 − α) + e(1 − β) − E < 0 i.e., d(1 − α) + e(1 − β) < E. a Theorem 3 The equilibrium point G5 (0, y5 , z5 ) is locally asymptotically stable if (y5 − 1)z5 < − (1−α) .
SC
The characteristic equation of the above matrix MG5 (0,y5 ,z5 ) is given by det MG5 (0,y5 ,z5 ) − λI = 0
RI PT
Proof. The variational matrix at the equilibrium point G5 (0, y5 , z5 ) is given by a + (1 − α)(y5 − 1)z5 0 0 (1 − β)y5 z5 −by5 −(1 − β)y5 MG5 (0,y5 ,z5 ) = d(1 − α)z5 e(1 − β)z5 −cz5
(3)
One root of the characteristic equation (3) i.e. one eigenvalue is given by a + (1 − α)(y5 − 1)z5 . This eigenvalue is a . The other two eigenvalues are the root of the quadratic equation negative if (y5 − 1)z5 < − (1−α) (4)
M AN U
λ2 + λ(cz5 + by5 ) + y5 z5 {cb + e(1 − β)2 } = 0
Now, sum of the roots of the equation (4) is −(cz5 + by5 ), which is negative. The product of the roots of the equation (4) is y5 z5 {cb + e(1 − β)2 }, which is always positive. Hence, the roots of the quadratic equation (4) are real and negative or complex conjugates with negative real parts. Therefore, G5 (0, y5 , z5 ) is locally asymptotically stable if a . (y5 − 1)z5 < − (1−α) b . Theorem 4 The equilibrium point G6 (x6 , 0, z6 ) is locally asymptotically stable if (x6 − 1)z6 < − (1−β)
TE D
Proof. The variational matrix at the equilibrium point G6 (x6 , 0, z6 ) is given by −ax6 (1 − α)x6 z6 −(1 − α)x6 MG6 (x6 ,0,z6 ) = 0 b + (1 − β)(x6 − 1)z6 0 d(1 − α)z6 e(1 − β)z6 −cz6
EP
The characteristic equation of the above matrix MG6 (x6 ,0,z6 ) is given by det MG6 (x6 ,0,z6 ) − λI = 0
(5)
One root of the characteristic equation (5) is given by b + (1 − β)(x6 − 1)z6 . This eigenvalue is negative if (x6 − 1)z6 < b .The other two eigenvalues are the root of the quadratic equation − (1−β)
AC C
λ2 + λ(cz6 + ax6 ) + x6 z6 {ca + d(1 − α)2 } = 0
(6)
The sum of the roots of the equation (6) is−(cz6 + ax6 ) < 0. The product of the roots of the equation (6) is x6 z6 {ca + d(1 − α)2 }, which is always positive. Hence, the roots of the above quadratic equation (6) are real and negative or complex conjugates with negative real parts. Therefore, G6 (x6 , 0, z6 ) is locally asymptotically stable if b (x6 − 1)z6 < − (1−β) . Theorem 5 The equilibrium point G7 (1, 1, z7 ) is locally asymptotically stable if ab > (1 − α)(1 − β)z72 . Proof. The variational matrix at the equilibrium point G7 (1, 1, z7) is given by −a (1 − α)z7 0 −b 0 MG7 (1,1,z7 ) = (1 − β)z7 d(1 − α)z7 e(1 − β)z7 −cz7
The characteristic equation of the matrix MG6 (x6 ,0,z6 ) is given by
det MG7 (1,1,z7 ) − λI = 0 4
(7)
ACCEPTED MANUSCRIPT One root of the characteristic equation (7) is −cz7 , which is negative. The other two eigenvalues are the root of the quadratic equation λ2 + λ(a + b) + {ab − (1 − α)(1 − β)z72 } = 0 (8) The sum of the roots of the equation (8) is −(a + b) < 0. The product of the roots of the equation (8) is{ab − (1 − α)(1 − β)z72 }. The product of the roots is positive if ab > (1 − α)(1 − β)z72 . Therefore, if ab > (1 − α)(1 − β)z72 then the roots of the above quadratic equation (8) are real and negative or complex conjugates with negative real parts. Hence, G7 (1, 1, z7 ) is locally asymptotically stable if ab > (1 − α)(1 − β)z72 .
RI PT
Theorem 6 The equilibrium point G8 (x8 , y8 , z8 ) is locally asymptotically stable if {ac + d(1 − α)2 }x8 + {bc + e(1 − β)2 }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2(y8 − 1)} + (1 − α)(1 − β)z8{d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0.
SC
Proof. The variational matrix at the equilibrium point G8 (x8 , y8 , z8 ) is given by −ax8 (1 − α)x8 z8 (1 − α)x8 (y8 − 1) −by8 (1 − β)y8 (x8 − 1) MG8 (x8 ,y8 ,z8 ) = (1 − β)y8 z8 d(1 − α)z8 e(1 − β)z8 −cz8 The characteristic equation of the the above matrix MG8 (x8 ,y8 ,z8 ) is given by
(9)
M AN U
f (λ) = 0
Global stability analysis of predator–prey system
EP
5.2
TE D
where f (λ) = λ3 + λ2 ω1 + λω2 + ω3 , ω1 = (ax8 + by8 + cz8 ), ω2 = [{ac + d(1 − α)2 }x8 z8 + {bc + e(1 − β)2 }y8 z8 − {e(1 − β)2 + d(1 − α)2 }x8 y8 z8 ], ω3 = −[x8 y8 z8 {ae(1 − β)2 (x8 − 1) + bd(1 − α)2 (y8 − 1)} + (1 − α)(1 − β)x8 y8 z82 {d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}]. Now ω1 > 0, ω2 > 0 if {ac + d(1 − α)2 }x8 + {bc + e(1 − β)2 }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and ω3 > 0 if [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2 (y8 − 1)} + (1 − α)(1 − β)z8 {d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0. The roots λi (i = 1, 2, 3), of the polynomial equation (9)satisfies the following conditions (i) λ1 + λ2 + λ3 = −ω1 < 0 (ii)λ1 λ2 + λ2 λ3 + λ3 λ1 = ω2 > 0 (iii)λ1 λ2 λ3 = −ω3 < 0. Since ω1 , ω2 and ω3 are all positive, so there is no change of sign in the polynomial f (λ) for λ > 0. Therefore, f (λ) = 0 does not possess any positive root. Again in f (−λ), there are three changes of sign. Therefore, the nature of the roots f (λ) = 0 has three possibilities such as (i) no positive root (ii) one negative root and two complex roots with negative real parts or (iii) three negative roots. Therefore, the equilibrium point G8 (x8 , y8 , z8 ) is locally asymptotically stable if {ac + d(1 − α)2 }x8 + {bc + e(1 − 2 β) }y8 > {e(1 − β)2 + d(1 − α)2 }x8 y8 and [{ae(1 − β)2 (x8 − 1) + bd(1 − α)2(y8 − 1)} + (1 − α)(1 − β)z8{d(1 − α)(x8 − 1) + e(1 − β)(y8 − 1)}] < 0.
AC C
Here, we study the global stability of the interior equilibrium points G7 (1, 1, z7 ) and G8 (x8 , y8 , z8 ) by considering suitable Lyapunov function. Theorem 7 Equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable. Proof. In order to test whether the equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable or not, we construct a suitable Lyapunov function as follows z V = [(x − 1) − log x] + l1 [(y − 1) − log y] + l2 (z − z7 ) − z7 log (10) z7 It is obvious that V is positive definite. Differentiating both sides of (10) with respect to time t along the solutions of (1) we get, dV dt
⇒
= dV dt
dy z−z7 dz + l1 y−1 y dt + l2 z dt 2 = −a(x − 1) + (1 − α)(x − 1)(yz − z7 ) − l1 b(y − 1)2 + l1 (1 − β)(y − 1)(xz − z7 ) − l2 c(z − z7 )2
x−1 dx x dt
The above equation can be written in real quadratic form dV = −X T P X dt 5
ACCEPTED MANUSCRIPT where X T = [(x − 1) , (y − 1) , (z − z7 )] and 7) a −(1 − α) (yz−z x−1 P = 0 bl1 0 0
0 7) −l1 (1 − β) (xz−z y−1 cl2
Theorem 8 Equilibrium point G8 (x8 , y8 , z8 ) is globally asymptotically stable.
RI PT
Now dV dt < 0, if the matrix P is positive definite. The matrix P is positive definite as all the principal minors are positive and hence dV dt < 0. Therefore, the interior equilibrium point G7 (1, 1, z7 ) is globally asymptotically stable.
Proof. we construct a suitable Lyapunov function as follows: x y z W = (x − x8 ) − x8 log + l1 (y − y8 ) − y8 log + l2 (z − z8 ) − z8 log x8 y8 z8
(11)
dW dt
dW dt
z−z8 dz 8 dy + l1 y−y y dt + l2 z dt = −a(x − x8 )2 + (1 − α)(x − x8 )(yz − y8 z8 ) − l1 b(y − y8 )2 + l1 (1 − β)(y − y8 )(xz − x8 z8 ) − l2 c(z − z8 )2
x−x8 dx x dt
M AN U
⇒
=
SC
It is obvious that W is positive definite. Differentiating both sides of (11) with respect to time t along the solutions of (1) we get,
The above equation can be written in real quadratic form
dW = − − Y T QY dt
TE D
where Y T = [(x − x8 ) , (y − y8 ) , (z − z8 )] and 8 z8 ) a −(1 − α) (yz−y x−x8 Q= 0 bl1 0 0
0 8 z8 ) −l1 (1 − β) (xz−x y−y8 cl2
Now dV dt < 0, if the matrix Q is positive definite. The matrix Q is positive definite as all the principal minors are positive and hence dW dt < 0. Therefore, the equilibrium point G8 (x8 , y8 , z8 ) is globally asymptotically stable.
Numerical Simulation
EP
6
AC C
Analytical studies can never be completed without numerical verification of the derived results. In this section we present computer simulation of some solutions of the systems (1). Beside verification of our analytical findings, these numerical solutions are very important from practical point of view. For the purpose of simulation experiments we mainly use the software MATLAB R2008a and Wolfram Mathematica 8. As the problem is not a case study for a particular species, here some hypothetical data are taken for the sole purpose of illustrating the analytical results of the previous sections. Moreover, it may be noted that as the parameters of the model are not based on real world observations, the main features described by the simulations presented in this section should be considered from a qualitative, rather than a quantitative point of view. We consider the case when E = 0.2 with the biological parameters values given in Table 1 and and the initial conditions (x(0), y(0), z(0)) = (0.2, 0.4, 0.6). Table 1: Parameter values for the proposed Prey-Predator system Parameter Value
0.2 0.3 0.2 0.5 0.5 0.5 0.5
a b c d e α β 6
ACCEPTED MANUSCRIPT Considering, biological parameter values given in Table 1, the equilibrium points, corresponding eigenvalues and the nature of the equilibrium points and the stability of the proposed prey-predator model are presented in Table 2.
RI PT
Table 2: Equilibrium points and their stability of the proposed prey-predator model Equilibrium point Eigenvalues Stable/Unstable (0.85, 0.0, 0.60) 0.2550, −0.3337 + 0.2290i, −0.3337 + 0.2290i Unstable (1, 1, 0) −0.3, −0.2, 0.3 Unstable (0, 0.86, 0.81) 0.1433, −0.4650 + 0.2505i, −0.4650 − 0.2505i Unstable (1, 1, 1.5) −0.3, −1.0017, 0.5017 Unstable (0.56, 0.64, 0.49) −0.2997, −0.0502 + 0.1619i, −0.0502 − 0.1619i Asymptotically stable
EP
TE D
M AN U
SC
From Table 1 it is observed that the different steady state points behave in a different way. All the equilibrium points are unstable except the nontrivial equilibrium point (0.56, 0.64, 0.49). The equilibrium point (0.56, 0.64, 0.49) is locally asymptotically stable as the eigenvalues corresponding to (0.56, 0.64, 0.49) are negative and complex with negative real parts (see Theorem 6). The variations of the three species with respect to time for different α (0.1, 0.5, 0.9) and β (0.1, 0.5, 0.9) are depicted in the Figure 1 starting with initial condition (x(0), y(0), z(0)) = (0.2, 0.4, 0.6).
Figure 1: Variation of the prey and predator species with time for different α and β
AC C
From Figure 1 we observe a natural phenomenon that when refuge population of both the prey groups increases, the stability value of both prey groups also increases. But an interesting case arises when the refuge population of both prey groups increases, then the stability value of predator also increases up to a certain value of α and β (0.1, 0.5). Again, when α and β are equal to 0.9, the predator is going to extinct. Again, the variations of the three species with respect to time for the same values of all the parameters given in Table 1, together with E = 0.1, E = 0.3 and E = 0.5 beginning with initial condition (x(0), y(0), z(0)) = (0.2, 0.4, 0.6) is presented in Figure 2.
7
SC
RI PT
ACCEPTED MANUSCRIPT
M AN U
Figure 2: Variation of the prey and predator species with time for different E
AC C
EP
TE D
Figure 2. shows that when harvesting effort E of the predator species increases then the stability value of both the prey species also increase. On the contary, when the harvesting effort E of the predator increases then stability value of the predator species gradually decreases. Also, when E = 0.5 then the predator species goes to extinct. The phase space trajectories of the proposed prey–predator model corresponding to the artificially chosen data for different α and β (0.1, 0.5 and 0.9) are presented in Figure 3. It is clear from Figure 3, using different initial population sizes of the species the trajectories for all the species indicate that the nontrivial equilibrium point are spiral for α and β equal to 0.1 and 0.5. But when α and β equal to 0.9 the trajectories for all the species converge at a point. Spiral figure concludes that the stability comes after the oscillation of population size of the three species. It is an interesting fact that when refuge is low then oscillation occurs but when refuge level is high then there does not exist any oscillation.
Figure 3: Phase space trajectories corresponding to different α and β Next, the phase space trajectories of the proposed prey–predator model corresponding to different E (0.1, 0.3 and 0.5) are presented in Figure 4.
8
SC
RI PT
ACCEPTED MANUSCRIPT
Figure 4: Phase portrait for different E and initial points
7
M AN U
Figure 4 shows that using different initial population sizes of the species the trajectories for all the species indicate that the nontrivial equilibrium point are spiral for E equals to 0.1 and 0.3. But when E = 0.5 the trajectories for all the species converge at a point. In this case also spiral figure concludes that the stability comes after the oscillation of population size of three species. It is also an interesting fact that when harvesting E effort is low (0.1 and 0.3) then oscillation occurs but when harvesting effort level is high (E = 0.5) then there does not exist any oscillation.
Conclusion
TE D
From our studies, we see that harvesting and refuge affect the stability of the model. Refuge help prey to increase their population but predator faces a problem of food for this region. When refuge crosses a certain level then predator will be extinct. Similar case arises during harvesting. Here when harvesting increases, predator decreases and prey increases. We also see from the numerical simulation that when refuge and harvesting are low then oscillation in population size occur but when refuge and harvesting are high then there exists no oscillation.
References
EP
Acknowledgments We are grateful to the anonymous referees and the Editor for their careful reading, valuable comments and helpful suggestions which have helped us to improve the presentation of this work significantly.
AC C
[1] A. J. Lotka, Elements of Physical Biology (Williams and Wilkins, Baltimore, 1925). [2] V. Volterra, Variations and fluctuations of a number of individuals in animal species living together, Translation by R. N. Chapman in Animal Ecology (McGraw Hill, NewYork, 1931), pp. 409-448.
[3] W. Ding, H. Finotti, S. Lenhart, Y. Louc, Q. Yed, Optimal control of growth coefficient on a steady-state population model, Nonlinear Anal: Real World Appl. 11 (2010) 688-704.
[4] D. Pal, G. S. Mahapatra, G. P. Samanta, Stability and bionomic analysis of fuzzy prey–predator harvesting model in presence of toxicity: a dynamic approach, Bull Math Biol 78 (2016) 1493–1519.
[5] D. Pal, G. S. Mahapatra, G. P. Samanta, Stability and bionomic analysis of fuzzy parameter based prey predator harvesting model using UFM, Nonlinear Dyn 79 (2014)1939–1955.
[6] D. Pal, G. S. Mahapatra, G. P. Samanta, Optimal harvesting of prey–predator system with interval biological parameters: a bioeconomic model, Math. Biosci. 241 (2013) 181–187.
[7] D. Pal, G. S. Mahapatra, A bioeconomic modeling of two-prey and one-predator fishery model with optimal harvesting policy through hybridization approach, Appl. Math. Comput. 242 (2014) 748–763.
9
ACCEPTED MANUSCRIPT [8] R. P. Gupta, P. Chandra, Bifurcation analysis of modified Leslie–Gower predator–prey model with Michaelis–Menten type prey harvesting, J. Math. Anal. Appl. 398 (2013) 278–295.
[9] L. F. Murphy, S. J. Smith, Optimal harvesting of an age-structured population, J. Math. Biol. 29 (1990) 77–90. [10] A. R. Palma, E. G. Olivares, Optimal harvesting in a predator-prey model with Allee effect and sigmoid functional response, Appl. Math. Model. 5 (2012) 1864–1874.
[11] Z. Luo, X. Fan, Optimal control for an age-dependent competitive species model in a polluted environment, Appl. Math.
RI PT
Comput. 228 (2014) 91–101.
[12] B. Rai, H.L. Freedman and J.F. Addicott, Analysis of three species models of mutualism in predator-prey and competitive systems, Math. Biosci. 65(1983) 13-50.
[13] G. D. Ruxton, Short term refuge use and stability of predator–prey models, Theor. Popul. Biol. 47 (1995) 1–17. [14] H. Huang, F. Chen, L. Zhong, Stability analysis of a prey–predator model with Holling type III response function
SC
incorporating a prey refuge, Appl. Math. Comput. 182 (2006) 672–683.
[15] G. B. Frisvold, Jeanne M. Reeves, The costs and benefits of refuge requirements: The case of Bt cotton, Ecological Economics, 65(1) (2008) 87-97.
[16] Z. Ma, W. Li, Y. Zhao, W. Wang, H. Zhang, Z. Li, Effects of prey refuges on a predator–prey model with a class of
M AN U
function responses: The role of refuges, Math. Biosci. 218 (2009) 73–79.
[17] C. Ross, G. J´ozsef, A predator–prey refuge system: Evolutionary stability in ecological systems, Theor. Popul. Biol. 76(4) (2009) 248–257.
[18] 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.
[19] L. Ji, C. Wu, Qualitative analysis of a predator–prey model with constant-rate prey harvesting incorporating a constant prey refuge, Nonlinear Anal. Real World Appl. 11 (2010) 2285–2295.
TE D
[20] Y. Tao, X. Wang and X. Song, Effects of prey refuge on a harvested predator–prey model with generalized functional response, Commun. Nonlinear Sci. Numer. Simul.16 (2011) 1052–1059.
[21] X. M. Wu, J. W. Li and Z. C. Wang, Existence of positive solutions for a generalized predator–prey model with harvesting term, Comput. Math. Appl. 55 (2008) 1895–1905.
[22] F. Mostafa, H. Mahmoud, Periodic solutions for predator–prey systems with Beddington–DeAngelis functional response
EP
on time scals, Nonlinear Anal. Real World Appl. 9 (2008) 1244–1235.
[23] M. Haque, A detailed study of the Beddington–DeAngelis predator–prey model, Math. Biosci. 234 (2011) 1–16.
AC C
[24] Y. Jia, H. Xu, R. P. Agarwal, Existence of positive solutions for a prey-predator model with refuge and diffusion, Appl. Math. Comput. 217(21) (2011) 8264-8276.
[25] B. Leard, J. Rebaza, Analysis of predator–prey models with continuous threshold harvesting, Appl. Math. Comput. 217 (2011) 5265–5278.
[26] F. D. Chen, Y. M. Wu, Z. Z. Ma, Stability property for the predator-free equilibrium point of predator–prey systems with a class of functional response and prey refuges, Discrete Dynamics in Nature and Society, 2012 (2012) 5 (Article ID 148942).
[27] C. S. Holling, The functional response of predators to prey density and its role in mimicry and population regulation, Mem. Entomol. Soc. Canada 45 (1965) 3-60.
[28] M. F. Elettreby and H. El-Metwally, Multi-team prey-predator model, Int. J. Mod. Phys. C. 18(10) (2007) 1609-1617. [29] T. Park, M. Lloyd, Natural Selection and the Outcome of Competition, The American Naturalist. 89 (1955) 235-240. [30] J. A. Prevedello, C. R. Dickman, M. V. Vieira, E. M. Vieira, Population responses of small mammals to food supply and predators: a global meta-analysis. J. Animal Ecology. 82 (2013) 927-936.
10
ACCEPTED MANUSCRIPT [31] E. Sciubba, F. Zullo, Is sustainability a thermodynamic concept?, Int. J. Exergy. 8 (2011) 68-85. [32] E. Sciubba, F. Zullo, An exergy-based analysis of the co-evolution of different species sharing common resource, Ecol. Model. 273 (2014) 277–283.
[33] Y. M. Svirezhev, Nonlinearities in mathematical ecology: Phenomena and models. Would we live in Volterra’s world?, Ecol. Model. 216 (2008) 89–101.
AC C
EP
TE D
M AN U
SC
RI PT
[34] J. K. Hale, Theory of functional differential equations. Springer, New York, 1977.
11