# The Dynamical Analysis of a Prey-Predator Model with a Refuge-Stage Structure Prey Population

Authors: Raid Kamel Naji; Salam Jasim Majeed
Journal: International Journal of Differential Equations
 (2016)

---

## Abstract

We proposed and analyzed a mathematical model dealing with two species of prey-predator system. It is assumed that the prey is a stage structure population consisting of two compartments known as immature prey and mature prey. It has a refuge capability as a defensive property against the predation. The existence, uniqueness, and boundedness of the solution of the proposed model are discussed. All the feasible equilibrium points are determined. The local and global stability analysis of them are investigated. The occurrence of local bifurcation (such as saddle node, transcritical, and pitchfork) near each of the equilibrium points is studied. Finally, numerical simulations are given to support the analytic results.


---

## Body

## 1. Introduction

The development of the qualitative analysis of ordinary differential equations is deriving to study many problems in mathematical biology. The modeling for the population dynamics of a prey-predator system is one of the important and interesting goals in mathematical biology, which has received wide attention by several authors [1–6]. In the natural world many kinds of prey and predator species have a life history that is composed of at least two stages: immature and mature, and each stage has different behavioral properties. So, some works of stage structure prey-predator models have been provided in a good number of papers in the literatures [7–12]. Zhang et al. in [9] and Cui and Takeuchi in [11] proposed two mathematical models of prey stage structure; in these models the predator species consumes exclusively the immature prey. Indeed, there are many factors that impact the dynamics of prey-predator interactions such as disease, harvesting, prey refuge, delay, and many other factors. Several prey species have gone to extinction, and this extinction must be caused by external effects such as overutilization, overpredation, and environmental factors (pollution, famine). Prey may avoid becoming attacked by predators either by protecting themselves or by living in a refuge where it will be out of sight of predators. Some theoretical and empirical studies have shown and tested the effects of prey refuges and making an opinion that the refuges used by prey have a stabilizing influence on prey-predator interactions; also prey extinction can be prevented by the addition of refuges; for instance, we can refer to [13–23]. Therefore, it is important and worthwhile to study the effects of a refuge on the prey population with prey stage structure. Consequently, in this paper, we proposed and analyzed the prey-predator model involving a stage structure in prey population together with a preys refuge property as a defensive property against the predation.


## 2. Mathematical Model

In this section a prey-predator model with a refuge-stage structure prey population is proposed for study. Letx(T) represent the population size of the immature prey at time T; y(T) represents the population size of the mature prey at time T, while z(T) denotes the population size of the predator species at time T. Therefore in order to describe the dynamics of this model mathematically the following hypotheses are adopted:(1)
The immature prey grows exponentially depending completely on its parents with growth rater>0. There is an intraspecific competition between their individuals with intraspecific competition rate δ1>0. The immature prey individual becomes mature with grownup rate β>0 and faces natural death with a rate d1>0.(2)
There is an intraspecific competition between the individuals of mature prey population with intraspecific competition rateδ2>0. Further the mature prey species faces natural death rate too with a rate d2>0.(3)
The environment provides partial protection of prey species against the predation with a refuge rate0<m<1; therefore there is 1-m of prey species available for predation.(4)
There is an intraspecific competition between the individuals of predator population with intraspecific competition rateδ3>0. Further the predator species faces natural death rate too with a rate d3>0.The predator consumes the prey in both the compartments according to the mass action law represented by Lotka–Volterra type of functional response with conversion rates 0<e1<1 and 0<e2<1 for immature and mature prey, respectively. Consequently, the dynamics of this model can be represented mathematically with the following set of differential equations.(1)x˙=ry-δ1x2-d1x-βx-γ11-mxzy˙=βx-δ2y2-d2y-γ21-myzz˙=e1γ11-mxz+e2γ21-myz-δ3z2-d3zwith initial conditions, x(0)>0, y(0)>0, and z(0)>0.In order to simplify the analysis of system (1) the number of parameters in system (1) is reduced using the following dimensionless variables in system (1):(2)y1=e1γ11-md2x,y2=e2γ21-md2y,y3=δ3d2z,t=d2T.Accordingly, the dimensionless form of system (1) can be written as(3)y˙1=a1y2-a2y12-a3y1-a4y1y3y˙2=b1y1-b2y22-y2-b3y2y3y˙3=y1y3+y2y3-y32-b4y3,where the dimensionless parameters are given by(4)a1=e1γ1re2γ2d2,a2=δ1e1γ11-m,a3=d1+βd2,a4=γ11-mδ3,b1=e2γ2βe1γ1d2,b2=δ2e2γ21-m,b3=γ21-mδ3,b4=d3d2.According to the equations given in system (3), all the interaction functions are continuous and have a continuous partial derivatives. Therefore they are Lipschitzain and hence the solution of system (3) exists and is unique. Moreover the solution of system (3) is bounded as shown in the following theorem.Theorem 1.
All the solutions of system (3) that initiate in the positive octant are uniformly bounded.Proof.
LetW=y1+y2+y3 be solutions of system (3) with initial conditions, y1(0)>0, y2(0)>0, and y3(0)>0. Then by differentiation W with respect to t we get(5)dWdt≤a1y21-y2a1/b2+b1y11-y1b1/a2-σ1W≤σ2-σ1W;here σ1=min1,a3,b4 and σ2=a12/4b2+b12/4a2. Now by using comparison theorem, we get(6)0<W≤W0e-σ1t+σ2σ11-e-σ1t.Thus for t→∞ we obtain 0<W≤σ2/σ1. Hence, all solutions of system (3) in R+3 are uniformly bounded and therefore we have finished the proof.


## 3. Local Stability Analysis

It is observed that system (3) has at most three biologically feasible equilibrium points; namely, Ei, i=0,1,2. The existence conditions for each of these equilibrium points are discussed below:(1)
The trivial equilibrium pointsE0=0,0,0 exist always.(2)
The predator free equilibrium point is denoted byE1=y~1,y~2,0, where(7)y~1=1b1b2y~22+y~2whiley~2 is a positive root of the following third-order polynomial(8)A1y23+A2y22+A3y2+A4=0;here,A1=a2b22/b12, A2=2a2b2/b12, A3=a2+a3b1b2/b12, and A4=a3-a1b1/b1 exist uniquely in the interior of y1y2-plane if and only if the following condition holds:(9)a3<a1b1.(3)
The interior (positive) equilibrium point is given byE2=y1∗,y2∗,y3∗, where(10)y1∗=b2+b3y2∗+1-b3b4y2∗b1-b3y2∗,y3∗=y1∗+y2∗-b4whiley2∗ is a positive root of the following third-order polynomial:(11)B1y23+B2y22+B3y2+B4=0;here,(12)B1=a2b2+b32+a4b2b2+b3>0,B2=2a2b2+b3+2a4b2+a4b31-b3b4-a3b2b3+b32a3+a1,B3=2a1b1b3+a2+a41-b3b42+b1a3-a4b4b2+b3+b1a4-b3a3-a4b41-b3b4,B4=b1a3-a1b1-a3b3b4-a4b41-b3b4.Clearly, (11) has a unique positive root represented by y2∗ if the following set of conditions hold:(13)B4<0withB2>0orB3<0Therefore,E2 exists uniquely in int.  R+3 if in addition to condition (13) the following conditions are satisfied.(14)y1∗+y2∗>b4b3b4-1b2+b3<y2∗<b1b3orb1b3<y2∗<b3b4-1b2+b3.Now to study the local stability of these equilibrium points, the Jacobian matrix J(y1,y2,y3) for the system (3) at any point (y1,y2,y3) is determined as(15)-2a2y1-a3-a4y3a1-a4y1b1-2b2y2-1-b3y3-b3y2y3y3y1+y2-2y3-b4.Thus, system (3) has the following Jacobian matrix near E0=(0,0,0).(16)JE0=-a3a10b1-1000-b4.Then the characteristic equation of J(E0) is given by(17)λ+b4λ2+a3+1λ+a3-a1b1=0.Clearly, all roots of (17) have negative real parts if and only if the following condition holds:(18)a3>a1b1.So, E0 is locally asymptotically stable under condition (18) and saddle point otherwise. Therefore, E0 is locally asymptotically stable whenever E1 does not exist and unstable whenever E1 exists.The Jacobian matrix of system (3) around the equilibrium point E1=(y~1,y~2,0) reduced to(19)JE1=-2a2y~1-a3a1-a4y~1b1-2b2y~2-1-b3y~200y~1+y~2-b4=aij3×3.Then the characteristic equation of J(E1) is written by(20)λ-a33λ2-a11+a22λ+a11a22-a12a21=0.Straightforward computation shows that all roots of (20) have negative real part provided that(21)y~1+y~2<b4(22)2a2y~1+a32b2y~2+1>a1b1.So, E1 is locally asymptotically stable if the above two conditions hold.Finally the Jacobian matrix of system (3) around the interior equilibrium point E2=(y1∗,y2∗,y3∗) is written(23)JE2=bij3×3;here(24)b11=-2a2y1∗+a4y3∗+a3,b12=a1,b13=-a4y1∗,b21=b1,b22=-2b2y2∗+b3y3∗+1,b23=-b3y2∗,b31=y3∗,b32=y3∗,b33=-y3∗.Hence, the characteristic equation of J(E2) becomes(25)λ3+D1λ2+D2λ+D3=0with(26)D1=-b11+b22+b33>0,D2=b11b22-b12b21+b11b33-b13b31+b22b33-b23b32,D3=b33b12b21-b11b22+b11b23b32-b12b23b31+b22b13b31-b21b13b32.Consequently, Δ=D1D2-D3 can be written as (27)Δ=b11+b22b12b21-b11b22-b11b22b33-b11+b22b11b33-b13b31.-b222b33-b11b332-b22b332+b13b31b33+b21+b23b32b33+b12Since D1>0, then according to Routh-Hurwitz criterion E2 is locally asymptotically stable if and only if D3>0 and Δ=D1D2-D3>0. According to the form of D3 and the signs of Jacobian elements the last four terms are positive, while the first term will be positive under the sufficient condition (28) below. However Δ becomes positive if and only if in addition to condition (28) the second sufficient condition given by (29) holds.(28)2a2y1∗+a4y3∗+a32b2y2∗+b3y3∗+1>a1b1(29)y3∗>b1,y3∗>a1.Therefore under these two sufficient conditions E2 is locally asymptotically stable.


## 4. Global Stability

In this section the global stability for the equilibrium points of system (3) is investigated by using the Lyapunov method as shown in the following theorems.Theorem 2.
Assume that the vanishing equilibrium pointE0 is locally asymptotically stable; then it is globally asymptotically stable in R+3 if and only if the following condition holds:(30)a4b1a3<b3<a4a1.Proof.
Consider the following positive definite real valued function:(31)V0y1,y2,y3=1a4y1+1b3y2+y3.Straightforward computation shows that the derivative of V0 with respect to t is given by(32)dV0dt<a4b1-a3b3a4b3y1+a1b3-a4a4b3y2-b4y3.Therefore, by using condition (30), we obtain dV0/dt which is negative definite in R+3, and then V0 is a Lyapunov function with respect to E0. Hence E0 is globally asymptotically stable in R+3 and the proof is complete.Theorem 3.
Assume that the predator free equilibrium pointE1=(y~1,y~2,0) is locally asymptotically stable; then it is globally asymptotically stable in R+3 if the following condition holds:(33)a1a4y1+b1b3y22<4a1y~2a4y1y~1+a2a4b1y~1b3y2y~2+b2b3.Proof.
Consider the following positive definite real valued function:(34)V1y1,y2,y3=1a4y1-y~1-y~1ln⁡y1y~1+1b3y2-y~2-y~2ln⁡y2y~2+y3.Straightforward computation shows that the derivative of V1 with respect to t is given by(35)dV1dt=-a1y~2a4y1y~1+a2a4y1-y~12-b1y~1b3y2y~2+b2b3y2-y~22+a1a4y1+b1b3y2y1-y~1y2-y~2-y32-b4-y~1-y~2y3.Now using condition (33) gives us that(36)dV1dt<-a1y~2a4y1y~1+a2a4y1-y~1-b1y~1b3y2y~2+b2b3y2-y~22-b4-y~1-y~2y3.Clearly dV1/dt is negative definite due to local stability condition (21). Hence V1 is a Lyapunov function with respect to E1, and then E1 is globally asymptotically stable, which completes the proof.Theorem 4.
Assume that the interior equilibrium pointE2=(y1∗,y2∗,y3∗) is locally asymptotically stable in R+3; then it is globally asymptotically stable if and only if the following condition holds:(37)y2∗a1b3y1∗b1-a42<4y2∗a1a2b3b4y1∗b1.Proof.
Consider the following positive definite real valued function aroundE2:(38)V1y1,y2,y3=y1-y1∗-y1∗ln⁡y1y1∗+y2∗a1y1∗b1y2-y2∗-y2∗ln⁡y2y2∗+y2∗a1b3y1∗b1y3-y3∗-y3∗ln⁡y3y3∗.Our computation for the derivative of V2 with respect to t gives that(39)dV2dt=-a1y1y2y1∗y2y1∗-y1y2∗2-a2y1-y1∗2-y2∗a1b2y1∗b1y2-y2∗2+y2∗a1b3y1∗b1-a4y1-y1∗y3-y3∗-y2∗a1b3b4y1∗b1y3-y3∗2.Now by using the condition (37) we obtain that(40)dV2dt<-a1y1y2y1∗y2y1∗-y1y2∗2-y2∗a1b2y1∗b1y2-y2∗2-a2y1-y1∗-y2∗a1b3b4y1∗b1y3-y3∗2.According to the above inequality we have dV2/dt which is negative definite; therefore, E2 is globally asymptotically stable in R+3 and hence the proof is complete.


## 5. Local Bifurcation

In this section the local bifurcation near the equilibrium points of system (3) is investigated using Sotomayor’s theorem for local bifurcation [24]. It is well known that the existence of nonhyperbolic equilibrium point is a necessary but not sufficient condition for bifurcation to occur. Now rewrite system (3) in the form(41)Y˙=fY,where(42)Y=y1,y2,y3T,f=f1,f2,f3T.Then according to Jacobian matrix of system (3) given in (15), it is simple to verify that for any nonzero vector U=(u1,u2,u3)T we have(43)D2fy1,y2,y3U,U=∂2f∂y12u12+∂2f∂y1∂y2u1u2+∂2f∂y2∂y1u2u1+∂2f∂y22u22+∂2f∂y1∂y3u1u3+∂2f∂y3∂y1u3u1+∂2f∂y2∂y3u2u3+∂2f∂y3∂y2u3u2+∂2f∂y32u32.Consequently, we obtain that(44)D2fy1,y2,y3U,U=-2a2u12-2a4u1u3-2b2u22-2b3u2u32u1u3+2u2u3-2u32and therefore(45)D3fy1,y2,y3U,U,U=000.Thus system (3) has no pitchfork bifurcation due to (45). Moreover, the local bifurcation near the equilibrium points is investigated in the following theorems:Theorem 5.
System (3) undergoes a transcritical bifurcation near the vanishing equilibrium point, but saddle node bifurcation cannot occur, when the parameter a3 passes through the bifurcation value a3∗=a1b1.Proof.
According to the Jacobian matrixJ(E0) given by (16), system (3) at the equilibrium point E0 with a3=a3∗ has zero eigenvalue, say λ0∗=0, and the Jacobian matrix J(E0,a3∗) becomes(46)JE0,a3∗=J0∗=-a3∗a10b1-1000-b4.Now let U[0]=(u1[0],u2[0],u3[0])T be the eigenvector corresponding to the eigenvalue λ0∗=0. Thus J0∗U[0]=0 gives U[0]=(u1[0],b1u1[0],0)T, where u1[0] represents any nonzero real number. Also, let W[0]=(w1[0],w2[0],w3[0])T represents the eigenvector corresponding to the eigenvalue λ0∗=0 of J0∗T. Hence J0∗TW[0]=0 gives that W[0]=w1[0],a1w1[0],0T, where w1[0] denotes any nonzero real number. Now, since(47)dfda3=fa3Y,a3=df1da3,df2da3,df3da3T=-y1,0,0Tthus fa3E0,a3∗=(0,0,0)T, which gives W[0]Tfa3(E0,a3∗)=0. So, according to Sotomayor’s theorem for local bifurcation, system (3) has no saddle node bifurcation at a3=a3∗. Also, since(48)Dfa3E0,a3=-100000000then,(49)W0TDfa3E0,a3∗U0=w10,a1w10,0-u10,0,0T=-u10w10≠0.Moreover, by substituting E0,a3∗, and U[0] in (44) we get that(50)D2fE0,a3∗U0,U0=-2a2u102,-2b2b12u102,0T.Hence, it is obtain that(51)W0TD2fE0,a3∗U0,U0=-2a2+a1b2b12u102w10≠0.Thus, according to Sotomayor’s theorem system (3) has a transcritical bifurcation at E0 as the parameter a3 passes through the value a3∗; thus the proof is complete.Theorem 6.
Assume that condition (22) holds; then system (3) undergoes a transcritical bifurcation near the predator free equilibrium point E1, but saddle node bifurcation cannot occur, when the parameter b4 passes through the bifurcation value b4∗=y~1+y~2.Proof.
According to the Jacobian matrixJ(E1) given by (19), system (3) at the equilibrium point E1 with b4=b4∗ has zero eigenvalue, say λ1∗=0, and the Jacobian matrix J(E1,b4∗) becomes(52)JE1,b4∗=J1∗=aij∗3×3,where aij∗=aij;∀i,j=1,2,3, with a33∗=0. Now let, U[1]=(u1[1],u2[1],u3[1])T be the eigenvector corresponding to the eigenvalue λ1∗=0. Thus J1∗U[1]=0 gives U[1]=(Λ1u3[1],Λ2u3[1],u3[1])T, where Λ1=a12a23-a22a13/a11a22-a12a21 and Λ2=a21a13-a11a31/a11a22-a12a21 are negative according to the sign of the Jacobian elements and u3[1] represents any nonzero real numbers. Also, let W[1]=(w1[1],w2[1],w3[1])Trepresent the eigenvector corresponding to eigenvalue λ1∗=0 of J1∗T. Hence J1∗TW[1]=0 gives that W[1]=(0,0,w3[1])T, where w3[1] stands for any nonzero real numbers. Now, since(53)dfdb4=fb4Y,b4=df1db4,df2db4,df3db4T=0,0,-y3Tthus fb4(E1,b4∗)=(0,0,0)T, which gives W[1]Tfb4(E1,b4∗)=0. So, according to Sotomayor’s theorem for local bifurcation, system (3) has no saddle node bifurcation at b4=b4∗. Also, since(54)Dfb4E1,b4=00000000-1then, we can have(55)W1TDfb4E1,b4∗U1=0,0,w310,0,-u31T=-u31w31≠0.Moreover, substituting E1,b4∗, and U[1] in (44) gives(56)D2fE1,b4∗U1,U1=2u312-a2Λ12-a4Λ1,-b2Λ22-b3Λ2,Λ1+Λ2-1T.Hence, it is obtained that(57)W1TD2fE1,b4∗U1,U1=2Λ1+Λ2-1u312w31≠0.Thus, according to Sotomayor’s theorem system (3) has a saddle node bifurcation at E1 as the parameter b4 passes through the value b4∗; thus the proof is complete.Theorem 7.
Assume that condition (21) holds; then system (3) undergoes a saddle node bifurcation near the predator free equilibrium point E1 when the parameter a1 passes through the bifurcation value a1∗=(2a2y~1+a3)(2b2y~2+1)/b1.Proof.
According to the Jacobian matrixJ(E1) given by (19), system (3) at the equilibrium point E1 with a1=a1∗ has zero eigenvalue, say λ1∗∗=0, and the Jacobian matrix J(E1,a1∗) becomes(58)JE1,a1∗=J1∗∗=aij∗∗3×3,where aij∗∗=aij;∀i,j=1,2,3, with a12∗∗=a1∗. Now let U[11]=(u1[11],u2[11],u3[11])T be the eigenvector corresponding to the eigenvalue λ1∗∗=0. Thus J1∗∗U[11]=0 gives U[11]=(-a1∗/a11u2[11],u2[11],0)T, where u2[11] represents any nonzero real numbers. Also, let W[11]=(w1[11],w2[11],w3[11])T represent the eigenvector corresponding to eigenvalue λ1∗∗=0 of J1∗∗T. Hence J1∗∗TW[11]=0 gives that W[11]=(-a21/a11w2[11],w2[11],Ψ/a11a33w2[11])T, where Ψ=a13a21-a11a23 is negative due to the sign of the Jacobian elements and w2[11] denotes any nonzero real numbers. Now, since(59)dfda1=fa1Y,a1=df1da1,df2da1,df3da1T=y2,0,0Tthus fa1(E1,a1∗)=(y~2,0,0)T; hence W[11]Tfa1(E1,a1∗)=-a21/a11w2[11]y~2≠0. So, according to Sotomayor’s theorem for local bifurcation the first condition of saddle node bifurcation is satisfied in system (3) at a1=a1∗. Moreover, substituting E1,a1∗, and U[11] in (44) gives(60)D2fE1,a1∗U11,U11=u2112-2a2a1∗2a112,-2b2,0T.Hence, it is obtained that(61)W11TD2fE1,a1∗U11,U11=2a2a1∗2a21a113-b2u2112w211≠0.Thus, according to Sotomayor’s theorem system (3) has a transcritical bifurcation at E1 as the parameter a1 passes through the value a1∗; thus the proof is complete.Theorem 8.
Assume that(62)2a2y1∗+a4y3∗+a32b2y2∗+b3y3∗+1<a1b1(63)y1∗<a1a4.Then system (3) undergoes a saddle node bifurcation near the interior equilibrium point E2, as the parameter b1 passes through the bifurcation value b1∗=Γ1/Γ2, where Γ1 and Γ2 are given in the proof.Proof.
According to the determinant of the Jacobian matrixJ(E2) given by D3 in (25), condition (62) represents a necessary condition to have nonpositive determinant for J(E2). Now rewrite the form of the determinant as follows:(64)D3=Γ1-b1Γ2.Here Γ1=b11b23b32+b22b31b13-b11b22b33-b12b23b31 and Γ2=b13b32-b33b12. Obviously, Γ1 is positive always, while Γ2 is positive under the condition (63). Thus it is easy to verify that D3=0 and hence J(E2) has zero eigenvalue, say λ2∗=0, as b1 passes through the value b1∗=Γ1/Γ2, which means that E2 becomes a nonhyperbolic point. Let now the Jacobian matrix of system (3) at E2 with b1=b1∗ be given by(65)JE2,b1∗=J2∗=bij∗3×3,where bij∗=bij∀i,j=1,2,3, and b21∗=b1∗.
LetU[∗]=(u1[∗],u2[∗],u3[∗])T be the eigenvector corresponding to the eigenvalue λ2∗=0. Thus J2∗U[∗]=0 gives U[∗]=Φ1u3[∗],Φ2u3[∗],u3[∗]T, where Φ1=b12b23-b22b13/b11b22-b12b1∗ and Φ2=b1∗b13-b11b23/b11b22-b12b1∗ are positive due to the Jacobian elements and u3[∗] represents any nonzero real number. Also, let W[∗]=(w1[∗],w2[∗],w3[∗])T represent the eigenvector corresponding to eigenvalue λ2∗=0 of J2∗T. Hence J2∗TW[∗]=0 gives that W[∗]=Θ1w3[∗],Θ2w3[∗],w3[∗]T, where Θ1=b21b32-b22b31/b11b22-b12b1∗ and Θ2=b12b31-b11b32/b11b22-b12b1∗ are negative due to the Jacobian elements and w3[∗] denotes any nonzero real numbers. Now, since(66)dfdb1=fb1Y,b1=df1db1,df2db1,df3db1T=0,y1,0Tthus fb1(E2,b1∗)=(0,y1∗,0)T, which gives W[∗]Tfb1(E2,b1∗)=Θ2y1∗w3[∗]≠0. Consequently the first condition of saddle node bifurcation is satisfied. Moreover, by substituting E2,b1∗, and U[∗] in (44) we get that(67)D2fE2,b1∗U∗,U∗=2u3∗2-a2Φ12-a4Φ1,-b2Φ22-b3Φ2,Φ1+Φ2-1T.Hence, it is obtained that(68)W∗TD2fE2,b1∗U∗,U∗=2u3∗2w3∗×-a2Φ12+a4Φ1Θ1-b2Φ22+b3Φ2Θ2+Φ1+Φ2-1≠0.So, according to Sotomayor’s theorem, system (3) has a saddle node bifurcation as b1 passes through the value b1∗ and hence the proof is complete.


## 6. Numerical Simulations

In this section, the global dynamics of system (3) is investigated numerically. The objectives first confirm our obtained analytical results and second specify the control set of parameters that control the dynamics of the system. Consequently, system (3) is solved numerically using the following biologically feasible set of hypothetical parameters with different sets of initial points and then the resulting trajectories are drawn in the form of phase portrait and time series figures.(69)a1=2,a2=0.1,a3=0.4,a4=0.5,b1=0.4,b2=0.1,b3=0.5,b4=0.2.Clearly, Figure1 shows the asymptotic approach of the solutions, which started from different initial points to a positive equilibrium point (0.46, 0.15, 0.42), for the data given by (69). This confirms our obtained result regarding the existence of globally asymptotically stable positive point of system (3) provided that certain conditions hold.Figure 1
3D phase plot of the system (3) for the data given by (69) starting from different initial values in which the solution approaches asymptotically to 0.46, 0.15, and 0.42.Now in order to discuss the effect of the parameters values of system (3) on the dynamical behavior of the system, the system is solved numerically for the data given in (69) with varying one parameter each time. It is observed that varying parameters values a2,a4,b2,b3, and b4 have no qualitative effect on the dynamical behavior of system (3) and the system still approaches to a positive equilibrium point. On the other hand, when a1 decreases in the range (a1≤1.05) keeping other parameters fixed as given in (69) the dynamical behavior of system (3) approaches asymptotically to the vanished equilibrium point as shown in the typical figure given by Figure 2. Similar observations have been obtained on the behavior of system (3) in case of increasing the parameter a3 in the range (a3≥0.8) or decreasing the parameter b1 in the range (b1≤0.2), with keeping other parameters fixed as given in (69), and then the solution of system (3) is depicted in Figures 3 and 4, respectively. Finally, for the parameters a2=0.9 and b4=0.7 with other parameters fixed as given in (69), the solution of system (3) approaches asymptotically to the predator free equilibrium point as shown in the Figure 5.Figure 2
Time series of system (3) approaches asymptotically to the vanishing equilibrium point for a1=0.8 with other parameters given by (69).Figure 3
Time series of system (3) approaches asymptotically to the vanishing equilibrium point for a3=0.85 with other parameters given by (69).Figure 4
Time series of system (3) approaches asymptotically to the vanishing equilibrium point for b1=0.15 with other parameters given by (69).Figure 5
Time series of system (3) approaches asymptotically to the predator free equilibrium point E1=(0.42,0.16,0) for a2=0.9 and b4=0.7 with other parameters given by (69).


## 7. Discussion

In this paper, a model that describes the prey-predator system having a refuge and stage structure properties in the prey population has been proposed and studied analytically as well as numerically. Sufficient conditions which ensure the stability of equilibria and the existence of local bifurcation are obtained. The effect of each parameter on the dynamical behavior of system (3) is studied numerically and the trajectories of the system are drowned in the typical figures. According to these figures, which represent the solution of system (3) for the data given by (69), the following conclusions are obtained.(1)
System (3) has no periodic dynamics rather than the fact that the system approaches asymptotically to one of their equilibrium points depending on the set of parameter data and the stability conditions that are satisfied.(2)
Although the position of the positive equilibrium point in the interior ofR+3 changed as varying in the parameters values a2,a4,b2,b3, and b4, there is no qualitative change in the dynamical behavior of system (3) and the system still approaches to a positive equilibrium point. Accordingly adding the refuge factor which is included implicitly in these parameters plays a vital role in the stabilizing of the system at the positive equilibrium point.(3)
Decreasing in the value of growth rate of immature prey or in the value of conversion rate from immature prey to mature prey keeping the rest of parameter as in (69) leads to destabilizing of the positive equilibrium point and the system approaches asymptotically to the vanishing equilibrium point, which means losing the persistence of system (3).(4)
Increasing in the value of grownup rate of immature prey keeping the rest of parameter as in (69) leads to destabilizing of the positive equilibrium point and the system approaches asymptotically to the vanishing equilibrium point too, which means losing the persistence of system (3).(5)
Finally, for the data given by (69), increasing intraspecific competition of immature prey and natural death rate of the predator leads to destabilizing of the positive equilibrium point and the solution approaches instead asymptotically to the predator free equilibrium point, which confirm our obtained analytical results represented by conditions (21)-(22).



---

*Source: 2010464-2016-11-24.xml*