International Journal of Mathematics and Mathematical Sciences
Volume 2019 (2019), Article ID 1573516, 17 pages
https://doi.org/10.1155/2019/1573516

Research Article

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

Huda Abdul Satar and Raid Kamel Naji

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

Correspondence should be addressed to Raid Kamel Naji; rknaji@gmail.com

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 region

Here 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 that

while

The predator free equilibrium point located in the interior of positive quadrant of -plane and denoted by exists uniquely provided that

while

Finally, 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 as

Therefore the eigenvalues are given by

Hence, 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 by

Clearly, one of the eigenvalues of is given by

where 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 by

Consequently, the eigenvalues of are given by

where 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 by

where

with

Consequently, 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:

and

Obviously, 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 have

Clearly 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 obtain

where are given in the right-hand side of system (3); therefore

Now, 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 have

Obviously, 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 as

Hence by using system (3) with some algebraic manipulations we obtain that

Therefore, by choosing the values of , so that , , and , we get

Consequently, 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 get

Now by using the supremum value of we obtain that

Clearly, 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 as

Hence by using system (3) with some algebraic manipulations we obtain that

Therefore, by choosing the values of , so that , , and , and then using the supremum value of , we obtain that

Clearly, 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 as

Further simplification leads to

Therefore, by choosing the values of , so that , , and , and then by using some algebraic manipulations we obtain that

Clearly, 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 as

So 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 yields

Thus 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 that

Now 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 that

Then 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 as

Clearly 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 by

where

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 that

Here

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:

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.

Figure 1: The phase portrait of system (3) using data (76).

(a)

(b)

(c)

Figure 2: Time series of the attractor in Figure 1. (a) Prey as a function of time. (b) Predator as a function of time. (c) Scavenger as a function of time.

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 the solution of system (3) still approaches asymptotically the positive equilibrium point and for the solution of system (3) approaches asymptotically scavenger free equilibrium point , while it approaches asymptotically the axial equilibrium point when 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 ; see Figure 3(c). Note that the parameter exists in the form of and as shown in Equation (9) and Equation (10b), respectively. Moreover the stability conditions of and given in (17) and (19) depend on and . Also the positivity of the rate of change of the prey species represented by 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.

(a)

(b)

(c)

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 when . (b) The system approaches axial equilibrium point when . (c) The system approaches trivial equilibrium point when .

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 Figure 3, 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 .

‎It is observed that, for , 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 . The system approaches .

On the other hand, for the parameter data (76) with or 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 . The system approaches .

Obviously the numerical result in Figure 5 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 range 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 or , respectively, as seen in the following typical figures given by Figure 6.

(a)

(b)

Figure 6: Time series for the solution of system (3) using data (76) with typical values of . (a) The system approaches for . (b) The system approaches for .

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 , system (3) still approaches asymptotically the positive equilibrium point; however it approaches scavenger free equilibrium point when as shown in Figure 7.

Figure 7: Time series for the solution of system (3) using data (76) with . The system approaches .

Moreover, it is observed that for the parameter data (76) with or 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 . The system approaches .

Note that the numerical result in Figure 8 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.

(a)

(b)

Figure 9: Time series for the solution of system (3) using data (76) with varying some parameters. (a) The system approaches positive equilibrium point for , . (b) The system approaches scavenger free equilibrium point for .

(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 for , , . (b) The system approaches predator free equilibrium point for , , .

The effect of toxic substance on the existence of species in the environment is well known; Figures 9 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 range leads to losing the persistence of the system and the solution approaches scavenger free equilibrium point. However increasing this parameter further in the range 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 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 or the ratio of harvesting of predator to intrinsic growth rate of prey species 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 cases extinction in predator species and the system approaches predator free equilibrium point, while decreasing it in the range 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 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 or the ratio of harvesting of scavenger to intrinsic growth rate of prey species 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.

Data Availability

No data were used to support this study.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors are very much thankful to the reviewers for their valuable suggestions and constructive comments in the first and second rounds of review. Their useful comments have contributed to the improvement of the authors’ work.

References

R. Arditi and L. R. Ginzburg, “Coupling in predator-prey dynamics: ratio-dependence,” Journal of Theoretical Biology, vol. 139, no. 3, pp. 311–326, 1989.

Y. Kuang and E. Beretta, “Global qualitative analysis of a ratio-dependent predator-prey system,” Journal of Mathematical Biology, vol. 36, no. 4, pp. 389–406, 1998.

K. Saleh, “A ratio-dependent predator-prey system with quadratic predator harvesting,” Asian Transactions on Basic & Applied Sciences, vol. 02, no. 4, pp. 21–25, 2013.

R. P. Gupta and P. Chandra, “Bifurcation analysis of modified Leslie-Gower predator-prey model with Michaelis-Menten type prey harvesting,” Journal of Mathematical Analysis and Applications, vol. 398, no. 1, pp. 278–295, 2013.

P. Feng, “Analysis of a delayed predator-prey model with ratio-dependent functional response and quadratic harvesting,” Applied Mathematics and Computation, vol. 44, no. 1-2, pp. 251–262, 2014.

S. Gakkhar and R. K. Naji, “On a food web consisting of a specialist and a generalist predator,” Journal of Biological Systems, vol. 11, no. 4, pp. 365–376, 2003.

T. K. Kar and K. S. Chaudhuri, “Harvesting in a two-prey one-predator fishery: a bioeconomic model,” The ANZIAM Journal, vol. 45, no. 3, pp. 443–456, 2004.

K. Tanabe and T. Namba, “Omnivory creates chaos in simple food web models,” Ecology, vol. 86, no. 12, pp. 3411–3414, 2005.

R. K. Naji and I. H. Kasim, “The dynamics of food web model with defensive switching property,” Nonlinear Analysis: Modelling and Control, vol. 13, no. 2, pp. 225–240, 2008.

A. Al Basheer, R. D. Parshad, E. Quansah, S. Yu, and R. K. Upadhyay, “Exploring the dynamics of a Holling-Tanner model with cannibalism in both predator and prey population,” International Journal of Biomathematics, vol. 11, no. 1, Article ID 1850010, 29 pages, 2018.

J. P. Previte and K. A. Hoffman, “Period doubling cascades in a predator-prey model with a scavenger,” SIAM Review, vol. 55, no. 3, pp. 523–546, 2013.

R. P. Gupta and P. Chandra, “Dynamical properties of a prey-predator-scavenger model with quadratic harvesting,” Communications in Nonlinear Science and Numerical Simulation, vol. 49, pp. 202–214, 2017.

J. E. Jansen and R. A. Van Gorder, “Dynamics from a predator-prey-quarry-resource-scavenger model,” Theoretical Ecology, vol. 11, no. 1, pp. 19–38, 2018.

H. Abdul Satar and R. K. Naji, “Stability and bifurcation in a prey–predator–scavenger system with michaelis–menten type of harvesting function,” Differential Equations and Dynamical Systems, pp. 1–24, 2019.

Y. F. Lv, R. Yuan, and Y. Z. Pei, “A prey-predator model with harvesting for fishery resource with reserve area,” Applied Mathematical Modelling, vol. 37, no. 5, pp. 3048–3062, 2013.

M. Sen, P. D. Srinivasu, and M. Banerjee, “Global dynamics of an additional food provided predator-prey system with constant harvest in predators,” Applied Mathematics and Computation, vol. 250, pp. 193–211, 2015.

K. Belkhodja, A. Moussaoui, and M. A. Aziz Alaoui, “Optimal harvesting and stability for a prey-predator model,” Nonlinear Analysis: Real World Applications, vol. 39, pp. 321–336, 2018.

T. K. Kar and K. S. Chaudhuri, “On non-selective harvesting of two competing fish species in the presence of toxicity,” Ecological Modelling, vol. 161, no. 1-2, pp. 125–137, 2003.

S. Sinha, O. P. Misra, and J. Dhar, “Study of a prey-predator dynamics under the simultaneous effect of toxicant and disease,” Journal of Nonlinear Sciences and Applications. JNSA, vol. 1, no. 2, pp. 102–117, 2008.

T. Das, R. N. Mukherjee, and K. S. Chaudhuri, “Harvesting of a prey-predator fishery in the presence of toxicity,” Applied Mathematical Modelling, vol. 33, no. 5, pp. 2282–2292, 2009.

R. K. Naji and A. N. Mustafa, “The dynamics of a prey- predator model with the existence of disease and pollution,” Journal of Mathematical and Computational Science, vol. 3, no. 1, pp. 94–123, 2013.

Q. Huang, H. Wang, and M. A. Lewis, “The impact of environmental toxins on predator-prey dynamics,” Journal of Theoretical Biology, vol. 378, pp. 12–30, 2015.

L. Perko, Differential Equations and Dynamical Systems, Springer, New York, NY, USA, 3rd edition, 2001.