International Journal of Mathematics and Mathematical Sciences

Volume 2019, Article ID 1573516, 17 pages

https://doi.org/10.1155/2019/1573516

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

Department of Mathematics, College of Science, University of Baghdad, Baghdad, Iraq

Correspondence should be addressed to Raid Kamel Naji; moc.liamg@ijankr

Received 28 October 2018; Revised 14 January 2019; Accepted 5 March 2019; Published 26 March 2019

Academic Editor: Rodica D. Costin

Copyright © 2019 Huda Abdul Satar and Raid Kamel Naji. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 rate and carrying capacity .(2)All the species are assumed to be subjected to a combined harvesting effort , with positive catchability coefficients for the prey, predator, and scavenger respectively.(3)Lotka-Volterra types of functional responses are used for describing the predation and scavenge process with as positive maximum attack rates, while the positive constants and 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 rates and , 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 represent the coefficients of toxicity of prey, predator, and scavenger, respectively, with 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, when , , and represent the density of prey, predator, and scavenger at time , the dynamics of the prey-predator-scavenger model can be described using the following set of differential equations:with , , and . 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 () is used in the first equation while quadratic term ( and ) 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:Here 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:Therefore, system (3) has the following domain:Clearly, 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 Otherwise 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 are uniformly bounded in the regionHere and are given in the proof.*

*Proof. *From the first equation we have ; then it is easy to verify that for we get . Now let ; then the derivative of with respect to time can be written as Here . Then, according to the above differential inequality, it is easy to verify that for ; we obtain .

Now assume that ; then the derivative of with respect to time can be written as Here . Then direct computation shows that, for , we obtain . 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 . 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 point and the axial equilibrium point always exist; here The scavenger free equilibrium point located in the interior of positive quadrant of -plane and denoted by exists uniquely provided thatwhileThe predator free equilibrium point located in the interior of positive quadrant of -plane and denoted by exists uniquely provided thatwhileFinally, the coexistence or positive equilibrium point is located at the interior of positive octant and denoted by , where

Here is a positive root of the following quadratic polynomial equation:

where

Clearly, exists uniquely provided that the following conditions hold: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 has the following eigenvalues: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 point can be written asTherefore the eigenvalues are given byHence, the axial equilibrium point is locally asymptotically stable if the following condition holds:Now the Jacobian matrix evaluated at the scavenger free equilibrium point is given byClearly, one of the eigenvalues of is given bywhere and . Therefore these eigenvalues have negative real parts and hence is locally asymptotically stable if and only if the following condition holds:Further, the Jacobian matrix evaluated at the predator free equilibrium point is given byConsequently, the eigenvalues of are given bywhere and . Hence all these eigenvalues have negative real parts and hence is locally asymptotically stable if and only if the following condition holds:Finally, the Jacobian matrix evaluated at the coexistence equilibrium point is given by The characteristic equation associated with is given bywherewithConsequently, 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 provided that the following sufficient condition holds:*

*Proof. *According to the Routh-Hurwitz criterion Equation (22b) has three roots with negative real parts provided that ; and . Straightforward computation shows that Clearly the given sufficient condition (25) leads to . Hence all the eigenvalues of 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 (), 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 plane while the second subsystem is obtained in the absence of predator and located in the positive quadrant of plane. These subsystems can be written as follows:andObviously, subsystem (27) has a unique positive equilibrium point which belongs to interior of positive quadrant of plane given by . However, subsystem (28) has a unique positive equilibrium point which belongs to interior of positive quadrant of plane given by .

Theorem 3. *There is no periodic dynamics in the interior of positive quadrant of plane surrounding a unique equilibrium point .*

*Proof. *Define that obviously satisfies and continuously differentiable function in the interior of positive quadrant of plane. Now, we haveClearly is not equal to zero and it is not change sign in the interior of positive quadrant of plane. So according to the Dulac-Bendixson’s criterion there is no closed curve surrounding a unique equilibrium point .

Now according to the above theorem along with the Poincare-Bendixson theorem, the equilibrium point is a globally asymptotically stable in the interior of positive quadrant of the 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 of plane surrounding a unique equilibrium point .*

Again we can easily conclude that is globally asymptotically stable in the interior of positive quadrant of the 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 provided that the following conditions hold:*

*Proof. *Consider the function , where are undetermined positive constants. Obviously is a positive function defined in the interior of .

Consequently we obtainwhere are given in the right-hand side of system (3); thereforeNow, since there are no periodic attractors in the boundary planes, then the only possible omega limit sets of system (3) are the equilibrium points , , . Thus the system is uniformly persistent if we can proof that at each of these points. Also we haveObviously, for all values of is sufficiently large enough with respect to

for any positive constants provided that the condition (30a) is satisfied. However, 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 point be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient condition holds:*

*Proof. *Consider the following positive definite real valued function:where are positive constants to be determined. Then the derivative of this function with respect to time can be written asHence by using system (3) with some algebraic manipulations we obtain thatTherefore, by choosing the values of , so that , , and , we getConsequently, is negative definite under the local stability condition (17) and condition (34). Moreover it is clear that the function is radially unbounded; then according to the Lyapunov first theorem is a globally asymptotically stable point. Hence the proof is complete.

Theorem 7. *Let the scavenger free equilibrium point be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient condition holds:*

*Proof. *Consider the following positive definite real valued function:Therefore, by choosing the values of , so that , , and , we getNow by using the supremum value of we obtain thatClearly, condition (39) guarantees that . Moreover as the function is radially unbounded then according to the Lyapunov first theorem is a globally asymptotically stable point. Hence the proof is complete.

Theorem 8. *Let the predator free equilibrium point be locally asymptotically stable. Then it is globally asymptotically stable provided that the following sufficient conditions hold:*

*Proof. *Consider the following positive definite real valued function:where are positive constants to be determined. Then the derivative of this function with respect to time can be written asHence by using system (3) with some algebraic manipulations we obtain thatTherefore, by choosing the values of , so that , , and , and then using the supremum value of , we obtain thatClearly, condition (44) guarantees that . Again, as the function is radially unbounded, then according to the Lyapunov first theorem 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 . Then it is globally asymptotically stable provided that the following sufficient condition holds:*

*Proof. *Consider the following positive definite real valued function:where are positive constants to be determined. Then the derivative of this function with respect to time can be written asFurther simplification leads toTherefore, by choosing the values of , so that , , and , and then by using some algebraic manipulations we obtain thatClearly, condition (49) guarantees that is negative definite. Furthermore, because the function is radially unbounded then according to the Lyapunov first theorem the coexistence equilibrium point 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: Also the second derivative of with respect to can be written as Here is a general vector.

Theorem 10. *Assume that ; then, as the parameter passes through the value , system (3) undergoes a transcritical bifurcation but neither saddle node nor pitchfork bifurcation can occur.*

*Proof. *Note that when , then the Jacobian matrix of system (3) at can be written asSo has the following eigenvalues , , and and hence the necessary but not sufficient condition for bifurcation is satisfied and is a nonhyperbolic point.

Let be the eigenvector of corresponding to the eigenvalue . Then straightforward computation gives that , where represents any nonzero real number and .

Also, let be the eigenvector of that corresponds to the eigenvalue . Then direct calculation shows that , where is any nonzero real number.

Because , hence we obtain that , which yieldsThus system (3) at with does not experience saddle-node bifurcation in view of Sotomayor theorem. Moreover, we have where represents the derivative of with respect to . Also by using Equation (55) at with the eigenvector we obtain thatNow as then and hence . Accordingly by Sotomayor theorem [23], system (3) near the equilibrium point with undergoes a transcritical bifurcation but pitchfork cannot occur.

Theorem 11. *Assume thatThen near the scavenger free equilibrium point , 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 . So the Jacobian matrix of system (3) at when can be written asClearly has zero eigenvalues with two negative real parts eigenvalues given by Equation (18b).

Let be the eigenvectors of corresponding to the eigenvalue . So, direct computation shows that , where represents any nonzero real number, while and .

Let be the eigenvector of corresponding to the zero eigenvalue . Then straightforward calculation shows that , where is any nonzero real number.

Because , hence we obtain that . Therefore Thus system (3) at the scavenger free equilibrium point with does not undergo saddle-node bifurcation in view of Sotomayor theorem.

Now, we have where represents the derivative of with respect to ; moreover where is the second derivative of with respect to given in Equation (55). Clearly due to the sign of and hence system (3) undergoes a transcritical bifurcation near when while pitchfork cannot occur.

Theorem 12. *Assume that Then, near the predator free equilibrium point , 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 . Further the Jacobian matrix of system (3) at with can be written as Therefore has zero eigenvalues with two negative real parts eigenvalues given by Equation (20b).

Let be the eigenvector of corresponding to the zero eigenvalue . Then, direct computation shows that , where represents any nonzero real number and and .

Let be the eigenvector of corresponding to the eigenvalue . Then simple calculation shows that , where is any nonzero real number.

Now due to , hence . Therefore we obtain that Consequently system (3) at the predator free equilibrium point with does not undergo saddle-node bifurcation in view of Sotomayor theorem.

Now, we have where represents the derivative of with respect to ; moreover where is the second derivative of with respect to given in Equation (55). Clearly due to the sign of and hence system (3) undergoes a transcritical bifurcation near when while pitchfork cannot occur.

Finally, because the determinant of , which is given by 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:here are given in Eq. (22a), then it is easy to verify that as the parameter passes through the positive value given bywhere then we obtain that , and hence the roots of the characteristic Equation (22b) can be written as and . Therefore system (3) has two pure imaginary eigenvalues with the third one real and negative at . 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 , then by substituting one of them in Equation (22b) in case of 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 , we obtain thatHereRecall 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:Clearly condition (75) is equivalent to at . Straightforward computation shows that ; hence Hopf bifurcation is not possible to occur around the 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: It is observed that system (3) approaches asymptotically the positive equilibrium point 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.