# Stability and Bifurcation of a Prey-Predator-Scavenger Model in the Existence of Toxicant and Harvesting

Authors: Huda Abdul Satar; Raid Kamel Naji
Journal: International Journal of Mathematics and Mathematical Sciences
 (2019)

---

## Abstract

In this paper a prey-predator-scavenger food web model is proposed and studied. It is assumed that the model considered the effect of harvesting and all the species are infected by some toxicants released by some other species. The stability analysis of all possible equilibrium points is discussed. The persistence conditions of the system are established. The occurrence of local bifurcation around the equilibrium points is investigated. Numerical simulation is used and the obtained solution curves are drawn to illustrate the results of the model. Finally, the nonexistence of periodic dynamics is discussed analytically as well as numerically.


---

## Body

## 1. Introduction

Food web model is one of natural ecological systems that describe the interaction among ‎the species for food. It is well known that in the real world environment there is no species living separately from other species, rather all the species living together interact with each other in different ways. The interaction of species may be competition for food, mutualism, prey-predator, etc. Different types of prey-predator models including different biological factors were studied by many researchers. See [1–5] and the references cited therein. Later on several three species food web models have been proposed in which a third species is included within prey-predator model so that it stands for second prey, second predator, omnivore, cannibalism, etc. See [6–10] and the references cited therein. All these models were formulated according to the Lotka-Volterra framework model.On the other hand, few works have considered in their formulation a scavenger species [11–13]. A scavenger is an animal that consumed carcasses of other animals, those which are dead naturally or killed by other animals. Accordingly study of such type of animals is important due to their support of the environment’s cleaning. Recently, Abdul Satar and Naji [14] proposed and studied an ecological model consisting of prey, predator, and scavenger involving Michaelis–Menten type of harvesting function. They studied the stability, persistence, and local bifurcation of the proposed model analytically as well as numerically. It is observed that the proposed model is very sensitive for varying in their parameters values especially those related with scavenger and undergoes different types of local bifurcation.Clearly, food and energy are very important things for human populations; therefore exploitation of these resources has increased significantly with time. Consequently the harvesting problem for multispecies food web systems becomes an important subject for many researchers nowadays [15–17]. On the other hand, the effects of toxicants on ecological communities have become problems of major environmental concern due to the direct effects of food web systems, which in turn affects the human life. Accordingly, the effects of toxicants on the ecological systems are studied by many researches [18–22]. Therefore problems related to harvesting and toxicants in ecological systems have been drawing attention of many researchers and they are considered the main objective of this research work too.‎Keeping the above in view, in this paper, a food web model consisting of prey, predator, and scavenger is proposed and studied. It is assumed that the effects of harvesting and toxicant on the model species are also included. The next section is concerned with the model formulation and discusses the existence, uniqueness, and boundedness of the solution. The stability and persistence conditions are established in Section Three. Section Four is concerned with the study of the bifurcation in the system. Numerical simulation is carried out in Section Five and the last section included the discussion and conclusions.


## 2. Model Formulation

In this section a mathematical model formulation is carried out according to the following set of hypotheses:(1)
The prey species in the absence of predation, harvesting, and toxicant is assumed to be growth logistically with intrinsic growth rater>0 and carrying capacity K>0.(2)
All the species are assumed to be subjected to a combined harvesting effortE>0, with positive catchability coefficients ci;i=1,2,3 for the prey, predator, and scavenger respectively.(3)
Lotka-Volterra types of functional responses are used for describing the predation and scavenge process withαi;1,2,3 as positive maximum attack rates, while the positive constants e1 and e2 stand for the conversion rates of prey to predator and scavenger, respectively.(4)
It is assumed that the predator and scavenger decay exponentially in the absence of prey species with natural death ratesd1>0 and d2>0, respectively.(5)
Finally it is assumed that the prey is directly infected by some external toxic substance such as industrial wastes, while the individuals of predator and scavengers are indirectly affected by the toxic substance due to feeding on this infected prey or carcass of predator. Hence the positive parametersγi;i=1,2,3 represent the coefficients of toxicity of prey, predator, and scavenger, respectively, with γ2<γ3<γ1 due to the fact that the effects of toxic substance are direct on the prey and indirect of the other species; also the scavenger will be under the effects of the toxic more than the predator due to their feeding on both the infected prey and the carcass of infected predator.According to these hypotheses, whenXT, Y(T), and Z(T) represent the density of prey, predator, and scavenger at time T, the dynamics of the prey-predator-scavenger model can be described using the following set of differential equations:(1)dXdT=rX1-XK-α1XY-α2XZ-c1EX-γ1X3dYdT=e1α1XY-d1Y-c2EY-γ2Y2dZdT=e2α2XZ+α3YZ-d2Z-c3EZ-γ3Z2with X0≥0, Y0≥0, and Z0≥0. Note that since the prey is assumed to be directly affected by toxic substance, the other species are indirectly affected by this substance due to their feeding on the prey. Therefore, a cubic term (γ1X3) is used in the first equation while quadratic term (γ2Y2 and γ3Z2) is used in the other equations of the system and that is due to the mathematical fact that there is an accelerated growth in the production of toxic substance to the prey species more than that of the other species in the system.Now, in order to study the above system of equations more generally, we drop all the units from it by using the following dimensionless variables and constants:(2)x=XK,y=α1rY,z=α2rZ,t=rT,α=c1Er,β=γ1K2r,γ=e1α1Kr,δ=d1r,η=c2Er,μ=γ2α1,ν=e2α2Kr,ω=α3α1,ζ=d2r,ε=c3Er,τ=γ3α2Here the dimensionless parameters α,β,γ,δ,η,μ,ν,ω,ζ,ε, and τ stand for the ratio of harvesting of prey to intrinsic growth rate of prey species, ratio of toxicity of prey to intrinsic growth rate of prey, ratio of prey’s attack from predator to intrinsic growth rate of prey, ratio of predator’s death rate to the intrinsic growth rate of the prey, ratio of harvesting of predator to intrinsic growth rate of prey species, ratio of toxicity coefficient of predator to prey’s attack rate from predator, ratio of prey’s attack rate from scavenger to intrinsic growth rate of prey, ratio of predator’s attack rate from scavenger to the prey’s attack rate from predator, ratio of scavenger death rate to the intrinsic growth rate of the prey, ratio of harvesting of scavenger to intrinsic growth rate of prey species, and the ratio of toxicity coefficient of scavenger to prey’s attack rate from scavenger, respectively. Accordingly, the following dimensionless system is obtained:(3)dxdt=x1-x-y-z-α-βx2=xf1x,y,zdydt=yγx-δ-η-μy=yf2x,y,zdzdt=zνx+ωy-ζ-ε-τz=zf3x,y,zTherefore, system (3) has the following domain:(4)R+3=x,y,z∈R3∣x≥0,y≥0,z≥0Clearly, all the interaction functions in the right-hand side are continuous and have continuous partial derivatives, and hence the solution of system (3) exists and is unique. Note that according to the first equation of system (3) we have always that (5)1-α>0Otherwise all the species will go extinct and hence from now onward condition (5) is assumed to be always satisfied. Moreover, in order to guarantee the convergent of the solution to an attractor, the solution of system (3) is proved to be uniformly bounded as shown in the following theorem.Theorem 1.
All solutions of the system (3) with the initial condition belonging to R+3 are uniformly bounded in the region(6)Λ=x,y,z∈R+3:x≤1;x+yγ≤l1;x+yγ+zν≤l2Here l1 and l2 are given in the proof.Proof.
From the first equation we havedx/dt≤x1-x; then it is easy to verify that for t→∞ we get x≤1. Now let W1=x+1/γy; then the derivative of W1 with respect to time can be written as (7)dW1dt≤2-α-M1W1Here M1=min.1,δ+η. Then, according to the above differential inequality, it is easy to verify that for t→∞; we obtain W1≤2-α/M1=l1.
Now assume thatW2=x+y/γ+z/ν; then the derivative of W2 with respect to time can be written as (8)dW2dt≤2-α-M2W2Here M2=min.1,δ+η,ζ+ε. Then direct computation shows that, for t→∞, we obtain W2≤2-α/M2=l2. Hence the proof is complete.Keeping the above in view, system (3) describes the dynamics of a food web model consisting of prey, predator, and scavenger involving the effects of two ecological factors harvesting and toxicant. It has a unique solution moving in the interior of bounded region Λ belonging to the positive octant R+3. Therefore it is important to investigate the asymptotic behavior of the solution of system (3) and its attracting sets as time increases and this is the main objective of Section 3. Furthermore the effects of varying the parameters values on the asymptotic behavior of the system are important from the biological point of view to complete the study of any mathematical system and that will be the main objective of Section 4. Finally in Section 5 we will try to confirm our obtained analytical results with the help of numerical simulation using hypothetical set of parameters values.


## 3. Stability Analysis with Persistence

In the following the conditions of all possible equilibrium points of system (3) are established and then the local stability analysis of each of them is studied. Straightforward computation shows that system (3) has at most the following equilibrium points.The trivial equilibrium pointE0=0,0,0 and the axial equilibrium point E1=x-,0,0 always exist; here (9)x-=-12β+12β1-4βα-1The scavenger free equilibrium point located in the interior of positive quadrant of xy-plane and denoted by E2=(x~,y~,0) exists uniquely provided that(10a)x~>δ+ηγwhile(10b)x~=-μ+γ2μβ+12μβμ+γ2-4μβμα-μ+δ+ηy~=γx~-δ+ημThe predator free equilibrium point located in the interior of positive quadrant of xz-plane and denoted by E3=(x^,0,z^) exists uniquely provided that(11a)x^>ζ+ενwhile(11b)x^=-τ+ν2τβ+12τβτ+ν2-4τβτα-τ+ζ+εz^=νx^-ζ+ετFinally, the coexistence or positive equilibrium point is located at the interior of positive octant R+3 and denoted by E4=(x∗,y∗,z∗), where(12a)x∗=-A22A1+12A1A22-4A1A3y∗=γx∗-δ+ημz∗=μν+ωγx∗-ωδ+η+μζ+ετμHerex∗ is a positive root of the following quadratic polynomial equation:(12b)A1x2+A2x+A3=0where(13)A1=βτμ>0A2=τμ+γ+μν+ωγ>0A3=τμα-1-μζ+ε-δ+ητ+ω<0Clearly,E4=(x∗,y∗,z∗) exists uniquely provided that the following conditions hold:(14a)x∗>δ+ηγ(14b)μν+γx∗>ωδ+η+μζ+εNow the local stability analysis of these equilibrium points is investigated using the linearization technique by computing the Jacobian matrix of the system at each of these points.Straightforward computation shows that the Jacobian matrix of system (3) around the trivial equilibrium point E0=(0,0,0) has the following eigenvalues:(15)λ01=1-α>0,λ02=-δ+η,λ03=-ζ+εAccordingly, the origin is always unstable saddle point. Note that if we violate condition (5) then the trivial equilibrium point transfers to asymptotically stable point.The Jacobian matrix evaluated at the axial equilibrium pointE1=(x-,0,0) can be written as(16a)JE1=-x--2βx-2-x--x-0γx--δ+η000νx--ζ+εTherefore the eigenvalues are given by(16b)λ11=-x-+2βx-2,λ12=γx--δ+η,λ13=νx--ζ+εHence, the axial equilibrium point E1=(x-,0,0) is locally asymptotically stable if the following condition holds:(17)x-<min⁡.δ+ηγ,ζ+ενNow the Jacobian matrix evaluated at the scavenger free equilibrium point E2=(x~,y~,0) is given by(18a)JE2=-x~-2βx~2-x~-x~γy~-μy~000νx~+ωy~-ζ+εClearly, one of the eigenvalues of J(E2) is given by(18b)λ21=T22+T22-4D22;λ22=T22-T22-4D22;λ23=νx~+ωy~-ζ+εwhere T2=-(x~+2βx~2+μy~)<0 and D2=x~+2βx~2μy~+γx~y~>0. Therefore these eigenvalues have negative real parts and hence E2 is locally asymptotically stable if and only if the following condition holds:(19)ζ+ε>νx~+ωy~Further, the Jacobian matrix evaluated at the predator free equilibrium point E3=(x^,0,z^) is given by(20a)JE3=-x^-2βx^2-x^-x^0γx^-δ+η0νz^ωz^-τz^Consequently, the eigenvalues of J(E3) are given by(20b)λ31=T32+T32-4D32;λ32=γx^-δ+η;λ33=T32-T32-4D32where T3=-(x^+2βx^2+τz^)<0 and D3=x^+2βx^2τz^+νx^z^>0. Hence all these eigenvalues have negative real parts and hence E3 is locally asymptotically stable if and only if the following condition holds:(21)δ+η>γx^Finally, the Jacobian matrix evaluated at the coexistence equilibrium point E4=(x∗,y∗,z∗) is given by (22a)JE4=-x∗-2βx∗2-x∗-x∗γy∗-μy∗0νz∗ωz∗-τz∗=aijThe characteristic equation associated with J(E4) is given by(22b)λ43+A1λ42+A2λ4+A3=0where(23)A1=-a11+a22+a33A2=a11a22-a12a21+a22a33+a11a33-a13a31A3=-a33a11a22-a12a21-a13a22a31-a21a32with(24)∆=A1A2-A3=-a11+a22a11a22-a12a21-a11a33a11+a22+a33-a22a33a22+a33-a11a22a33-a12a21+a13a31a33+a21a32Consequently, the local stability of the coexistence equilibrium point of system (3) can be discussed in the following theorem.Theorem 2.
The coexistence equilibrium point of system (3) is locally asymptotically stable in the interior of R+3 provided that the following sufficient condition holds:(25)ντz∗>γωy∗Proof.
According to the Routh-Hurwitz criterion Equation (22b) has three roots with negative real parts provided that A1>0; A3>0 and ∆=A1A2-A3>0. Straightforward computation shows that (26)A1=x∗+2βx∗2+μy∗+τz∗>0A3=x∗y∗z∗τμ+2μβx∗+γ+μν+γω>0∆=x∗y∗x∗1+2βx∗+μy∗μ1+2βx∗+γ+τx∗z∗1+2βx∗x∗1+2βx∗+μy∗+τz∗+μτy∗z∗μy∗+τz∗+x∗y∗1+2βx∗μτz∗+γx∗+x∗z∗ντz∗-γωy∗Clearly the given sufficient condition (25) leads to ∆>0. Hence all the eigenvalues of JE4 have negative real parts and hence the proof is followed.Now before we go further to study the global stability of system (3), the persistence of the system is investigated. Recall that the system is persistent if and only if all the species coexist for all the future time. This is satisfied provided that the trajectories of system (3), those started in the positive cone (Int.R+3), are eventually bounded away from the boundary planes. Therefore the persistence of the continuous biological systems means the survival of all the interacting populations. Moreover, since the omega limit set in the boundary planes of any orbit started at a point in the interior of positive cone can be either an equilibrium point or periodic dynamic that belongs to one of boundary planes, therefore the investigation of the possible subsystems is important. Consequently, in order to investigate the persistence of system (3), the existence of periodic dynamics in the boundary planes is investigated as shown in the following theorems.Clearly system (3) has two subsystems; the first subsystem is obtained in the absence of scavenger and located in the positive quadrant of xy-plane while the second subsystem is obtained in the absence of predator and located in the positive quadrant of xz-plane. These subsystems can be written as follows:(27)dxdt=x1-x-y-α-βx2=xf11x,ydydt=yγx-δ-η-μy=yf12x,yand(28)dxdt=x1-x-z-α-βx2=xf21x,zdzdt=zνx-ζ-ε-τz=zf22x,zObviously, subsystem (27) has a unique positive equilibrium point which belongs to interior of positive quadrant of xy-plane given by e2=(x~,y~)≡E2. However, subsystem (28) has a unique positive equilibrium point which belongs to interior of positive quadrant of xz-plane given by e3=(x^,z^)≡E3.Theorem 3.
There is no periodic dynamics in the interior of positive quadrant ofxy-plane surrounding a unique equilibrium point e2=(x~,y~).Proof.
DefineBx,y=1/xy that obviously satisfies Bx,y>0 and continuously differentiable function in the interior of positive quadrant of xy-plane. Now, we have(29)∆=∂Bf11∂x+∂Bf12∂y=-1y-2βxy-μx<0Clearly ∆ is not equal to zero and it is not change sign in the interior of positive quadrant of xy-plane. So according to the Dulac-Bendixson’s criterion there is no closed curve surrounding a unique equilibrium point e2=(x~,y~).Now according to the above theorem along with the Poincare-Bendixson theorem, the equilibrium pointe2=(x~,y~)≡E2 is a globally asymptotically stable in the interior of positive quadrant of the xy-plane if and only if it is locally asymptotically stable. Similarly we can proof the following theorem.Theorem 4.
There is no periodic dynamics in the interior of positive quadrant ofxz-plane surrounding a unique equilibrium point e3=(x^,z^).Again we can easily conclude thate3=(x^,z^)≡E3 is globally asymptotically stable in the interior of positive quadrant of the xz-plane if and only if it is locally asymptotically stable. Keeping the above in view we can go further to give the following theorem.Theorem 5.
The system (3) is uniformly persistent in the interior of R+3 provided that the following conditions hold:(30a)x->max⁡.δ+ηγ,ζ+εν(30b)ζ+ε<νx~+ωy~(30c)δ+η<γx^Proof.
Consider the functionσx,y,z=xp1yp2zp3, where pi;i=1,2,3 are undetermined positive constants. Obviously σx,y,z is a C1 positive function defined in the interior of R+3.
Consequently we obtain(31)Φx,y,z=σ′x,y,zσx,y,z=p1f1+p2f2+p3f3where fi;i=1,2,3 are given in the right-hand side of system (3); therefore(32)Φx,y,z=p11-x-y-z-α-βx2+p2γx-δ-η-μy+p3νx+ωy-ζ-ε-τzNow, since there are no periodic attractors in the boundary planes, then the only possible omega limit sets of system (3) are the equilibrium points E0, E1, E2,andE3. Thus the system is uniformly persistent if we can proof that Φ(∙)>0 at each of these points. Also we have(33)ΦE0=p11-α-p2δ+η-p3ζ+εΦE1=p2γx--δ+η+p3νx--ζ+εΦE2=p3νx~+ωy~-ζ+εΦE3=p2γx^-δ+ηObviously, ΦE0>0 for all values of p1 is sufficiently large enough with respect to pi;i=2,3.
Φ E 1 > 0 for any positive constants pi;i=1,2,3 provided that the condition (30a) is satisfied. However, ΦE2andΦE3 are positive provided that the conditions (30b) and (30c) are satisfied, respectively. Then the system (3) is uniformly persistent.In the following theorems, the global stability of all the locally stable equilibrium points is studied with the help of Lyapunov method and uniformly bounded of the solution of system (3), which guarantees the existence of supremum value for each variable in the system.Theorem 6.
Let the axial equilibrium pointE1=(x-,0,0) be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient condition holds:(34)ω2ν<min.μγ,τνProof.
Consider the following positive definite real valued function:(35)L1=b1∫x-xS-x-Sds+b2y+b3zwhere bi;i=1,2,3 are positive constants to be determined. Then the derivative of this function with respect to time can be written as(36)dL1dt=b1x-x-xdxdt+b2dydt+b3dzdtHence by using system (3) with some algebraic manipulations we obtain that(37)dL1dt=-b11+βx+x-x-x-2-b1-b2γxy-b2δ+η-b1x-y-b1-b3νxz-b3ζ+ε-b1x-z-b3ω2y-b3ω2z2-b2μ-b3ω2y2-b3τ-b3ω2z2Therefore, by choosing the values of bi;i=1,2,3, so that b1=1, b2=1/γ, and b3=1/ν, we get(38)dL1dt=-1+βx+x-x-x-2-δ+ηγ-x-y-ζ+εν-x-z-ω2νy-ω2νz2-μγ-ω2νy2-τν-ω2νz2Consequently, dL1/dt is negative definite under the local stability condition (17) and condition (34). Moreover it is clear that the function L1 is radially unbounded; then according to the Lyapunov first theorem E1 is a globally asymptotically stable point. Hence the proof is complete.Theorem 7.
Let the scavenger free equilibrium pointE2=(x~,y~,0) be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient condition holds:(39)ζ+εν≥x~+ωsup.yνProof.
Consider the following positive definite real valued function:(40)L2=c1∫x~xS-x~Sds+c2∫y~yS-y~Sds+c3z(41)dL2dt=-c11+βx+x~x-x~2-c1-c2γx-x~y-y~-c2μy-y~2-c1-c3νxz-c3ζ+ε-c1x~z+c3ωyz-c3τz2Therefore, by choosing the values of ci;i=1,2,3, so that c1=1, c2=1/γ, and c3=1/ν, we get(42)dL2dt=-1+βx+x~x-x~2-μγy-y~2-ζ+εν-x~z+ωνyz-τνz2Now by using the supremum value of y we obtain that(43)dL2dt≤-1+βx+x~x-x~2-μγy-y~2-ζ+εν-x~-ωνsup.⁡yz-τνz2Clearly, condition (39) guarantees that dL2/dt<0. Moreover as the function L2 is radially unbounded then according to the Lyapunov first theorem E2 is a globally asymptotically stable point. Hence the proof is complete.Theorem 8.
Let the predator free equilibrium pointE3=(x^,0,z^) be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient conditions hold:(44)δ+ηγ+ωνz^>x^+ωsup.zνProof.
Consider the following positive definite real valued function:(45)L3=a1∫x^xS-x^Sds+a2y+a3∫z^zS-z^Sdswhere ai;i=1,2,3 are positive constants to be determined. Then the derivative of this function with respect to time can be written as(46)dL3dt=a1x-x^xdxdt+a2dydt+a3z-z^zdzdtHence by using system (3) with some algebraic manipulations we obtain that(47)dL3dt=-a11+βx+x^x-x^2-a1-a3νx-x^z-z^-a1-a2γxy-a2δ+η-a1x^y-a2μy2+a3ωyz-a3ωz^y-a3τz-z^2Therefore, by choosing the values of ai;i=1,2,3, so that a1=1, a2=1/γ, and a3=1/ν, and then using the supremum value of z, we obtain that(48)dL3dt≤-1+βx+x^x-x^2-μγy2-τνz-z^2-δ+ηγ-x^+ωνz^-ωνsup.zyClearly, condition (44) guarantees that dL3/dt<0. Again, as the function L3 is radially unbounded, then according to the Lyapunov first theorem E3 is a globally asymptotically stable point. Hence the proof is complete.Theorem 9.
Let the coexistence equilibrium point of system (3) be locally asymptotically stable in the interior of R+3. Then it is globally asymptotically stable provided that the following sufficient condition holds:(49)ω2ν≤min.μγ,τνProof.
Consider the following positive definite real valued function:(50)L4=d1∫x∗xS-x∗Sdx+d2∫y∗yS-y∗Sdydt+d3∫z∗zS-z∗Sdzdtwhere di;i=1,2,3 are positive constants to be determined. Then the derivative of this function with respect to time can be written as(51)dL4dt=d1x-x∗-x-x∗-y-y∗-z-z∗-βx2-x∗2+d2y-y∗γx-x∗-μy-y∗+d3z-z∗νx-x∗+ωy-y∗-τz-z∗Further simplification leads to(52)dL4dt=-d11+βx+x∗x-x∗2-d1-d2γx-x∗y-y∗-d1-d3νx-x∗z-z∗-d2μy-y∗2-d3τz-z∗2+d3ωy-y∗z-z∗Therefore, by choosing the values of di;i=1,2,3, so that d1=1, d2=1/γ, and d3=1/ν, and then by using some algebraic manipulations we obtain that(53)dL4dt≤-1+βx+x∗x-x∗2-ω2νy-y∗-ω2νz-z∗2-μγ-ω2νy-y∗2-τν-ω2νz-z∗2Clearly, condition (49) guarantees that dL4/dt is negative definite. Furthermore, because the function L4 is radially unbounded then according to the Lyapunov first theorem the coexistence equilibrium point E4 is a globally asymptotically stable point. Hence the proof is complete.


## 4. Bifurcation Analysis

In this section an investigation of a qualitative change in the dynamical behavior of system (3) under the effect of varying a specific parameter is carried out. As the necessary but not sufficient condition for the bifurcation to occur is the nonhyperbolic property of the equilibrium point, then the parameter is chosen so that the equilibrium point is a nonhyperbolic equilibrium point. An application to the Sotomayor’s bifurcation theorem is used throughout this section to investigate the occurrence and specify the type of bifurcation.Indeed the purpose of this section is to understand the effect and specify the parameters set that has a qualitative change on the dynamical behavior (bifurcation) of the system. In order to satisfy this purpose different parameters are selected so that the Jacobian matrix at an equilibrium point has zero eigenvalue and then study the effect of varying such parameter on the dynamical behavior of the system.Now for simplifying the notations rewrite system (3) in the vector form as follows: (54)dXdt=FX,withX=x,y,zTandF=F1,F2,F3T=xf1,yf2,zf3TAlso the second derivative of F with respect to X can be written as (55)D2FV,V=-21+3βxu12-2u1u2-2u1u32γu1u2-2μu222νu1u3+2ωu2u3-2τu32Here V=u1,u2,u3T is a general vector.Theorem 10.
Assume thatx-<(δ+η); then, as the parameter ν passes through the value ζ+ε/x-≡ν∗, system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occur.Proof.
Note that whenν=ζ+ε/x-≡ν∗, then the Jacobian matrix of system (3) at E1 can be written as(56)J1=JE1,ν∗=-x--2βx-2-x--x-0γx--δ+η0000So J1 has the following eigenvalues λ11∗=-x-+2βx-2<0, λ12∗=γx--δ+η<0, and λ13∗=0 and hence the necessary but not sufficient condition for bifurcation is satisfied and E1 is a nonhyperbolic point.
LetV1=v1,v2,v3T be the eigenvector of J1 corresponding to the eigenvalue λ13∗=0. Then straightforward computation gives that V1=γ1v3,0,v3T, where v3 represents any nonzero real number and γ1=-1/1+2βx-<0.
Also, letΨ1=ψ1,ψ2,ψ3T be the eigenvector of J1T that corresponds to the eigenvalue λ13∗=0. Then direct calculation shows that Ψ1=0,0,ψ3T, where ψ3 is any nonzero real number.
Because∂F/∂ν=Fν=0,0,xzT, hence we obtain that FνE1,ν∗=0,0,0T, which yields(57)Ψ1TFνE1,ν∗=0Thus system (3) at E1 with ν=ν∗ does not experience saddle-node bifurcation in view of Sotomayor theorem. Moreover, we have (58)Ψ1TDFνE1,ν∗V1=x-v3ψ3≠0where DFν represents the derivative of Fν with respect to X. Also by using Equation (55) at E1,ν∗ with the eigenvector V1 we obtain that(59)Ψ1TD2FE1,ν∗V1,V1=2ψ3v32ν∗γ1-2τNow as γ1<0 then ν∗γ1-2τ<0 and hence Ψ1TD2FE1,ν∗V1,V1≠0. Accordingly by Sotomayor theorem [23], system (3) near the equilibrium point E1 with ν=ν∗ undergoes a transcritical bifurcation but pitchfork cannot occur.Theorem 11.
Assume that(60)ζ+ε=νx~+ωy~Then near the scavenger free equilibrium point E2, system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occur.Proof.
Note that condition (60) is equivalent to the case when the parameter that stands for the death rate of scavenger ζ passes through the positive value ζ∗≡νx~+ωy~-ε. So the Jacobian matrix of system (3) at E2 when ζ=ζ∗ can be written as(61)J2=JE2,ζ∗=-x~-2βx~2-x~-x~γy~-μy~0000Clearly J2 has zero eigenvalues λ23∗=0 with two negative real parts eigenvalues given by Equation (18b).
LetV2=v1,v2,v3T be the eigenvectors of J2 corresponding to the eigenvalue λ23∗=0. So, direct computation shows that V2=γ1v3,γ2v3,v3T, where v3 represents any nonzero real number, while γ1=-μ/μ(1+2βx~)+γ<0 and γ2=-γ/μ1+2βx~+γ<0.
LetΨ2=ψ1,ψ2,ψ3T be the eigenvector of J2T corresponding to the zero eigenvalue λ23∗=0. Then straightforward calculation shows that Ψ2=0,0,ψ3T, where ψ3 is any nonzero real number.
Because∂F/∂ζ=Fζ=0,0,-zT, hence we obtain that FζE2,ζ∗=0,0,0T. Therefore (62)Ψ2TFζE2,ζ∗=0Thus system (3) at the scavenger free equilibrium point E2 with ζ=ζ∗ does not undergo saddle-node bifurcation in view of Sotomayor theorem.
Now, we have(63)Ψ2TDFζE2,ζ∗V2=-v3ψ3≠0where DFζ represents the derivative of Fζ with respect to X; moreover (64)Ψ2TD2FE2,ζ∗V2,V2=2ψ3v32νγ1+ωγ2-τwhere D2F is the second derivative of F with respect to X given in Equation (55). Clearly Ψ2TD2FE2,ζ∗V2,V2≠0 due to the sign of γi;i=1,2 and hence system (3) undergoes a transcritical bifurcation near E2 when ζ=ζ∗ while pitchfork cannot occur.Theorem 12.
Assume that(65)γx^=δ+ηThen, near the predator free equilibrium point E3, system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occur.Proof.
Observe that condition (65) is equivalent to the case when the predator’s death rate parameter δ passes through the positive value δ∗=γx^-η. Further the Jacobian matrix of system (3) at E3 with δ=δ∗ can be written as (66)J3=JE3,δ∗=-x^-2βx^2-x^-x^000νz^ωz^-τz^Therefore J3 has zero eigenvalues λ32∗=0 with two negative real parts eigenvalues given by Equation (20b).
LetV3=v1,v2,v3T be the eigenvector of J3 corresponding to the zero eigenvalue λ32∗=0. Then, direct computation shows that V3=δ1v2,v2,δ2v2T, where v2 represents any nonzero real number and δ1=-(τ+ω)/τ1+2βx^+ν<0 and δ2=ω1+2βx^-ν/τ1+2βx^+ν.
LetΨ3=ψ1,ψ2,ψ3T be the eigenvector of J3T corresponding to the eigenvalue λ32∗=0. Then simple calculation shows that Ψ3=0,ψ2,0T, where ψ2 is any nonzero real number.
Now due to∂F/∂δ=Fδ=0,-y,0T, hence FδE3,δ∗=0,0,0T. Therefore we obtain that (67)Ψ3TFδE3,δ∗=0Consequently system (3) at the predator free equilibrium point E3 with δ=δ∗ does not undergo saddle-node bifurcation in view of Sotomayor theorem.
Now, we have(68)Ψ3TDFδE3,δ∗V3=-v2ψ2≠0where DFδ represents the derivative of Fδ with respect to X; moreover (69)Ψ3TD2FE3,δ∗V3,V3=2ψ2v22γγ1-2μwhere D2F is the second derivative of F with respect to X given in Equation (55). Clearly Ψ2TD2FE3,δ∗V3,V3≠0 due to the sign of γ1 and hence system (3) undergoes a transcritical bifurcation near E3 when δ=δ∗ while pitchfork cannot occur.Finally, because the determinant ofJE4, which is given by A3 in Equation (22b), is always positive, then there is no possibility for the positive equilibrium point to be a nonhyperbolic point and hence there is no possibility for local bifurcation to occur. Moreover if we violate condition (25) so that the following condition holds:(70)a31a33+a21a32>0⇛γωy∗>ντz∗here aij;i,j=1,2,3 are given in Eq. (22a), then it is easy to verify that as the parameter ω passes through the positive value ω∗ given by(71)a32=Λ+a13a31a33-a13a21⇛ω∗=Λ+ντx∗z∗2γx∗y∗z∗where (72)Λ=-a11+a22a11a22-a12a21-a11a33a11+a22+a33-a22a33a22+a33-a11a22a33-a12a21then we obtain that ∆=0≈A1A2=A3, and hence the roots of the characteristic Equation (22b) can be written as λ1ω∗=-A1<0 and λ2,3ω∗=∓iA2. Therefore system (3) has two pure imaginary eigenvalues with the third one real and negative at (E4,ω∗). Consequently, there is possibility of occurrence of Hopf bifurcation if and only if the tranversilaty condition holds [23].As the general form of the complex eigenvalues in the neighborhood‎ ofω∗ can be written in the form λ2,3ω=σ1(ω)∓iσ2(ω), then by substituting one of them in Equation (22b) in case of A1A2=A3 and taking the derivative with respect to ω and then simplifying the resulting equation by comparing the real and imaginary parts and solving the resulting system for σ1′(ω), we obtain that(73)σ1′ω=-ΘωΨω+ΓωΦωΨω2+Φω2Here(74)Θω=A1′ωσ1ω2-A1′ωσ2ω2+A2′ωσ1ω+A3′ωΨω=3σ1ω2+2A1ωσ1ω-3σ2ω2+A2ωΓω=2A1′ωσ1ωσ2ω+A2′ωσ2ωΦω=6σ1ωσ2ω+2A1ωσ2ωRecall that, according to the Hopf bifurcation theorem, if the following transversality condition holds at ω=ω∗ then the system undergoes a Hopf bifurcation; otherwise there is no such type of bifurcation too:(75)σ1′ωω=ω∗≠0Clearly condition (75) is equivalent to ΘωΨω+Γ(ω)Φ(ω)≠0 at ω=ω∗. Straightforward computation shows that Θω∗Ψω∗+Γω∗Φω∗=0; hence Hopf bifurcation is not possible to occur around the E4 for system (3).


## 5. Numerical Simulation

In this section an investigation for the global dynamics of system (3) is carried out using numerical simulation. The objective is to understand the effects of varying the parameters of the system and confirm our obtained analytical results. For that purpose the following set of hypothetical parameters values is adopted: (76)α=0.2,β=0.3,γ=0.7,δ=0.1,η=0.1μ=0.1,ν=0.5,ω=0.2,ζ=0.1,ε=0.1,τ=0.1It is observed that system (3) approaches asymptotically the positive equilibrium point E4=(0.32,0.27,0.17) starting from three different sets of initial values as shown in phase portrait given in Figure 1 and their time series represented by Figure 2. This coincides with the analytical result obtained in Theorem 9, which determined the sufficient condition for global stability of positive equilibrium point and is given by condition (49), as the data (76) fails to satisfy condition (49) but the trajectory still approaches asymptotically the positive point.Figure 1
The phase portrait of system (3) using data (76).Figure 2
Time series of the attractor in Figure1. (a) Prey as a function of time. (b) Predator as a function of time. (c) Scavenger as a function of time.
(a) (b) (c)Further numerical simulations have been done for other sets of parameters values satisfying the sufficient condition for global stability of positive equilibrium point given by Equation (49); in all these cases the numerical findings coincide with analytical findings and the trajectories approach the positive equilibrium point.It is observed further that forα<0.45 the solution of system (3) still approaches asymptotically the positive equilibrium point and for 0.45≤α<0.69 the solution of system (3) approaches asymptotically scavenger free equilibrium point E2, while it approaches asymptotically the axial equilibrium point E1 when 0.69≤α<1 as shown in the typical figures given by Figures 3(a)-3(b). Otherwise the trajectory of system (3) approaches asymptotically the trivial equilibrium point for α>1; see Figure 3(c). Note that the parameter α exists in the form of x- and x~ as shown in Equation (9) and Equation (10b), respectively. Moreover the stability conditions of E1 and E2 given in (17) and (19) depend on x- and x~. Also the positivity of the rate of change of the prey species represented by dx/dt depends on condition (5) that leads to the collapse of the system totally. Therefore the parameter α plays a vital role in changing the dynamical behavior of system (3) as shown in Figure 3.Figure 3
Time series for the solution of system (3) using data (76) with typical values of α. (a) The system approaches scavenger free equilibrium point E2=(0.3,0.12,0) when α=0.55. (b) The system approaches axial equilibrium point E1=(0.19,0,0) when α=0.8. (c) The system approaches trivial equilibrium point E0=(0,0,0) when α=1.01.
(a) (b) (c)From the real world life, it is well known that when the ratio of harvesting ‎of prey to ‎intrinsic growth rate of prey increases then the prey population starts decreasing and ‎hence ‎all the predator species which feed on it will be extinct; however when this ratio ‎exceeds 1 ‎then the decay in the prey population will be greater than growth and hence the whole ‎system will collapse. ‎Accordingly, the numerical result, obtained in Figure3, agree with the real ‎world life in the ‎environment. Moreover this is coinciding with obtained analytical results given by Equation (15) that represents the existence of three negative eigenvalues around the trivial equilibrium point provided that α>1.‎It is observed that, forγ>0.43, system (3) approaches asymptotically the positive equilibrium point; however it approaches the predator free equilibrium point otherwise as shown in the typical figure given by Figure 4. This is biologically feasible due to the fact that when the ratio of prey’s attack rate by predator to the ‎intrinsic growth rate of prey increases then the predator will be growing; otherwise it will be ‎decreasing reaching extinction. In addition, this is coinciding with the stability condition of predator free equilibrium point given by Equation (21).Figure 4
Time series for the solution of system (3) using data (76) with γ=0.35. The system approaches E3=(0.45,0,0.28).On the other hand, for the parameter data (76) with δ≥0.23 or η≥0.23 the solution of system (3) approaches asymptotically the predator free equilibrium point as shown in the typical figure represented by Figure 5. However it approaches asymptotically the positive equilibrium point otherwise.Figure 5
Time series for the solution of system (3) using data (76) with δ=0.3. The system approaches E3=(0.45,0,0.28).Obviously the numerical result in Figure5 agree with the real life in the environment due to the fact that as the ratio of predator’s death rate to the intrinsic growth rate of the prey or ratio of harvesting of predator to intrinsic growth rate of prey species or both increases then the predator faces extinction; otherwise the system will persist. Again this is coinciding with the stability condition of predator free equilibrium point given by Equation (21).Now for the range0.33<ν<0.88 with the rest of parameters as given in data (76) the solution of system (3) still approaches the positive equilibrium point, while it approaches asymptotically scavenger free equilibrium point or predator free equilibrium point for ν≤0.33 or ν≥0.88, respectively, as seen in the following typical figures given by Figure 6.Figure 6
Time series for the solution of system (3) using data (76) with typical values of ν. (a) The system approaches E3=(0.26,0,0.51) for ν=0.95. (b) The system approaches E2=(0.34,0.42,0) for ν=0.25.
(a) (b)Note that as the ratio of prey’s attack rate by scavenger to the intrinsic growth rate of prey decreases, the availability of main source of food for scavenger will be decreasing too leading to extinction in scavenger in the environment, while if this ratio increases then consuming of prey population becomes more and this is causing lack in availability of food for predator due to competition with scavenger and hence the predator will face extinction this time. This is exactly shown in our obtained numerical results which in turn coincide with our analytical results given by existence of scavenger free equilibrium point represented by Equations (11a) and (11b) and stability condition of predator free equilibrium point represented by Equation (19).Similarly it is observed that, for the data (76) with ω>0.06, system (3) still approaches asymptotically the positive equilibrium point; however it approaches scavenger free equilibrium point when ω≤0.06 as shown in Figure 7.Figure 7
Time series for the solution of system (3) using data (76) with ω=0.02. The system approaches E2=(0.34,0.42,0).Moreover, it is observed that for the parameter data (76) with ζ≥0.16 or ε≥0.16 the solution of system (3) approaches asymptotically the scavenger free equilibrium point as shown in the typical figure represented by Figure 8. However it approaches asymptotically the positive equilibrium point otherwise.Figure 8
Time series for the solution of system (3) using data (76) with ζ=0.2. The system approaches E2=(0.34,0.42,0).Note that the numerical result in Figure8 agrees with the real life in the environment due to the fact that as the ratio of scavenger death rate to the intrinsic growth rate of the prey or the ratio of harvesting of scavenger to intrinsic growth rate of prey species or both increases then the scavenger facing extinction; otherwise the system will persist. Again this is coinciding with the stability condition of scavenger free equilibrium point given by Equation (19).Finally the effect of varying toxicant parameters (β,μ,τ) on the dynamics of system (3) is also investigated and it is obtained that increasing these parameters or some of them causes extinction in some of the system species as shown in the typical figures given by Figures 9 and 10.Figure 9
Time series for the solution of system (3) using data (76) with varying some parameters. (a) The system approaches positive equilibrium point E4=(0.31,0.22,0.03) for α=0.4, β=0.3. (b) The system approaches scavenger free equilibrium point E2=(0.31,0.19,0) for =0.4,β=0.9.
(a) (b)Figure 10
Time series for the solution of system (3) using data (76) with varying some parameters. (a) The system approaches positive equilibrium point E4=(0.43,0.02,0.19) for α=0.3, β=0.3, η=0.2. (b) The system approaches predator free equilibrium point E2=(0.42,0,0.11) for α=0.3, β=0.9, η=0.2.
(a) (b)The effect of toxic substance on the existence of species in the environment is well known; Figures9 and 10 show clearly the same effect, as the ‎increasing of toxic rates in the system causes extinction in the predator or scavenger populations. ‎According to the obtained numerical results above, which describe the global dynamics of system (3), for the set of parameters values given by (76), under the varying of specific parameter each time, the effects of varying each parameter on the dynamical behavior of the system can be summarized as follows:(1)
Increasing the ratio of harvesting of prey to intrinsic growth rate of prey species in the range0.45≤α<0.69 leads to losing the persistence of the system and the solution approaches scavenger free equilibrium point. However increasing this parameter further in the range 0.69≤α<1 cases extinction in predator species too and the solution approaches axial equilibrium point. Finally, violating the condition (5) leads to the collapse of the system and the solution approaches trivial equilibrium point.(2)
Decreasing the ratio of prey’s attack from predator to intrinsic growth rate of prey in the rangeγ≤0.43 leads to extinction in predator species and the system approaches asymptotically predator free equilibrium point.(3)
Increasing the ratio of predator’s death rate to the intrinsic growth rate of the prey in the rangeδ≥0.23 or the ratio of harvesting of predator to intrinsic growth rate of prey species η≥0.23 leads to extinction in predator species too and the system approaches asymptotically predator free equilibrium point.(4)
Increasing the ratio of prey’s attack rate from scavenger to intrinsic growth rate of prey in the rangeν≥0.88 cases extinction in predator species and the system approaches predator free equilibrium point, while decreasing it in the range ν≤0.33 leads to extinction in the scavenger species and the system approaches asymptotically scavenger free equilibrium point. Otherwise the system still persist at the positive equilibrium point.(5)
Decreasing the ratio of predator’s attack rate from scavenger to the prey’s attack rate from predator in the rangeω≤0.06 leads to extinction in scavenger species and the system approaches asymptotically scavenger free equilibrium point.(6)
Increasing the ratio of scavenger death rate to the intrinsic growth rate of the prey in the rangeζ≥0.16 or the ratio of harvesting of scavenger to intrinsic growth rate of prey species ε≥0.16 leads to extinction in scavenger species and the system approaches asymptotically scavenger free equilibrium point.(7)
Finally, increasing the toxic rates above a specific value leads to extinction in of either predator or scavenger species and the system approaches asymptotically one of the boundary equilibrium points.Finally, according to the above numerical results, the dynamical behavior of system (3) appears to be sensitive to varying in the parameters values and the bifurcation may occur at specific values of parameter so that the stability of the system changed from positive equilibrium point to other boundary points which means losing the persistence of the system. The intensity of sensitivity of system to varying in the parameters values is different from one parameter to another. In fact it is observed that system (3) has three bifurcation points at the parameter α, while it has two bifurcation points at the parameter ν. This makes these parameters (α and ν) more important than others in the system due to their effects on the dynamical behavior of the system.


## 6. Discussion and Conclusions

In this paper, a food web model consisting of prey-predator-scavenger is proposed and analyzed. It is assumed that two types of biological factors are considered in this system; these are represented by linear type of harvest and nonlinear type of environment toxicant. The investigation of the dynamical behavior of this system shows that the system has at most five equilibrium points and their solution exists and is unique and bounded. The stability conditions of all these possible equilibrium points are established. The conditions that make the system persist are obtained. Suitable Lyapunov functions are used to investigate the global dynamics of each equilibrium point. Moreover, an investigation of the occurrence of local bifurcation in view of Sotomayor bifurcation theorem is carried out. It is observed that system (3) undergoes transcritical bifurcation; neither saddle node nor pitchfork bifurcation can occur near the boundary equilibrium points when a specific parameter passes through a critical value. However, the system does not approach a periodic dynamics in the positive octant or in the boundary planes; rather, it approaches asymptotically one of their possible equilibrium points. Finally, the global dynamics of system (3) is also investigated numerically using the hypothetical set of parameters given by Equation (76). The objective is to understand the effects of varying the system parameters on the global dynamics of the system in addition to confirming our obtained analytical results. The obtained numerical results are represented in the form of some time series figures that show the behavior of the solution clearly.This type of study is important due to the availability of different types of species in the environment which contains at the same time different factors, such as harvesting and toxicant substance and others, affecting the persistence of all species. Therefore, mathematical modeling is used to understand the dynamics of such food web model and the effects of varying the parameters values that represent different factors in the environment.In order to understand the dynamics behavior of such a food web model the stability analysis with the help of different mathematical tools, such as linearization, Lyapunov functions, and persistence, is used to investigate the dynamical behavior of the proposed model and establish the persistence condition. However bifurcation analysis is used to understand the effects of varying the system parameters. Finally numerical simulation is also included to complete the analysis and verify the obtained results.Keeping the above in view, it is observed that the proposed food web model in this paper is very sensitive to varying in the parameters values and undergoes a transcritical bifurcation at different equilibrium points.



---

*Source: 1573516-2019-03-26.xml*