Bifurcation and Stability Analysis of a Food Web in a Chemostat

E-mail: srana.mthdu@gmail.com Abstract: In this study, a classical model describing a food web in a chemostat involving three species competing for non-reproducing, growth rate-limiting nutrient in which one of the competitors predates on one of the other competitors is considered. Quantitative analyses of non-negativity and boundedness of solution trajectories, dissipativity, behavior around equilibria, global stability and persistence of the model equations are analyzed. We present the global stability of equilibria by constructing a Lyapunov function. Hopf bifurcation theory is applied.


Introduction
In microbiology and population biology, the laboratory device chemostat extensively uses as research technique to culture microorganisms continuously under nutrient limitation in a controlled environment in order to study the general properties of its growth and interaction. Since one can measure the control parameter easily, the device has various applications in ecology and population biology. It can be viewed as a simple lake system in ecology while it serves as a laboratory bioreactor in chemical engineering used for investigations in genetically altered cell. As for example, the prey (bacteria) consumes nutrient (waste) while the predator (ciliates) feeds on the prey in waste water treatment process. It is of mathematical interest to construct models with chemostat. The dynamics of chemostat model with nutrient uptake is of Monod kinetics play an important role in population ecology. After the first introduction of chemostat the researchers have paid their attention to develop mathematical theories of models in it. Qualitative analyses of predator-prey models in chemostat both from the experimental and the modeling aspect describe by set of differential equations were studied by many authors (Aris and Humphrey, 1977;Butler and Wolkowicz, 1985;Drake and Tsuchiya, 1976;Huang and Zhu, 2005;Hsu et al., 1977;Jost et al., 1973a;Li, 1999;Tsuchiya et al., 1972;Wolkowicz and Lu, 1992;Waltman et al., 1980). For detailed descriptions and the complete mathematical theories for understanding models in a chemostat, we refer Smith and Waltman (1995). Moreover, using chemostat the model describing predation-prey interaction like food chain and food web and their analyses are more interesting problems in the field of mathematical biology (Alhumazi and Ajbar, 2005;Butler and Wolkowicz, 1986;El-Owaidy and Moniem, 2003;Li and Kuang, 2000;Pavlou, 1985;Zhu et al., 2002).
The interactions between species or organisms by using chemostat studied in ecology are predation and competition. In predation, one of the individuals serves as a prey for other individuals but in competition, the individuals compete for the same resources (food supply, water, etc). Competition can occur between individuals of different species or among individuals of the same species. The system comprise of the two interactions (predation and competition) is a food chain or food web. A chemostat type food chain model had been studied in (Tsuchiya et al., 1972;Jost, et al., 1973a;Drake and Tsuchiya, 1976;Butler et al., 1983) and related experiments are described in (Jost et al., 1973b;Drake and Tsuchiya, 1976). They used Michaelis-Menten type functional responses and discharge rates for the both population are equal to the washout rate of the chemostat. In the study (Butler et al., 1983), the authors obtained a global result of their model that show all competing predators coexist which circumvented competitive exclusion principle. But the authors (Li and Kuang, 2000) relaxed the assumptions that the functional responses are general monotone functions and discharge rates for the both population are distinct. They studied a predator-prey food chain model in a chemostat, where predator feeds on the prey and the prey consumes the nutrient. The authors (El-Owaidy and Moniem, 2003) also studied the same model with the extension that the predator feeds on the prey and the nutrient simultaneously. The paper (Alhumazi and Ajbar, 2005) examined a model describing the dynamic behavior of predator-prey interaction with general growth rates having no restriction on its forms. In general, it is difficult to analyze and study a complete dynamics for tri (or higher)-tropic level food chain models concern with its solutions. Consequently, the persistence or uniformly persistence analysis has been applied to describe the situation when all interacting individuals coexist in a given system. A tri-tropic level food-chain model is considered in (Chiu and Hsu, 1998). They restricted their attentions to study when the prey grows logistically in the absence of predator and the functional responses of predators are of Holling's type II and established the criterion where the top-predator goes extinct globally. The dynamics of a three-level food chain models that incorporate either Michaelis-Menten or general monotone response functions for all trophic levels had been examined by many researchers (Boer et al., 1998;El-Sheikh and Mahrouf, 2005;Freedman and Ruan, 1992;Freedman and Waltman, 1977;Hastings and Powell, 1991;Klebanoff and Hastings, 1994;Ruan and Freedman, 1991;Al-Sheikh, 2008). Bifurcation techniques also used in several studies to explore whether it is possible the coexistence of all trophic levels (Boer et al., 1998;El-Sheikh and Mahrouf, 2005;Freedman and Ruan, 1992;Al-Sheikh, 2008). Persistence analyses of models were carried out in Freedman and Waltman, 1984;Ruan and Freedman, 1991). A food web model describing predator -prey interaction consist of one predator and two competitors with a growth rate-limiting nutrient in which the predator feeds on one of the preys only had been theoretically investigated by (Butler and Wolkowicz, 1986;Wolkowicz, 2006). In (Butler and Wolkowicz, 1986), they considered the general monotone response functions, but in (Wolkowicz, 2006), the response functions that satisfy the law of mass action were assumed. On the other hand, a food web in a chemostat that incorporates one predator and two preys in which the predator predates on both of the preys and the preys consume the nutrient had been studied theoretically by (Nasrin and Rana, 2011), experimentally in lab by (Jost et al., 1973b) and numerically by (Vayenas and Pavlou, 1999). In all of their models, they employed either general monotone functions or Michaelis-Menten kinetics to describe uptake function of nutrient and growth of competitor (prey). In this study, we consider a model in a chemostat comprises of a predator individual and two competitors one of which serves as prey. We restrict our attention to examine the dynamics of the model (so called food web) where the competitors consume the nutrient and the predator feeds on one of the competitors and nutrient in the chemostat that incorporates general response functions. Our system can be viewed as resulting from the food-chain system (El-Owaidy and Moniem, 2003) by adding another competitor or from the competition system (Li, 1999) by adding a competitor preying upon one of the competitors. In this study we confine our interest to explore the criterion for global stability and persistence of the model for which all the individuals coexists in a steady state for all future time. This model can be viewed as a system experimentally in lab by considering two bacteria as prey, a protozoan as predator and glucose as the limiting nutrient in the chemostat. This paper is organized as follows. The food web model with general monotonic functional responses is described in section 2. In section 3, the basic properties (boundedness and dissipativity of the system and the existence and the local stability of equilibrium points) are given. Section 4 deals with the global stability of the system around the equilibria. In section 5, Hopf bifurcation theory applied to the system is presented. Finally in section 6, we carried out a short discussion.

Model Equations
Suppose that the concentration of nutrient and populations are all continuous functions of time. Let the concentration of nutrient, the competing populations and the predator population be denoted by s(t), x 1 (t), x 2 (t)and y(t) respectively at time t. The model equations take the following form (a system of differential equations): In the system (1), s 0 denotes the input nutrient concentration, p i (s) represent the ith competitor's specific per capita growth rate, η i are the yield constants and D is the washout rate of the chemostat. In the model equations (1), the population x 1 (t) can be treated as both a prey and a competitor. Here, q(x 1 ) is the predator's growth rate on prey population and γ is the yield factor of predator feeding on prey.
The basic assumptions on the response functions p i and q of the model (1) are: , : where, p i , q are continuously differentiable; '( ) 0 i p s ≥ for all s + ∈ R ; q′(x 1 )≥0 for all 1 x + ∈ R ; p i (0) = 0 and q(0) = 0.
The model (1) describes four-level predator-prey food web in a chemostat with non-reproducing, growth rate-limiting nutrient, constant dilution rate and input nutrient concentration, negligible death rates and instantaneous growth rates adjustment due to changes in the nutrient concentration. Furthermore, we assume that the uptake rate of nutrient is proportional to conversion rate of prey and similarly uptake rate of prey is proportional to conversion rate of predator. In model (1), the two competitor populations compete exclusively for a single, growth rate-limiting nutrient and the predator population feeds exclusively on one of the competitor population as prey that would be the sole survivor in the absence of the predator population. When all population is present, the full system will be a classical food web. This model is similar to the model (Wolkowicz, 2006).
First we simplify the system (1) by choosing the dimensionless variables. Let: The aforesaid assumptions hold good for system (2) also. Therefore, without loss of generality, one can consider the system (2) instead of the system (1) and can always reinterpret their findings in terms of original variables.

Boundedness and Dissipativity
In this section, we shall prove that the solution trajectories of model (1) initiated in + 4 R are bounded, positively invariant and dissipative for all future time.

Theorem 1
Let the region Γ be defined by: For: where: ( ) For the y equation: Finally, we define the function Ψ(t) = s(t) + x 1 (t) + x 2 (t) + y(t) and taking time derivative along the trajectories of system (2) yields Therefore, the trajectories of the system (2) are bounded uniformly and hence the system (2) is dissipative (El-Sheikh and Mahrouf, 2005). This completes the proof.

Equilibria: Existence and their Local Stability
The equilibrium points of the system (2) are the solutions of the following system of equations: with the assumed hypotheses of functions p i , q, i = 1,2,3. The eight possible equilibria of the form E(s, x 1 , x 2 , y) are given as follows: • Extinction of all populations: E 1 (1,0,0,0) • Survival of population x 1 only: E 2 (λ 1 , 1-λ 1 , 0,0) x ɶ are defined as the unique solutions of 3 1 ( ) ( ) 1 p s q x + = and they satisfy the equation: x is defined as the unique solution of The values λ i represent the break-even concentration of nutrient. Now we say that an equilibrium point will exist if its components are non-negative. The washout equilibrium point E 1 (1,0,0,0) always exists. Since p i is increasing with p i (0) = 0, λ i exist and satisfying: In this case, there are equilibrium points E 2 , E 3 and E 5 respectively, otherwise, no such equilibrium points exist. Now E 4 is a equilibrium point provided 1 . Similar conditions hold for E 7 . In case of the equilibrium point E 6 , since p 3 ,q are increasing with x ɶ exist and satisfy: For the existence of E 6 , we have two cases as in (El-Owaidy and Moniem, 2003): , then no condition is necessary Case 2: Therefore, E 6 exists if and only if: Next, the local stability analysis of the eight equilibrium points will be investigated by finding the eigenvalues of the associated Jacobian matrices. The Jacobian matrix due to the linearization of (2) about an arbitrary equilibrium The dynamic behavior of plane equilibria E 1 , E 2 , E 3 , E 4 , E 5 and E 7 have been stated in the form of Theorems [2-7] respectively. We left the proofs as these directly follow from the Routh-Hurwitz criteria. Theorem 3 , then E 2 is hyperbolic saddle and is repelling locally in both directions of x 2 and y. In particular, the stable manifold W + and the unstable manifold Wboth are of dimensions two and given by  and is repelling locally in direction of x 1 . Next, we perform the linearized stability analyses for the rest point E 6 . The Jacobian matrix due to linearization of system (2) around E 6 is given by: where: The eigenvalues of the matrix J(E 6 ) are given by: We summarize above discussion as follows.

Theorem 8
Suppose (8) where: a TraceJ E n n a n n n n n n n n n a n n n n n n n n n n n n n n a J E n n n

Global Stability Analysis
In the last section, we showed the existence of all equilibria and their local stability. Here, we shall present the global stability analysis of system (2) around the equilibria.

Proof
The proof is very straightforward and similar to (Li and Kuang, 2000;Nasrin and Rana, 2011). Now we show that the system (2) It is clear that C is a real symmetric matrix and ( ) where t stands for transpose.
Define two sets as: Then it is obvious that the sets are identical, i.e., if the matrix C is a negative definite.

Theorem 11
The equilibrium point 2 1

Proof
The proof is very straightforward and similar to (Nani and Freedman, 2000). Now, we shall deduce the global asymptotic stability criteria of along the trajectories of the system (2) and set: ( ) Then: 2 2 2 11 1 22 2 33 3 12 1 2 13 1 3 It is clear that B is a real symmetric matrix such that Thus B is negative definite if: Follows from the Routh-Hurwitz criteria and Lemma 6.1 of (Nani and Freedman, 2000).

Theorem 12
The equilibrium point We left the proof as it is very straightforward and similar to (Nani and Freedman, 2000).

Hopf-Andronov-Poincare Bifurcation
In this section, we shall show that the system (2) undergoes a Hopf-Andronov-Poincare bifurcation by using µ as a bifurcation real parameter. Without loss of generality, suppose that p 1 is a function of s and µ. Then system (2) becomes: be the set of equilibria of system (29). Suppose that G is a sufficiently large open set containing E i ∈P µ such that F(P µ ) = 0 for µ ∈ R . We are now interested to vary µ in order to observe the qualitative changes of solution trajectories near P µ .

Proof
If Theorems [3-7] hold respectively, then E 2 , E 3 , E 4 , E 5 and E 7 are repeller and their stable manifolds lie along an axis. for any values of µ is given by: where, J µ (E 6 ) is defined as in (11) and its nonzero elements are defined as in (12), with 1 ( ) p s ɶ and 1 '( )  (13). By the Routh-Hurwitz criteria, all the roots of (31) will have negative real parts if: It is required for a Hopf bifurcation that either R 1 or R 1 must be violated. Suppose that b 1 >0, b 3 >0, then clearly (31) has two purely imaginary roots if and only if µ µ = , we have (Freedman and Ruan, 1992, p.80): The roots of (33)  For the application of Hopf's bifurcation theory to the system (29) (Marsden and Mckracken, 1976), it is required to satisfy the transversality condition: Substituting r 2 (µ) = α(µ)+ iβ(µ) and r 3 (µ) = α(µ)+ iβ(µ) into (31) and calculating the derivatives with respect to µ, we obtain: where: ( ) Hence there is a Hopf bifurcation at * 1 µ µ = . We have the following result:
where, J µ (E 8 ) and its nonzero elements n ij are defined as in (14) and (15)   and R 4 is violated. Then system (29) exhibits a Hopf bifurcation in the first orthant, leading to a family of periodic solutions trajectories that bifurcate from E 8 when µ varies in a small neighborhood of * 2 µ µ = .

Discussion and Conclusion
We showed a detailed mathematical analysis of a classical food web model describing three species interaction competing for a single nutrient in which one of the competitors predates on one of the other competitors. We assumed general monotone response functions for the nutrient uptake of the competitors and growth rate for predator. The considered model is more general and its analysis extends the works of the system (El-Owaidy and Moniem, 2003;Li and Kuang, 2000;Butler and Wolkowicz, 1986;Wolkowicz, 2006). We proved that the trajectories of system (2) are bounded for all future time. We showed that the equilibria for system (2) are LAS if exist by using Routh-Hurwitz criterion. We found that all the species cannot persist if break-even concentrations of the preys are close to one. This happens when E 1 is a global attractor of system (2). We proved that if exists the equilibria E 2 and E 6 are globally asymptotically stable by using Liapunov functions as in (Nani and Freedman, 2000) showing that the predator will be die out in the chemostat whatever be the initial concentrations of prey and predator. Next, when E 8 exists then all the species (predator-preys) coexist, that is, the system (2) is persistent uniformly and circumvented the conservation principle. We applied the similar technique as in (Freedman and Waltman, 1984) and presented that the system (2) passes though a Hopf bifurcation at E 6 leading to periodic solutions. Finally, we established theoretical criteria for which the system (2) undergoes a Hopf Andronov-Poincare bifurcation at the equilibrium point E 8 and possesses a family of periodic solutions that bifurcates from E 8 . Therefore, we can conclude that the results obtained in this study partially generalize and improve to the works done by (Butler and Wolkowicz, 1986;El-Owaidy and Moniem, 2003;Li and Kuang, 2000;Nasrin and Rana, 2011;Wolkowicz, 2006).