Abstract
The aim of this paper is to study the dynamics of predator-prey interaction in a chemostat to determine whether including a discrete delay to model the time between the capture of the prey and its conversion to viable biomass can introduce oscillatory dynamics even though there is a globally asymptotically stable equilibrium when the delay is ignored. Hence, Holling type I response functions are chosen so that no oscillatory behavior is possible when there is no delay. It is proven that unlike the analogous model for competition, as the parameter modeling the delay is increased, Hopf bifurcations can occur.
1. Introduction
The chemostat, also known as a continuous stir tank reactor (CSTR) in the engineering literature, is a basic piece of laboratory apparatus used for the continuous culture of microorganisms. It has potential applications for such processes as wastewater decomposition and water purification. Some ecologists consider it a lake in a laboratory. It can be thought of as three vessels, the feed bottle that contains fresh medium with all the necessary nutrients, the growth chamber where the microorganisms interact, and the collection vessel. The fresh medium from the feed bottle is continuously added to the growth chamber. The growth chamber is well stirred and its contents are then removed to the collection vessel at a rate that maintains constant volume. For a detailed description of the importance of the chemostat and its application in biology and ecology, one can refer to [1, 2].
The following system describes a food chain in the chemostat where a predator population feeds on a prey population of microorganisms that in turn consumes a nonreproducing nutrient that is assumed to be growth limiting at low concentrations
Here represents the concentration of the growth limiting nutrient, the density of the prey population, and the density of the predator population. Parameter denotes the concentration of the growth limiting nutrient in the feed vessel, the dilution rate, the growth yield constant, the sum of the dilution rate and the natural species specific death rate of the prey (predator) population, respectively. Here denotes the functional response of the prey population on the nutrient and denotes the functional response of the predator on the prey.
Butler et al. [3] considered the coexistence of two competing predators feeding on a single prey population growing in the chemostat. As a subsystem of their model, they studied the global stability of system (1.1) with both and taking the form of Holling type II. They proved that under certain conditions the interior equilibrium is globally asymptotically stable with respect to the interior of the positive cone. However, they also proved that for certain ranges of the parameters there is at least one nontrivial limit cycle and conjectured that the limit cycle is unique and would be a global attractor with respect to the noncritical orbits in the open positive octant. This conjecture was partially solved by Kuang [4]. He showed that there is a range of parameters for which a unique periodic orbit exists and roughly located the position of the limit cycle.
Bulter and Wolkowicz [5] studied predator mediated coexistence in the chemostat assuming . Model (1.1) was studied as a submodel. For general monotone response functions, Bulter and Wolkowicz showed that (1.1) is uniformly persistent if the sum of the break even concentrations of substrate and prey is less than the input rate of the nutrient . However they showed that it is necessary to specify the form of the response functions in order to discuss the global dynamics of the model. If is modelled by Holling type I or II and by Holling type I, Bulter and Wolkowicz proved that (1.1) could have up to three equilibrium points and that there is a transfer of global stability from one equilibrium point to another as different parameters are varied making conditions favorable enough for a new population to survive. In this case, there are no periodic solutions. However, even if is given by Holling type I, if is given by Holling type II, they showed that a Hopf bifurcation can occur in (1.1), and numerical simulations indicated that the bifurcating periodic solution was asymptotically stable.
We include a time delay in (1.1) to model the time between the capture of the prey and its conversion to viable biomass. Our aim is to show that such a delay can induce nontrivial periodic solutions in a model where there is always a globally asymptotically stable equilibrium when delay is ignored, and hence no such periodic solutions are possible otherwise. For this reason we select the response functions of the simplest form; that is, we choose the Holling type I form for both and , so that (1.2) always has a globally asymptotically stable equilibrium when the conversion process is assumed to occur instantaneously. It is interesting to note that in the analogous model of competition between two species in the chemostat, delay cannot induce oscillatory behavior for any reasonable monotone response functions (see Wolkowicz and Xia [6]).
With delay modelling the time required for the predator to process the prey after it has been captured, the model is given by
For ,
Here variables , , , and parameters , , , , , , and have the same interpretation as for model (1.1). Note therefore that , and . The additional parameter is a nonnegative constant modelling the time required for the conversion process. Hence, represents the concentration of the predator population in the growth chamber at time that were available at time to capture prey and were able to avoid death and washout during the units of time required to process the captured prey.
We analyze the stability of each equilibrium and prove that the coexistence equilibrium can undergo Hopf bifurcations. Numerical simulations appear to show that (1.2) can have a stable periodic solution bifurcating from the coexistence equilibrium as the delay parameter increases from zero. This periodic orbit can then disappear through a secondary Hopf bifurcation as the delay parameter increases further.
2. Scaling of the Model and Existence of Solutions
Suppose that functions and are of Holling type I form, that is, () and (). System (1.2) reduces to
Introducing the following change of variables gives:
With this change of variables, omitting the ’s for convenience, system (2.1) becomes
where and , with initial data given by (1.3). For biological significance, a point is assumed to be an equilibrium point of (2.3) only if all of its components are nonnegative.
Let . Model (2.3) reduces to a special case of the model considered in [7]. If , the model has only one equilibrium point and it is globally asymptotically stable. If and , the model has a second equilibrium point and it is globally asymptotically stable. When , the model has a third equilibrium point and it is the global attractor. Therefore, model (2.3) has no periodic solutions when the time delay is ignored. If is of Holling type II form, Butler and Wolkowicz [5] proved that a Hopf bifurcation is possible resulting in a periodic solution for a certain range of parameter values. We emphasize again here, that it is for this reason that in this paper we restrict our attention to the simplest case for both response functions, that is, Holling type I, in order to see whether delay can be responsible for periodic solutions in (1.2).
Theorem 2.1. Assuming , then there exists a unique solution of (2.3) passing through with , and for . The solution is bounded. In particular, given any , for all sufficiently large .
Proof. For , one has , and . System (2.3) becomes
a system of nonautonomous ordinary differential equations with initial conditions , , and . Since the right-hand side of (2.4) is differentiable in both and , by Theorems , , and Corollary in Miller and Michel [8], there exists a unique solution defined on satisfying (2.4). By using the method of steps in Bellman and Cooke [9], it can be shown that the solution through is defined for all .
Now we prove for all . From the first equation of (2.3),
Proceed using the method of contradiction. Suppose that there exists a first such that and for . Then . But from the first equation of (2.3)
a contradiction.
To prove for , assume there is a first such that , and for . Divide both sides of the second equation of (2.3) by and integrate from to , to obtain
contradicting .
To show that is positive on , suppose that there exists such that and for . Then From the third equation of (2.3), we have
a contradiction.
To prove the boundedness of solutions, define
It follows that
where the first inequality holds since , , and . It follows that
Therefore, the solution is bounded, and given any , for all sufficiently large .
3. Equilibria and Stability
Model (2.3) has three equilibrium points: , and
We call the washout equilibrium, the single species equilibrium, and the coexistence equilibrium. For the sake of biological significance, exists (distinct from ) if and only if its third coordinate , that is, , or equivalently, lies between and , where
Note that if , the equilibrium does not exist for any (0), and if , then .
The linearization of (2.3) about an equilibrium is given by
The associated characteristic equation is given by
Direct calculation of the left-hand side of (3.4) gives
For convenience, define
Theorem 3.1. Equilibrium is stable if and unstable if .
Proof. Evaluating the characteristic equation at gives The eigenvalues and are both negative. The third eigenvalue is . Therefore the equilibrium is stable if and unstable if .
Remark 3.2. If , then there is only one equilibrium, . If , equilibrium also exists.
Lemma 3.3. Assume . The characteristic equation evaluated at has two negative eigenvalues, and the remaining eigenvalues are solutions of In addition, the characteristic equation evaluated at has zero as an eigenvalue if and only if .
Proof. Assume . Equilibrium exists. Consider the characteristic equation at . Since at ,
where and . Therefore, and have negative real parts. The rest of the eigenvalues are roots of (3.8).
Assuming that is a root of (3.8), we have
Solving for gives
Theorem 3.4. Assume that , , , , and so that . Equilibrium is locally asymptotically stable if and unstable if . If , then equilibrium is globally asymptotically stable for .
Proof. Assume that . Assumptions , , and imply , or equivalently . By Lemma 3.3, to prove that equilibrium is locally asymptotically stable, one only needs to show that (3.8) admits no root with nonnegative real part.
Consider the real roots of (3.8) first. Note that . Equation (3.8) has no solution for . Otherwise the left-hand side would be less than zero, but the right-hand side would be greater than zero. Assume . The left-hand side of (3.8) is a monotone increasing function in both and , takes value at , and goes to positive infinity as or . By Lemma 3.3, when , then is a solution of (3.8). Thus for , any real root of (3.8) must satisfy .
For any , we have and . Therefore there exists at least one such that is a solution of (3.8). Equilibrium is unstable if .
In what follows, we prove that if all complex eigenvalues of (3.8) have negative real parts. Suppose that is a solution of (3.8). Using the Euler formula, we have
Equating the real parts and imaginary parts of the equation, we have
Squaring both equations, adding, and taking the square root on both sides give
The left-hand side of (3.14) is monotonically increasing in , , and provided that . Since (3.14) has solution , at , any roots of (3.14) must satisfy since . Hence . Therefore (3.8) has no complex eigenvalue with nonnegative real part and so is locally asymptotically stable for .
Assume that . Now we prove that is globally asymptotically stable when , or equivalently . In this case, choose small enough such that By Theorem 2.1, for such , there exists a so that for . Hence, for , In Example of Kuang ([10, page 32]), choose , , , and . We obtain . Therefore as . Let . Noting , from (2.3), we have . Multiply by the integrating factor , . Integrating both sides from to gives
If , then . Therefore . If , by L'Hôspital's rule,
since is bounded and . It again follows that . Hence
We show that and . First assume that the limits exist, that is, and . From (2.3), we know that and are uniformly continuous since , , and are bounded. By Theorem A.3, it follows that and . Note that . Letting in (2.3) gives
Either or . Assume that , that is, and . Note that . There exists such that . For such , there exists a sufficiently large so that and . Recalling that , by (2.3)
for all sufficiently large . Therefore it is impossible for to approach from above giving a contradiction. Therefore, we must have .
Now suppose that the limits do not exist. In particular if does not converge, then let and . By Lemma A.2 in the appendix, there exists and such that
From (2.3),
Noting that , we have . Since , . By (3.17), . Therefore . Similarly we can show that . This implies that , a contradiction.
Since converges and converges, then must also converge. Hence and . It follows that is globally asymptotically stable.
4. Hopf Bifurcations at Assuming
Now consider the stability of . The characteristic equation at is
By assumption , and so
where
The characteristic equation at has one eigenvalue equal to and the others are given by solutions of the equation
Lemma 4.1. Assuming , , and so that , then has no zero eigenvalue for .
Proof. Assume that . By the method of contradiction, suppose that there exists a zero root of (4.4). Therefore Noting that if and only if , for any , a contradiction.
Lemma 4.2. Assume , , . Equilibrium is asymptotically stable when .
Proof. For , (4.4) reduces to Both coefficients are positive, since and implies . Therefore, all the roots of the characteristic equation have negative real parts.
Lemma 4.3. As is increased from , a root of (4.4) with positive real part can only appear if a root with negative real part crosses the imaginary axis.
Proof. Taking and in Kuang [10, Theorem , page 66] gives
Therefore, no root of (4.4) with positive real part can enter from infinity as increases from . Hence roots with positive real part can only appear by crossing the imaginary axis.
For , assuming () is a root of ,
Substituting into (4.10) gives
Separating the real and imaginary parts, we obtain
Solving for and gives
Noting , squaring both sides of equations (4.13), adding, and rearranging gives
Solving for , we obtain two roots and :
Define conditions () and () as follows:
Lemma 4.4. If () holds for all in some interval , then (4.14) has two positive roots for all with when all the inequalities in () are strict. If () holds for all in some interval , then (4.14) has only one positive root, for all . If no interval exists where either () or () holds, then there are no positive real roots of (4.14).
Define the interval
When the end points of are real and , define
We prove that holds for any .
From ,
If , then . It follows that
Therefore,
Theorem 4.5. Assume and , then is not empty, and for any , but , condition holds and . If , then .
Proof. For any , we have , and therefore
Hence,
Therefore, . Since , it follows that
Hence,
From , we have . Therefore,
and so is not empty. Noting and , for any , but , we have .
In what follows, we intend to show that for any such , condition holds. From (4.3),
Since , to show that the first inequality in holds, it suffices to show that the factor on the right-hand side of the above expression is positive. Since , , and
Since for ,
For any ,
Hence,
Next consider the second inequality in . For , since , . Therefore, . For
Finally,
where
Since ,
It follows that . Again noting that ,
Hence, for any , we have . This leads to
Therefore, holds for any . By Lemma 4.4, both and .
If , we have . Noting (4.25), we obtain
By (4.15), it follows that and .
Now we define interval and prove that holds on
In the following theorem, we consider the case that parameters are chosen so that
Theorem 4.6. Assume and . Interval given by (4.41) is not empty. For any , holds and hence .
Proof. Assume . Letting
then
is an increasing function of and . implies that Therefore
This gives
By assumption , we obtain
Noting that and recalling the definition of given in (3.2),
Therefore, given by (4.41) is not empty. For any , noting and , we have .
In what follows, we intend to show for any , or equivalently , holds. For any , by (4.44), it follows that and so . Hence,
For any , we have , since and imply that
Therefore,
Condition holds. By Lemma 4.4, .
Next, to determine whether (4.2) has a pair of pure imaginary eigenvalues, we consider
We obtain
If there exists satisfying (4.52), then (4.2) has a pair of pure imaginary roots . A necessary condition for (4.52) to have solutions is . Otherwise, , and the second equation of (4.52) becomes . However, for any , we have , since
Hence, the second equation of (4.52) has no solution. Assume for and and denote the right-hand sides of (4.52) by
Define functions
Lemma 4.7. Assume and . For any given by (4.19), there exists and with (j=1,2) such that
Proof. For any , by Theorem 4.5, and . It is easy to see that and . There are two roots of , and Hence, for . For any , as shown in Theorem 4.5, . This implies . Therefore, and . The function is monotonically increasing for , since Since , Also, . Therefore, there exists a unique , such that . Solving for and noting that , it is easy to see that Then, for any . Since , . Therefore, for any . Since is a positive root of , we have , which implies that . Therefore, , and so . In fact, since Thus, is defined and . Since , Hence, satisfies (4.56). From (4.61), , and so . Since is continuous on the interval and is closed, there exists such that . Similarly we can prove the existence of .
Lemma 4.8. Assume and . For any given by (4.41), there exists and such that and satisfies (4.56) for .
Proof. For any , by Theorem 4.6, only .
As in Lemma 4.7, we have for . For any , as shown in Theorem 4.6, . Letting
we have
is an increasing function and . Since , it follows that . Therefore, for any , we obtain , or equivalently . Since
is monotonically increasing for any . For any ,
which implies that . Hence,
For , , since
As in the proof of Lemma 4.7, there exists a unique such that . Then . Therefore, for any . The rest of the proof is similar to the proof of Lemma 4.7. Furthermore, , since for any .
Theorem 4.9. Consider system (2.3) with . (1)Suppose , , and given by (4.22). For and , is nonnegative and there exists and such that and satisfies (4.56). If there exists such that intersects at some , then (4.4) has a pair of pure imaginary eigenvalues . System (2.3) undergoes a Hopf bifurcation at provided .(2)Suppose , , and given by (4.41). For , only is positive and there exists and such that and satisfies (4.56) for . If there exists such that intersects at some , then (4.4) has a pair of pure imaginary eigenvalues . System (2.3) undergoes a Hopf bifurcation at provided .
Proof. Assume in system (2.3).Case 1. Suppose . By Theorem 4.5, for . By Lemma 4.7, there exists and such that and satisfies (4.56). Assume that there exists a positive integer such that for some integer . Then system (4.52) has one solution . Equation (4.4) has a pair of pure imaginary eigenvalues .
In what follows, we show that the conditions required for a Hopf Bifurcation (see Theorem A.1 in the appendix) are satisfied by the linearization (3.3) of (2.3) at . In (A.1), choosing as the bifurcation parameter and letting
the linearization (3.3) of (2.3) at is of the form (A.1). Taking to be any positive real number and , hypothesis in the Hopf Bifurcation Theorem holds, since
for all and .
The characteristic equation (4.4) of (3.3) at has a pair of pure imaginary eigenvalues and no other root of (4.4) is an integral multiple of . Hence the hypothesis in the Hopf Bifurcation Theorem holds. Therefore, (2.3) undergoes a Hopf bifurcation at when provided .Case 2. Suppose . By Theorem 4.5, only . By Lemma 4.8, there exists and such that and satisfies (4.56). Assume there exists such that for some integer . Then system (4.52) has one solution . Equation (4.4) has a pair of pure imaginary eigenvalues . The rest of the proof is similar to that of Case 1 when .
Corollary 4.10. Consider system (2.3) with . (1)Suppose , , and given by (4.22). For , is nonnegative and there exists and such that and satisfies (4.56) for . If there exists a positive integer such that and , then intersects at least once at some . System (2.3) undergoes a Hopf bifurcation at provided .(2)Suppose , , and defined in (4.41). For , only is positive. There exists and such that and satisfies (4.56) for . If there exists a positive integer such that , then for any , intersects at least once at some . System (2.3) undergoes a Hopf bifurcation at provided .
Proof. Assume in system (2.3).Case 1. Suppose . By Theorem 4.5, for . By Lemma 4.7, there exists and such that and satisfies (4.56). Assume that there exists a positive integer such that and . For such , By the Mean Value Theorem, there exists such that . By Theorem 4.9. Case 1, the conclusion follows.Case 2. Suppose . By Theorem 4.5, only . By Lemma 4.8, there exists and such that and satisfies (4.56). Assume that there exists a positive integer such that . By (4.41), . Therefore . For , By the Mean Value Theorem, there exists such that . By Theorem 4.9. Case 2, the conclusion follows.
Corollary 4.11. Consider system (2.3) with . Assume and . If , then , where was defined in (4.22). For any , is nonnegative and there exists and such that and satisfies (4.56) for . If there exists a positive integer () such that , then for any , intersects at least once at some . System (2.3) undergoes a Hopf bifurcation at provided .
Proof. Assume . Then By (4.22), For any , by Theorem 4.5, for . By Lemma 4.7, there exists and such that and satisfies (4.56). Noting , . Assume there exists a positive integer () such that . For any , and . By Corollary 4.10, the conclusion follows.
5. Numerical Results
This section includes bifurcation diagrams involving the interior equilibrium and numerical simulations of periodic solutions of the predator-prey model in the chemostat.
5.1. Variation of Eigenvalues
To study the stability switches of , DDEBIFTOOL (see [11, 12]) was chosen to illustrate how the real part of the eigenvalues of (4.2) changes as parameters and vary.
First fix parameters , , and . Taking as the bifurcation parameter and varying it from to , the real part of the eigenvalues with largest real part of (4.2) was plotted in Figure 1. At and , there is either a zero eigenvalue or a pair of pure imaginary roots. For , all eigenvalues have negative real parts. For example, taking , Figure 2(a) shows that the eigenvalues of (4.2) with largest real parts (the ones in the circle) have negative real parts. Note that due to the scaling, the eigenvalues in the circle seem to be indistinguishable from zero. But in fact, they are a pair of complex eigenvalues with real parts slightly less than zero. DDEBIFTOOL can keep track of the occurrence of a pair of pure imaginary eigenvalues as varies in the neighborhood of . Figure 2(b) clearly shows that there is a pair of pure imaginary eigenvalues. Hence, Hopf bifurcation is possible. Note that by continuation, the pair of eigenvalues with largest real parts in Figure 2(a) for becomes the pair of pure imaginary eigenvalues in Figure 2(b) for .
Finally fix all parameters as before and vary both and . In Figure 3, we plot the Hopf bifurcation diagram in and parameter space. The curve at the left upper corner is . For any pair below that curve, a coexistence equilibrium exists (i.e., all components are positive). For any pair on the closed curve, there is a Hopf bifurcation. Inside the closed curve, there is a periodic solution surrounding . For any outside the closed curve and below , the coexistence equilibrium is stable.
5.2. Simulations Demonstrating Hopf Bifurcations
In this section, we illustrate Theorem 4.9 for system (2.3). Take and let vary. We choose parameters and for Case 1 (see Figures 4–12), and and for Case 2 (see Figures 13–18).
(a)
(b)
(a)
(b)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Case 1. Note that is given by (4.22). Since ,
Also, . Therefore . By Theorem 4.9, is positive and satisfies (4.56) for any and .
Figure 4 shows that intersects at some with and . We see that has no intersection with for and . By Theorem 4.9, (4.4) has two distinct pairs of pure imaginary eigenvalues . Next we need to check if .
As in Beretta and Kuang [13], we can define
Any zero of is an intersection of and and vice versa. By (4.10) in Beretta and Kuang [13] and noting that , we have the relation
where we take for and for .
From Figure 5(a), it is observed that has only one zero at with Hence, By Theorem 4.9, system (2.3) undergoes a Hopf bifurcation at . Similarly from Figure 5(b), has only one zero for and (2.3) undergoes a Hopf bifurcation at .
Next we used MATLAB to simulate solutions of model (2.3) for several values of . For each fixed delay , we chose initial data , , and for . From Figure 6, we can see that the equilibrium is stable if . As delay increases past , where a Hopf bifurcation occurs, a pair of complex eigenvalues of (4.4) enters the right-half plane. The equilibrium loses its stability and a periodic solution bifurcates from (see Figures 7 and 8). As we increase the delay further to , the periodic solution still exists and remains stable (see Figures 9 and 10). However, as the delay increases further, past , the stable periodic solution disappears in a second Hopf bifurcation, and regains stability (see Figure 11). We provide a bifurcation diagram illustrating the change in dynamics as varies (see Figure 12). For any , there is an orbitally asymptotically stable periodic solution.
Case 2. Take and . For such parameters, . By Theorem 4.9, is positive and satisfies (4.56) for and . Figure 13 shows that intersects twice. To distinguish these intersections, denote them as and . On the other hand, has no intersection with .
From Figure 14,
By (5.3),
By Theorem 4.5, system (2.3) undergoes a Hopf bifurcation at and at . For less than , the equilibrium is asymptotically stable (see Figure 15). For greater than , but less than , there is an orbitally asymptotically stable periodic solution surrounding the equilibrium (see Figures 16 and 17). At , there is a second Hopf bifurcation, where the periodic solution coalesces with . For , the periodic orbit no longer exists and regains stability (see Figure 18) until it disappears when .
Appendix
Preliminary Results
To establish the existence of periodic solutions in autonomous delay differential equations, one of the simplest ways is through Hopf Bifurcation. Below is a general Hopf Bifurcation theorem for delay differential equations due to De Oliveira [14]. Before stating the theorem we require some notation.
Consider a one parameter family of neutral delay differential equations:
where and are continuously differentiable in and ( is a constant), , , and are linear in , and
for . Assume , where , , and , , , and satisfy
It is easy to see that the characteristic matrix
is continuously differentiable in and is an entire function of . Making the following assumptions on (A.1).
() There exist constants , such that, for all complex values such that and all , the following inequalities hold: () The characteristic equation has, for , a simple purely imaginary root , , and no root of , other than , is an integral multiple of .().Now we are ready to state the Hopf bifurcation theorem for (A.1).
Theorem A.1 (Hopf bifurcation theorem, see Kuang [10, page 60]). In (A.1), assume that hold. Then there is an such that, for , , there are functions , , , , such that (A.1) has an -periodic solution , that is continuously differentiable in , and with . Furthermore, for , , every -periodic solution of (A.1) with must be of this type, except for a translation in phase; that is, there exists and such that for all .
The following lemma is usually called the Fluctuation Lemma. For a proof, see Hirsh et al. [15].
Lemma A.2. Let be a differentiable function. If , then there are sequences and such that for all
The proof of the following useful lemma can be found in [16].
Theorem A.3. Let and be a differentiable function. If exists (finite) and the derivative function is uniformly continuous on , then .
Acknowledgment
This research is partially supported by NSERC.