The Dynamics of a Food Web System: Role of a Prey Refuge Depending on Both Species

This paper aims to study the role of a prey refuge that depends on both prey and predator species on the dynamics of a food web model. It is assumed that the food transfer among the web levels occurs according to Lotka-Volterra functional response. The solution properties, such as existence, uniqueness, and uniform boundedness, are discussed. The local, as well as the global, stabilities of the solution of the system are investigated. The persistence of the system is studied with the assistance of average Lyapunov function. The local bifurcation conditions that may occur near the equilibrium points are established. Finally, numerical simulation is used to confirm our obtained results. It is observed that the system has only one type of attractors that is a stable point, while periodic dynamics do not exist even on the boundary planes.


Introduction
Food webs are established in ecosystems and are significant for all living organisms. Energy needs for all activities of animals are provided through food consumption. Furthermore, food webs represent the pathways by which energy flows from a level to another. They are fundamental and inescapable in any attempt to describe how the real world life in the environment is organized or complexes of species interact. The three types of food web samples are the bases of large measure ecosystems. Therefore, to understand the dynamical behavior of an ecosystem, it is substantial to realize the dynamics of the three kinds of food web models. Hence, attention has been given to the

ISSN: 0067-2904
Jaber and Bahlool Iraqi Journal of Science, 2021, Vol. 62, No. 2, pp: 639-657 640 dynamics with 3-species food chains and food web systems [1][2][3][4][5][6][7][8]. Later on, many researchers investigated the dynamical behavior of food webs with various types of functional responses and including many biological factors. Agarwal and Kumarb [9] proposed and studied the alternative resource effects for top predators using Holling type III model. They showed that an alternative resource has the ability to prohibit top predator annihilation. Xiao et al. [10] proposed and investigated a predator-prey model that associates time delay with a Holling type II and Allee impact in prey. They observed that if the Allee impact is high or the birth rate is low, then both species of the system are fading and hence the stability of the system is affected.
In fact, the refuge mechanism has a role in promoting the growth rate of prey and reducing both predator growth and the products lurking behavior of prey [11]. Therefore, the existence of prey refuges has substantial effects on the cohabitation of predators and their prey. There are three types of refuges; the first type is the supply of continual spatial protection for a small branch of the prey population, while the second type is the provision of temporary spatial preservation and, finally, the third type is the supply of a temporary refuge in numbers, which means the reduction of the venture of predation by rising the abundance of vulnerable prey [12]. The problem of prey-predator interactions under a prey refuge has been studied by some authors. Mukherjee [11] investigated a resource according to Holling type I and II functional responses and then showed the effect of constant prey refuge. Das et al. [13] described constant prey refuge and harvesting to both predator and prey. Ghosh et al. [14] studied the effects of extra food for predator with prey refuge. Santra et al. [15] investigated the dynamical actions with Crowley-Martin functional response-associated prey refuge. Finally, Molla et al. [16] suggested a model for Holling type-II prey-predator interactions in a prey refuge. In this paper, a Lotka-Volttera food web model with a prey refuge that consists of prey and predator is proposed and studied. The remaining sections are organized as follows. In section (2), the mathematical model with symbols of system is introduced. Section (3) deals with the analysis of local stability. In section (4), global stability analysis is investigated. Section (5) describes the persistence of the model. In section (6), the local bifurcation of the model is studied. The numerical simulation is given in section (7), while section (8) includes some discussion of the obtained results.

2.
The model formulation In this section, a mathematical modeling approach was used to study the role of a non-constant prey refuge in the dynamical behavior of a three-species food web system. It is assumed that the food web is consisting of prey at a lower level of density at time given by . The prey grows logistically in the absence of the predators and has a non-constant refuge property that is dependent on the density of their predator. The middle predator, denoted for density at time by , feeds on prey at the lower level and is preyed on by the predator at the higher level. It has another limited source of food in the absence of their preferred food. The top predator, denoted for density at time by , at the higher level, feeds on both the prey at the lower level and the middle predator at the second level. Accordingly, the dynamics of such a food web system can be described by using mathematical modeling with differential equations.
with and the model (1) becomes with . The description of the parameters for the system (2) is shown in Table-1. It is assumed that all the above parameters are positive constants and hence the domain of system (2) will be given by . Now to simplify the model analysis, the following parameters and dimensionless variables are used in system (2). (3) According to Eq. (3), the dimensionless system corresponding to system (2) will be where .The functions on the right hand side of system (4) are continuous and have continuous partial derivatives. Hence, they are Lipschitzian functions [17]. Therefore, the solution of system (1) exists and is unique. In addition, the following theorem shows the boundedness of the solution in . Theorem 1. The region is a uniformly bounded region that contains all the solutions of system (4). Proof. From system (4), we have that . By solving the inequality, it is simple to verify that as . Let , then after some algebraic manipulation we have . Then , where . By using Gronwall inequality [18], it is easy to obtain that as . Hence, all solutions are uniformly bounded in .

Local Stability Analysis
There are seven equilibrium points of system (4). Their existence conditions and local stability analyses are established.  The vanishing equilibrium point that is denoted by .  The first axial equilibrium point (FAEP) that is denoted by ( ̃ ) and always exists.  The second axial equilibrium point (SAEP) that is denoted by Clearly, the SAEP exists uniquely on axis if and only if it satisfies the following condition .
(5b)  The prey-free equilibrium point (PFEP) that is denoted by Obviously, the PFEP exists uniquely in the plane if and only if the following condition holds .
(6b)  The middle predator-free equilibrium point (MPFEP) that is denoted by is a root of the following equation and . Therefore, with the aid of the Descartes rule of signs, Eq. (7b) has either one or three positive roots but there are no negative roots, provided that the following condition holds .
(7c) By using Wolfram Mathematica 11.3, it is observed that Eq. (7b) has only one positive root denoted by ̅ and two other complex conjugate roots. Accordingly, the MPFEP remains uniquely in the plane, which stipulates that, in the extension to condition (7c), the following condition holds too ̅ ̅ .
(7d)  The top predator-free equilibrium point (TPFEP) that is denoted by Again, by utilizing the Descartes rule of signs, Eq. (8b) has either one or three positive roots but there is no negative roots, provided the following condition holds .
(8c) Similarly, Wolfram Mathematica 11.3 is used to compute the roots of Eq. (8b). It was found that Eq. (8b) has only one positive root denoted by ̿ and two other complex conjugate roots. Therefore, the TPFEP remains uniquely in the plane, which stipulates that, in the extension to condition (8c), the following condition holds too where is a positive intersection point for the following two isoclines Then, as we obtain that Clearly, each of Eq. (9d) and Eq. (9e) has a unique positive root on the axis, denoted by and respectively, provided that the next conditions hold Therefore, the two isoclines given by equations (9b) and (9c) have a unique positive intersection point in the interior of plane and then the PEP exists uniquely, provided that , , (10c) .
(10d) Now, to establish the local stability, the Jacobian matrix, that is denoted by of the system (4) at the is determined by It is clear that system (4) has the Jacobian at specified by The eigenvalues of are , , .
(12b) Therefore, is a saddle point. The Jacobian matrix at FAEP is Therefore, the eigenvalues of are given by , , .
(13b) Clearly, all the above eigenvalues are negative and is locally asymptotically stable in if the next conditions are satisfied , (13c) .
(13d) The Jacobian matrix at SAEP is Therefore, the eigenvalues can be written as The eigenvalues of will be negative and then is locally asymptotically stable if the next conditions are held The characteristic equation of is as follows Therefore, the eigenvalues of are Hence, if the next condition holds, will be locally asymptotically stable The characteristic of (16a) is shown as , Accordingly, all the above eigenvalues have real parts of less than zero and then is locally asymptotically stable if the next conditions are satisfied The Jacobian matrix at TPFEP is written as follows Therefore, the characteristic equation of is given by , Hence, the eigenvalues of are Consequently, all the above eigenvalues have real parts of less than zero and then is locally asymptotically stable if the next conditions are satisfied Next, the local stability conditions of the PEP are determined. Theorem 2: Suppose that the PEP of the system (4) exists, then it is a locally asymptotically stable with the following sufficient conditions are satisfied , , Proof: The Jacobian matrix at the PEP is as follows The characteristic equation of (19a) is , By using Routh-Hurwitz criterion [18], all the roots of equation (19b) have real parts of less than zero and then is locally asymptotically stable if , and are positive. It is found by computation that all Routh-Hurwitz constraints are satisfied under the sufficient conditions given by (18a) -(18d).

4.
Global tability analysis In t is section, the global dynamics of system (4) is investigated using Lyapnuov functions, as demonstrated in the next theorems [18]. Theorem 3: Suppose that the FAEP is locally asymptotically stable. T en is a globally asymptotically stable with the following sufficient conditions are satisfied ̃ , Using the above sufficient conditions (21a)-(21b) along with the local stability conditions, the derivative is a negative definite function. Hence, has a basin of attraction that satisfies the given sufficient conditions. The 5: Assume that the PFEP of system (4) is locally asymptotically stable. Then has a basin of attraction that satisfies the following sufficient conditions ̆ , Then, by using the above sufficient conditions (22a)-(22c), the derivative is a negative definite function. Hence, has a basin of attraction that satisfies the given sufficient conditions. Theorem 6: Suppose that the MPFEP is locally asymptotically stable. Then is a globally asymptotically stable, provided that Accordingly, by using the condition (23) along with the local stability condition (16f), the derivative becomes Clearly, is negative definite and hence is a globally asymptotically stable.

Theorem 7:
Assume that the TPFEP of system (4) is locally asymptotically stable. Then is globally asymptotically stable, provided that Proof: Define the following real valued function where the constants are greater than zero and must be determined. Obviously, is a positive definite function, which is defined for all and . Therefore, can be written as follows ] By choosing the positive constants as , and substituting these constants, we Therefore, by using the condition (24) along with the local stability condition (17f), the derivative becomes Clearly, is negative definite and hence is globally asymptotically stable. Obviously, is negative definite and is globally asymptotically stable.

5.
Persistence In this section the persistence of system (4) is investigated using the Gard technique [19]. First, we need to check the existence of periodic dynamics in the boundary planes. Obviously, the system (4) has three possible subsystems which are obtained in the absence of and and can be written, respectively, as Clearly, subsystems (26), (27) and (28) fall in the interior of plane, plane and plane, respectively. Now, define the function , which is a function in interior of of the plane. Then, determine the following quantity .
According to the Dulic-Bendixon criterion [20], the quantity does not have a changed sign and it is not identically zero, then there is no closed curve in the interior of of the plane. Further, since the subsystem (26) is bounded and has a unique positive equilibrium point in the interior of of the plane, that is given by ̆ ̆ , then, according to Poincare-Bendixon theorem [20], it is globally asymptotically stable, whenever it exists. and locally asymptotically stable. Similarly, by using the functions in the interior of of the and planes, which are defined as and respectively, it is easy to verify that there is no periodic dynamics in the interior of of their planes, and the positive equilibrium point of each subsystem ̅ ̅ and ̿ ̿ is globally asymptotically stable, whenever it exists, and locally asymptotically stable. Consequently, the necessary and sufficient conditions for the uniform persistence of system (4) are derived in the following theorem. Theorem 9: The system (4) is uniformly persistent if the following conditions hold ,

Bifurcation analysis
In this section, the appearance of the local bifurcation in system (4) is investigated. In a dynamical system, a bifurcation is carried out when a small smooth variation that is made to the parameter of a model gives a sudden or topical variation in its action. In general, the local stability properties of equilibrium point, periodic orbits, or other invariant sets, change at a bifurcation. Sotomayor's theorem [17] will be applied to study the occurrence of local bifurcation. We can write system (4)  . Theorem 10. Assume that condition (13c) holds with the parameter then the system (4) near FAEP undergoes a transcritical bifurcation provided that .
(32) However, it has a pitchfork bifurcation otherwise. Proof. Clearly, under the condition (13c) with , the Jacobian matrix that is given by equation (13a) has the following eigenvalues: and . So, the FAEP is a non-hyperbolic point and the Jacobian matrix can be written as Let be the eigenvectors of corresponding to . Then, we will get , where is any nonzero real number. In the same way, is the eigenvectors associated with the zero eigenvalue of . Hence, pitchfork bifurcation takes place. Theorem 11. Suppose that condition (14c) holds with the parameter ̂ then the system (4) near SAEP undergoes a transcritical bifurcation. Proof. Clearly, under the condition (14c) with , the Jacobian matrix that is given by equation (14a) has the following eigenvalues ̂ ̂ ̂ and . So, the SAEP is a non-hyperbolic point and the Jacobian matrix can be written as Similarly, as in theorem (10), direct computation gives that ( ) ; as an eigenvector corresponding to of .
as an eigenvector corresponding to of .
Hence, system (4) undergoes a transcritical bifurcation when . Theorem 12. Assume that ̆ ̆ ̆ ̆, then system (4) near PFEP undergoes a transcritical bifurcation, provided that ̆ ̆ , (33) where and are given in the proof. However, it has a pitchfork bifurcation otherwise. Proof. Clearly, when , the Jacobian matrix that is given by equation (15a) where and are given in the proof. However, it has a pitchfork bifurcation otherwise. Proof. Note that, when , the Jacobian matrix that is given by equation (16a)  where are given in equation (19a), then system (4) near PEP undergoes a saddle node bifurcation, provided that the following condition holds , (37) where all symbols are specified in the proof. Proof. Obviously, the Jacobian matrix at the PEP with can be written as Clearly, the elements of coincide with elements of that are given in equation (19a), except for the element , which is written in term of and denoted by . Straightforward computation shows that the determinant of that is denoted by in equation (19b) equals to zero ( ) when , which is positive under condition (36). Therefore, the characteristic equation that is given by equation (19b) has zero root with the other two negative real part roots, denoted by: where and are positive due conditions (18a) and (18b), and given by equation (19b). Therefore, the PEP is a non-hyperbolic equilibrium point. Now, to test the possibility of the occurrence of saddle node bifurcation, the following quantities are determined. as an eigenvector corresponding to of , where and .
Similar observation was obtained for varying . Now, increasing the parameter so that with other parameters as in set (38), leads to an extinction in the top predator too, and the solution approaches asymptotically to the TPFEP as shown typically by (Figure-3). Otherwise, it is still persistent at PEP. Further investigation of the dynamical behavior of system (4) using data set (38) is performed with varying one parameter each time to understand their effects on the solution and persistence of the system. It is observed that all the parameters have a quantitative change on the solution of system (4). This is due to the existence of more than on source of food in each level, which makes the extinction of any species difficult through using only one parameter. Therefore, in the following, we will investigate the system under the effects of varying of multi parameters simultaneously. For the data set (38) with and , it is noticed that the solution of system (4) approaches asymptotically to SAEP, as shown in (Figure-4). Note that it is simple to prove that the conditions (13c) and (13d) are held. Moreover, for the data set (38) with and , it is noticed that the solution of system (4) approaches asymptotically to SAEP, that is given by as shown in (Figure-5). Clearly, the data used in (Figure-5) satisfy the conditions (14c) and (14d). Again, for the data set (38) with and , the system approaches asymptotically to the PFEP, that is given by as shown in (Figure-6). Direct computation shows that the data used in (Figure-6) satisfy the condition (15f). Finally, the solution of system (4) approaches asymptotically to the MPFEP, that is given by using the data set (38) with and as shown in (Figure-7). Again, the data used in (Figure-7) satisfied the conditions (16f) and (16g).

Discussion
In this paper, a food web model incorporating a prey refuge that depends on both prey and predator species is proposed and studied. The food is consumed according to Lotka-Volterra functional response. Moreover, the intermediate predator grows logistically by the addition of favorite food at the lower level. The top predator behaves as a generalist predator and preys upon both species in the lower level and second level. All the properties of the solution of system (4) are investigated. It is observed that the system has seven nonnegative equilibrium points. The local stability conditions for each point are constructed. The global dynamics, whenever exists, is investigated too. The dynamics of all possible subsystems is also studied and it is observed that there is no periodic dynamics in the boundary planes. On the other hand, the persistence conditions of the food web system are established too. It is observed that the system has a wide range of persistence due to the existence of the prey refuge that is depending on both prey and predator species as well as the multisource of food for each species. The probability of occurrence of local bifurcation around the non-hyperbolic equilibrium point is also discussed. It is observed that the system has multi-types of bifurcations which may occur near the equilibrium points. On the other hand, numerical simulation is used to confirm our obtained results. It is observed that the system has only one type of attractors that is a stable point, while periodic dynamics does not exist even in the boundary planes. This indicates that the existence of the prey refuge that is depending on both prey and predator species is a stabilizing factor on the dynamics of food web model and extends the range of the parameters at which the system persists.