- Email: [email protected]
Simple Food Chain in a Chemostat with Distinct Removal Rates
Journal of Mathematical Analysis and Applications 242, 75᎐92 Ž2000. doi:10.1006rjmaa.1999.6655, available online at http:rrwww.idealibrary.com on
Simple Food Chain in a Chemostat with Distinct Removal Rates Bingtuan Li Department of Mathematics, Uni¨ ersity of Utah, Salt Lake City, Utah 84112
and Yang Kuang 1 Department of Mathematics, Arizona State Uni¨ ersity, Tempe, Arizona 85287-1804 Submitted by V. Lakshmikantham Received March 23, 1999
In this paper, we consider a model describing predator᎐prey interactions in a chemostat that incorporates both general response functions and distinct removal rates. In this case, the conservation law fails. To overcome this difficulty, we make use of a novel way of constructing a Lyapunov function in the study of the global stability of a predator-free steady state. Local and global stability of other steady states, persistence analysis, as well as numerical simulations are also presented. Our findings are largely in line with those of an identical removal rate case. 䊚 2000 Academic Press
Key Words: chemostat; predator; prey; food chain; persistence.
1. INTRODUCTION The chemostat is a laboratory apparatus used for the continuous culture of microorganisms. It can be used to study competition between different populations of microorganisms, and has the advantage that the parameters are readily measurable. See the monograph of Smith and Waltman w10x for a detailed description of a chemostat and for various mathematical methods for analyzing chemostat models. 1
Research partially supported by NSF Grant DMS-9306239. 75 0022-247Xr00 $35.00 Copyright 䊚 2000 by Academic Press All rights of reproduction in any form reserved.
76
LI AND KUANG
Consider a food chain in a chemostat with one predator and one prey. We assume that the predator feeds exclusively on the prey, and the prey consumes the nutrient in the chemostat. This is an interesting practical problem both mathematically and biologically. In a waste treatment process, the bacteria live on the waste Žor nutrient . while other organisms such as ciliates feed on the bacteria. The equations of interest are
S⬘ Ž t . s Ž S 0 y S Ž t . . D y
1
␥1
F1 Ž S Ž t . . x Ž t . ,
x⬘ Ž t . s x Ž t . Ž F1 Ž S Ž t . . y D 1 . y
1
␥2
F2 Ž x Ž t . . y Ž t . ,
Ž 1.1.
y⬘ Ž t . s y Ž t . Ž F2 Ž x Ž t . . y D 2 . , S Ž 0 . ) 0,
x Ž 0 . ) 0,
y Ž 0 . ) 0,
where SŽ t . denotes the concentration of nutrient at time t; x Ž t . denotes the concentration of prey at time t; y Ž t . denotes the concentration of predator at time t; S 0 denotes the input concentration of the nutrient; ␥ 1 and ␥ 2 denote the yield constants; F1 and F2 denote the specific per-capita growth rates of prey and predator organisms, respectively; D is the washout rate of the chemostat; each Di s D q i , i s 1, 2, where 1 and 2 denote the specific death rates of organisms x and y, respectively. The values D 1 s D and D 2 s D result from assuming that the death rates of x and y are negligible so that the only loss of organisms is due to ‘‘washout’’ at the same rate that the nutrient is lost. If an organism’s death rate is significant, the removal rate of this organism should be the sum of D and the death rate. We make the following assumptions on the response functions f i Fi : Rqª Rq ,
Ž 1.2.
Fi is continuously differentiable ,
Ž 1.3.
Fi Ž 0 . s 0,
Ž 1.4.
F1X Ž S . ) 0 for all S G 0
and
F2X Ž x . ) 0 for x G 0.
Ž 1.5.
By measuring concentrations of nutrient in units of S 0 , time in units of 1rD, x is units of ␥ 1 S 0 , and y in units of ␥ 1 ␥ 2 S 0 , one reduces the number
77
DISTINCT REMOVAL RATES
of parameters and obtains the following differential equations SX Ž t . s 1 y S Ž t . y f 1 Ž S Ž t . . x Ž t . , x⬘ Ž t . s x Ž t . Ž f 1 Ž S Ž t . . y D 1 . y f 2 Ž x Ž t . . y Ž t . , y⬘ Ž t . s y Ž t . Ž f 2 Ž x Ž t . . y D 2 . , S Ž 0 . ) 0,
x Ž 0 . ) 0,
Ž 1.6.
y Ž 0 . ) 0,
where D i s D irD, i s 1, 2; f 1Ž S . s D y 1 F1Ž S 0 S .; and f 2 Ž x . s Dy1 F2 Ž␥ 1 ␥ 2 S 0 x .. Clearly, f i satisfy Ž1.2. ᎐ Ž1.5.. If the functional response functions are of the Michaelis᎐Menten form and D 1 s D 2 s 1, then system Ž1.6. becomes S⬘ s 1 y S y x⬘ s x
ž
y⬘ s y
ž
S Ž 0 . ) 0,
m1 S a1 q S m2 x a2 q x
m1 S a1 q S
x,
y1 y
/
m2 x a2 q x
y,
Ž 1.7.
y1 ,
/
x Ž 0 . ) 0,
y Ž 0 . ) 0.
System Ž1.7. has been studied in w1, 4, 9, and 12x and related experiments are described in w2, 4x. In w1x, Butler et al. took advantage of the conservation principle and gave a global analysis of this model. The authors showed that the interior critical point is globally asymptotically stable if it is locally asymptotically stable, and if it is unstable then there exists a periodic orbit for Ž1.7.. The analysis in w1x critically depends on the assumption that removal rates for the prey and predator organisms are both equal to the washout rate of the chemostat. Consider the quantity T s S q x q y, where S, x, and y are solutions of Ž1.7.. Then T satisfies the differential equality T⬘ s 1 y T. In this special case, the differential equality implies that all solutions approach the plane S q x q y s 1 at an exponential rate. This in turn implies that system Ž1.7. can be reduced to a planar system by dropping the S equation and making the substitution S s 1 y x y y in the remain-
78
LI AND KUANG
ing two equations. This is often referred to as the conservation principle. However, a slight change in the removal rate of x or y destroys the form of the conservation principle and the reduction to a planar system is no longer possible. This paper is organized as follows: In the next section, we present results on the positivity and boundedness of solutions. Section 3 deals with the existence and local stability of steady states. In Section 4, we shall provide global analysis, including global stability of the boundary steady states and persistence analysis. This paper ends with a discussion section which consists of comments and numerical simulations.
2. PRELIMINARIES In this section, we shall present some preliminary results, including the positivity and boundedness of solutions. We consider first the positivity. LEMMA 2.1. The solutions SŽ t ., x Ž t ., y Ž t . of Ž1.6. are positi¨ e, and for large t, SŽ t . - 1. Proof. It is easy to see that SŽ t . stays positive. Assume the lemma is false. Let t 1 s min t : t ) 0, x Ž t . y Ž t . s 04 . Assume first that x Ž t 1 . s 0. Then y Ž t . G 0 for t g w0, t 1 x. Let A s min 0 F t F t 1 f 1Ž SŽ t .. y D 1 y Ž f 2 Ž x Ž t ..rx Ž t .. y Ž t .4 . Then, for t g w0, t 1 x, x⬘Ž t . G Ax Ž t ., which implies that x Ž t 1 . G x Ž0. e A t 1 ) 0, a contradiction. A similar argument shows that y Ž t 1 . s 0 is absurd. Finally, S⬘ - 0 for all S G 1. This proves the lemma. Define Dmax s max 1, D 1 , D 2 4
and
Dmin s min 1, D 1 , D 2 4 .
Adding the three equations in Ž1.6. yields
Ž S q x q y . ⬘ s 1 y Ž S q D1 x q D 2 y . . This leads to 1 y Dmax Ž S q x q y . F Ž S q x q y . ⬘ F 1 y Dmin Ž S q x q y . . Solving this inequality yields the following lemma. LEMMA 2.2.
For ) 0, the solutions SŽ t ., x Ž t ., y Ž t . of Ž1.6. satisfy 1 Dmax
for large t.
y F SŽ t . q xŽ t . q yŽ t . F
1 Dmin
q
Ž 2.1.
DISTINCT REMOVAL RATES
79
3. STEADY STATES AND THEIR STABILITY The washout steady state for system Ž1.6. is denoted by E1 s Ž1, 0, 0.. There is only one possible steady state involving prey organisms but not predator organisms, denoted by E2 s Ž s , Ž1 y s .rD 1 , 0. where s is defined as the unique solution of f 1Ž S . y D1 s 0 Žif it exists.. The interior steady state is denoted by Ec s Ž SU , x , x Ž f 1Ž SU . y D 1 .rD 2 . where x is defined as the unique solution of f 2 Ž x . y D2 s 0 Žif it exists., and S* is defined as the unique solution of 1 y S y f 1Ž S . x s 0
Ž 3.1.
with S* g Ž0, 1.. One can see that no steady state can exist where there is a predator species but no prey species. We say that E2 or Ec does not exist if any one of its components is negative. We now discuss the existence of steady states. The washout steady state E1 s Ž1, 0, 0. always exists. Since f 1 is increasing with f 1Ž0. s 0,
s exists, satisfying 0 - s - 1 and f 1 Ž s . s D 1 m f 1 Ž 1 . ) D 1 . Ž 3.2. In this case there is a predator-free steady state E2 s Ž s , Ž1 y s .rD 1 , 0.. Otherwise, no such steady state exists. In the case where f 1Ž S . - D 1 for all S ) 0, we regard s s q⬁. Next consider the mixed-culture Žinterior. steady state Ec . Since f 2 is increasing with f 2 Ž0. s 0,
x exists, satisfying f 2 Ž x . s D 2 m lim f 2 Ž x . ) D 2 . xª⬁
Ž 3.3.
For Ec to exist, f 1Ž S*. y D 1 must be positive or S* ) s . Note that F Ž S . s 1 y S y f 1Ž S . x is decreasing in S with F Ž0. s 1 ) 0, F Ž S*. s 0, and F Ž s . s 1 y s y D 1 x . So, S* ) s if and only if Ž1 y s .rD 1 ) x . In the case where f 2 Ž S . - D 2 for all S ) 0, we regard x s q⬁. Therefore E2 exists if and only if s - 1, and Ec exists if and only if s - 1 and Ž1 y s .rD 1 ) x . Next we investigate the local stability of these steady states by finding the eigenvalues of the associated Jacobian matrices.
80
LI AND KUANG
The Jacobian matrix of Ž1.6. takes the form y1 y xf 1X Ž S . Js
xf 1X
yf 1 Ž S . f 1Ž S . y D1 y
Ž S.
yf 2X
0
0 yf 2X
yf 2 Ž x .
Ž x.
. Ž 3.4.
f 2 Ž x . y D2
Ž x.
At E1 , y1
yf 1 Ž 1 .
0 0
f 1Ž 1. y D1 0
J Ž E1 . s
0 0 . yD 2
Ž 3.5.
The eigenvalues lie on the diagonal. They are all negative if and only if f 1Ž1. y D 1 - 0 or, equivalently, s ) 1. When E2 exists, the Jacobian matrix at E2 is y1 y J Ž E2 . s
1 y s D1
1 y s D1
f 1X Ž s .
yf 1 Ž s .
f 1X Ž s .
0 yf 2
0
0
f2
0
ž
ž
1 y s D1
1 y s D1
/
/
. Ž 3.6.
y D2
The determinant of the upper left-hand 2 = 2 matrix is positive and its trace is negative, so its eigenvalues have negative real parts. The third eigenvalue of J Ž E2 . is f 2 ŽŽ1 y s .rD 1 . y D 2 , the entry in the lower right-hand corner. Therefore E2 is asymptotically stable if and only if f 2 ŽŽ1 y s .rD 1 . y D 2 - 0 or Ž1 y s .rD 1 - x . When Ec exists, the Jacobian matrix at Ec takes the form y1 y f 1X Ž S* . x J Ž Ec . s
f 1X Ž S* . x 0
yf 1 Ž S* .
ž
1y
x D2
f 2X Ž x .
/Ž
f 1 Ž S* . y D 1 .
x Ž f 1 Ž S* . y D 1 . D2
0
f 2X Ž x .
yD 2
.
0
Ž 3.7.
81
DISTINCT REMOVAL RATES
The eigenvalues of EŽ Ec . satisfy
3 q a1 2 q a2 q a3 s 0,
Ž 3.8.
where a1 s 1 q f 1X Ž S* . x q a2 s Ž 1 q f 1X Ž S* . x .
ž
ž
x D2
x D2
f 2X Ž x . y 1
f 2X Ž x . y 1
/Ž
/Ž
f 1 Ž S* . y D 1 . ,
f 1 Ž S* . y D 1 .
q x Ž f 1 Ž S* . y D 1 . f 2X Ž c . q x f 1 Ž S* . f 1X Ž S* . , and a3 s x f 2X Ž x . Ž 1 q f 1X Ž S* . x .Ž f 1 Ž S* . y D 1 . . Note that the constant term a3 is positive, so the Routh᎐Hurwitz criterion says that Ec will be asymptotically stable if and only if a1 ) 0 and a1 a2 ) a3 . We summarize the above results in the following theorem. THEOREM 3.1. If s ) 1, then only E1 exists and E1 is locally asymptotically stable. If s - 1 and Ž1 y s .rD 1 - x , then only E1 and E2 exist, E1 is unstable, and E2 is locally asymptotically stable. If s - 1 and Ž1 y s .rD 1 ) x , then E1 , E2 , and Ec exist, and E1 and E2 are unstable. Ec is locally asymptotically stable if a1 ) 0 and a1 a2 ) a3 . If D 1 s D 2 s 1, then the limit plane of Ž1.6. is ⌺ : S q x q y s 1. If we drop the S equation, then the limit system of Ž1.6. takes the form x⬘ s x Ž f 1 Ž 1 y x y y . y 1 . y f 2 Ž x . y, y⬘ s y Ž f 2 Ž x . y 1 . , x Ž 0 . ) 0,
Ž 3.9.
y Ž 0 . ) 0.
For this reduced system, if Ec exists then the Jacobian matrix at Ec is J Ž Ec . s
y x f 1X Ž S* . y Ž x f 2X Ž x . y 1 .Ž f 1 Ž S* . y 1 . y x f 1X Ž S* . y f 2 Ž x .
x Ž f 1 Ž S* . y 1 . f 2X Ž x .
0
.
Ž 3.10.
82
LI AND KUANG
The determinant of this matrix is w x f 1X Ž S*. q f 2 Ž x .xw x Ž f 1Ž S* y 1. which is positive, and the trace is y x f 1X Ž S*. y Ž x f 2X Ž x . y 1. Ž f 1Ž S*. y 1.. Therefore if y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f Ž S*. y 1. - 0 then Ec is locally asymptotically stable; if y x f 1X Ž S*. y Ž x f 2X Ž x . y 1. Ž f 1Ž S*. y 1. ) 0 then Ec is a repeller in ⌺, and in this case there is a periodic solution in ⌺ Žby an application of the Poincare᎐Bendixson ´ teorem.. We summarize these in the following theorem. f 2U Ž x .x,
THEOREM 3.2. Assume D 1 s D 2 s 1 and Ec exists Ž i.e., x q s - 1.. If y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. - 0 then Ec is locally asymptotically stable. If y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. ) 0 then Ec is unstable, and there is a periodic solution in the plane ⌺ : S q x q y s 1.
4. GLOBAL ANALYSIS In the previous section, we showed that if only E1 exists then E1 is asymptotically stable, if E1 and E2 exist then E1 is unstable and E2 is Žlocally. asymptotically stable, if E1 , E2 , and Ec exist then E1 and E2 are unstable, and in this case Ec may or may not be stable. We shall show that E1 is globally asymptotically stable if only E1 exists. The proof is very straightforward. Most importantly, we shall show that if only E1 and E2 exist, under a reasonable additional assumption E2 is globally asymptotically stable. The proof involves the construction of a Lyapunov function and the application of the Lyapunov᎐LaSalle theorem. ŽWe shall use Theorem 2.1 in Wolkowicz and Lu w13x, which is a slightly modified version of the statements given in LaSalle w6x and Hale w3x.. We shall also show that system Ž1.6. is uniformly persistent if Ec exists. The following theorem states that E1 is a global attractor if it is the only steady state Ži.e., s ) 1.. THEOREM 4.1.
If s ) 1, then all solutions of Ž1.6. satisfy lim Ž S Ž t . , x Ž t . , y Ž t . . s Ž 1, 0, 0 . .
tª⬁
Proof. Since S Ž t . - 1 for large t and f 1Ž1. y D 1 - 0 Ži.e., s ) 1., there is ␣ ) 0 such that x⬘Ž t . - y␣ x Ž t . for t sufficiently large. This shows lim t ª⬁ x Ž t . s 0. It follows from the third equation of Ž1.6. that lim t ª⬁ y Ž t . s 0. Then the first equation of Ž1.6. yields lim t ª⬁ SŽ t . s 1. The proof is complete. Note that if Ž s q x . Dmin ) 1 then s q D 1 x ) 1. That is, 1rDmin s q x implies Ž1 y s .rD 1 - x .
83
DISTINCT REMOVAL RATES
THEOREM 4.2.
If s - 1 and 1
- s q x ,
Dmin
Ž 4.1.
then all solutions of Ž1.6. satisfy lim Ž S Ž t . , x Ž t . , y Ž t . . s s ,
ž
tª⬁
1 y s D1
/
,0 .
Proof. We choose d1 ) Dmax and d 2 - Dmin such that 1 d2
- s q x
Ž 4.2.
and that, for large t, 1 d1
- SŽ t . q xŽ t . q yŽ t . -
1 d2
.
Ž 4.3.
Let
␣s1q
max
0FxF Ž1y s .rD 1
½
f 2 Ž x . Ž Ž 1 y s . rD 1 y x . x Ž D2 y f 2 Ž x . .
5
and
s1q
max
xFxF1rDmin q1
½
␣ Ž f 2 Ž x . y D2 . D2
5
.
Let C Ž u. be a continuously differentiable function and C⬘Ž u. is given by
¡0, C⬘ Ž u . s
~
if u F
 x q s y 1rd 2
¢ ,
ž
u q s y
1 d2
/
,
if
1 d2
1 d2
y s - u - x ,
u G x .
C⬘Ž u. is simply linear on w1rd 2 y s , x x. In view of Ž4.3., xqyF
1 d2
y s
if S G s .
y s ,
84
LI AND KUANG
Therefore, if S G s , then C⬘Ž x q y . s 0. The graph of C⬘Ž x q y . is shown in Fig. 1. Define the Lyapunov function V Ž S, x, y . as follows Vs
S
Ž f 1Ž . y D1 . Ž 1 y s .
s
D1Ž 1 y .
H
q x y x* y x*ln
ž
x x*
/
d
q ␣ y q CŽ x q y.
Ž 4.4.
on the set ⌿ s Ž S, x, y . : S g Ž0, 1., x, y g Ž0, q⬁., S q x q y g Ž1rd1 , 1rd 2 .4 , where x* ' Ž1 y s .rD 1. Then the time derivative of V along solutions of the differential equation is V˙ s C⬘ Ž x q y . q 1 y q
f2 Ž x . x
ž
1 y s D1
Ž 1 y s . f 1Ž S . Ž f 1Ž S . y D1 . x D1Ž 1 y S . y x q ␣ Ž f 2 Ž x . y D 2 . y D 2 C⬘ Ž x q y . y.
/
First, note that w1 y Ž1 y 2 . f 1Ž S .rD 1Ž1 y S .xŽ f 1Ž S . y D 1 . x is nonpositive for 0 - S - 1 and equals 0 for S g w0, 1. if and only if S s s or x s 0. Since C⬘Ž x q y . s 0 for S G s and C⬘Ž u. G 0 for u G 0, C⬘Ž x q y . Ž f 1Ž S . y D 1 . x is nonpositive for S g w0, 1.. Therefore the first term in V˙ is always nonpositive and equals 0 for S g w0, 1. if and only if S s s or x s 0.
FIG. 1. A graphical depiction of C⬘Ž x q y . for S g w0, 1..
85
DISTINCT REMOVAL RATES
Define h Ž S, x, y . s
f2 Ž x . x
ž
1 y s D1
y x q ␣ Ž f 2 Ž x . y D 2 . y D 2 C⬘ Ž x q y . .
/
If 0 - x F 1 y srD 1 , then f2 Ž x . x
ž
1 y s D1
yx G0
/
and
f 2 Ž x . y D2 F 0
and they are equal to 0 if and only if x s Ž1 y s .rD 1. Note that C⬘Ž x q y . is always nonnegative. By the definition of ␣ , hŽ S, x, y . - 0 for 0 F x - Ž1 y s .rD 1 and possibly hŽ S, x, y . s 0 if x s Ž1 y s .rD 1. If Ž1 y s .rD 1 - x F x , all three terms in hŽ S, x, y . are nonpositive and one can easily see hŽ S, x, y . - 0. If x ) x , then x q y ) x and therefore C⬘Ž x q y . s  . Note that, if x ) x , only the second term in hŽ S, x, y . is nonnegative. According to the definition of  , hŽ S, x, y . - 0 if x ) x . Therefore hŽ S, x, y . - 0 for x G 0 and x / Ž1 y s .rD 1 , and it is possibly 0 if x s Ž1 y s .rD 1. By Lemma 2.1 every bounded solution of Ž1.6. is contained in ⌿, and hence by Theorem 2.1 in w13x every solution of Ž1.6. approaches the set ⌳, the largest invariant subset M of ⌽ s Ž S, x, y . g ⌿ : V˙ s 04 . ⌽ is made up of points of the following forms where S g w 0, 1 x ,
Ž S, 0, 0 . ,
ž
s ,
1 y s D1
/
,y ,
Ž s , x, 0 . ,
where y g w 0, ⬁ . , where x g w 0, ⬁ . .
Since V is bounded above, any point of the form Ž S, 0, 0. cannot be in the 3 Ž -limit set ⍀ of any solution initiating in the interior of Rq . s , x, 0. g M Ž . Ž . implies that S t s , which in turn leads to 0 s S⬘ t s 1 y s y f 1Ž s . x and hence x s Ž1 y s .rD 1. Ž s , Ž1 y s .rD 1 , y . g M implies that S Ž t . s s and x Ž t . s Ž1 y s .rD 1. The second equation of Ž1.6. implies that x⬘ s 0, which yields y s 0. Therefore M s E2 4 . This completes the proof. THEOREM 4.3. If s - 1 and Ž1 y s .rD 1 ) x , then system Ž1.6. is uniformly persistent; i.e., there exists a constant ) 0, independent of initial conditions, such that lim inf S Ž t . G , tª⬁
lim inf x Ž t . G , tª⬁
lim inf y Ž t . G . tª⬁
86
LI AND KUANG
Proof. Choose
½ ½ ½
X 1 s Ž S, x, y . ; 0 F S F 1, 0 - x F Y1 s Ž S, x, 0 . ; 0 F S F 1, 0 F x F Y2 s Ž S, 0, y . ; 0 F S F
1 Dmin
1 Dmin 1 Dmin
q 1, 0 - y F
1 Dmin
5
q1 ,
5
q1 ,
q 1, 0 F y F
1 Dmin
5
q1 ,
and X 2 s Y1 j Y2 .
FIG. 2. m1 s 3.6, a1 s 0.8, D 1 s 1.4, m 2 s 3, a2 s 0.6, D 2 s 1.2. Ž S Ž0., x Ž0., y Ž0.. s Ž0.6, 0.2, 0.7.. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. In this case, s q x - 1. Clearly, the solution approaches the predator-free steady state.
DISTINCT REMOVAL RATES
87
Then X 1 and X 2 are two disjoint subsets of R 3 , X 2 is compact, X s X 1 j X 2 is also compact, and X 1 and X 2 are positively invariant for Ž1.6.. By Lemma 2.2, X 2 and X are global attractors in the union of the S y x 3 plane and S y y plane and in Rq , respectively. We prove that X 2 is a uniformly strong repeller for X 1 Žfor the definitions of a uniformly strong repeller as well as a weak repeller, see Thieme w11x.. E1 and E2 are the only steady states in X 2 . E1 is a saddle in R 3 and its stable manifold is Ž S, 0, y .; y G 04 . E2 is also a saddle in R 3 and its stable manifold is Ž S, x, 0.; x ) 04 . Therefore E1 and E2 are weak repellers for X 1. The stable manifold structures of E1 and E2 Ž E2 is a global attractor in the S y x plane. imply that they are not cyclically chained to each other on the boundary X 2 . By Proposition 1.2 of Thieme w11x, X 2 is a uniform strong repeller for X 1; that is, there are ␦ 1 ) 0 and ␦ 2 ) 0 such that lim inf t ª⬁ x Ž t . ) ␦ 1 and lim inf t ª⬁ y Ž t . ) ␦ 2 with ␦ 1 and ␦ 2 not depend-
FIG. 3. m1 s 3.6, a1 s 0.8, D 1 s 1.4, m 2 s 3, a2 s 0.6, D 2 s 1. Ž SŽ0., x Ž0., y Ž0.. s Ž0.6, 0.2, 0.7.. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. Clearly, the solution approaches a positive steady state.
88
LI AND KUANG
ing on the initial values in X 1. Applying Proposition 2.2 of Thieme w11x to the first equation of Ž1.6. yields that there is ␦ 3 ) 0 such that lim inf t ª⬁ SŽ t . ) ␦ 3 with ␦ 3 not depending on the initial values of X 1. This completes the proof.
5. DISCUSSION In this paper, we considered a food chain with one prey and one predator in the chemostat. In this model, the prey consumes the nutrient and the predator consumes the prey but the predator does not consume the nutrient. We assumed that the functional response functions are general monotone response functions and the removal rates are different. The model we considered is more general and realistic than the models in w1, 4, 9, 12x.
FIG. 4. m1 s 8.5, a1 s 0.6, D 1 s D 2 s 1, m 2 s 6, a2 s 0.6. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. The top curve depicts SŽ t ., the middle one depicts y Ž t ., and the bottom one depicts x Ž t .. The solution appears to approach a periodic solution.
DISTINCT REMOVAL RATES
89
The main difficulty we faced was the lack of a conservation principle which was lost due to the different removal rates. In the case of different removal rates, the system cannot be reduced to a two-dimensional system and we therefore must look at the full system. We found that the washout steady state E1 is the global attractor if it is the only steady state Žthis happens when s ) 1.. This confirms the intuition that both prey and predator cannot persist if the removal rate of the prey is relatively large. When E1 and the predator-free steady state E2 are the only steady states, we found that E1 is unstable and E2 is locally asymptotically stable. By constructing a Lyapunov function, we were able to show that if E1 and E2 are the only steady states, under an additional assumption Žsee Ž4.1.., E2 is a global attractor. The construction of the Lyapunov function is rather novel and nontrivial. This novel idea has been used in w7, 8x. This condition does not depend on the specific properties of the functional response functions, and it becomes necessary if Dmin is close to both D 1 and 1. The
FIG. 5. m1 s 8.5, a1 s 0.6, D 1 s 1.1, m 2 s 6, a2 s 0.6, D 2 s 1. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. The top curve depicts SŽ t ., the middle one depicts y Ž t ., and the bottom one depicts x Ž t .. The solution oscillates but eventually approaches a positive steady state.
90
LI AND KUANG
global stability of E2 implies that the predator will be washed out in the chemostat regardless of the initial density levels of prey and predator. We also showed that, when Ec exists, the prey and predator coexist in the sense that the system is uniformly persistent. In this case, a switch of the stability of the interior steady state Ec may occur. E2 is a global attractor if it is locally asymptotically stable and Ž4.1. holds. Based on our extensive simulation work where the functional response functions take the Michaelis᎐Menten form f 1Ž S . s
m1 S a1 q S
and
f2 Ž x . s
m2 x a2 q x
,
we conjecture that this steady state remains globally asymptotically stable as long as it is locally stable Žsee Fig. 2.. As E2 becomes unstable, a locally asymptotically stable interior steady state Ec bifurcates from it. Our simulation work ŽFig. 3. suggests that Ec is a global attractor if it is locally asymptotically stable. As certain parameters increase or decrease further away, Ec loses its stability and oscillatory solutions appear. These oscillatory solutions Žsee Figs. 4 and 6. appear to be the results of Hopf bifurcations. Figure 4 shows a case in which D 1 s D 2 s 1 and system Ž1.6. possesses periodic solutions. Figure 5 indicates that perturbing D 1 Žwhile keeping other parameters in Fig. 4 fixed. leads to a bifurcation. For example, changing D 1 s 1 in Fig. 4 to D 1 s 1.1 in Fig. 5 seems to destroy the periodic solutions and possibly leads to the global stability of Ec . Therefore varying the values of D 1 and D 2 may affect the dynamics of Ž1.6. in a very surprising and significant way. Next, we see how the parameters D 1 and D 2 affect the dynamics of Ž1.6. if S 0 is fixed. D 1 s 1 q 1 and D 2 s 1 q 2 , where 1 and 2 denote the scaled specific death rates of the prey and predator, respectively. ŽNote that the analysis of the model requires no assumptions on the sign of i ’s, as long as the Di ’s all remain positive. This leaves the IDi ’s open to other interpretations.. Assume that D 1 and D 2 are large enough so that s ) 1, then both prey and predator populations will be washed out Ž E1 is stable . in the chemostat. As D 1 is gradually decreased, eventually there is a bifurcation when s - 1 and Ž1 y s .rD 1 - x . In this case, E1 loses its stability and the new bifurcated steady state E2 is asymptotically stable. As D 2 is gradually decreased, the next bifurcation occurs when s - 1 and Ž1 y s .rD 1 ) x hold. In this case E2 loses its stability, and a new interior steady state Ec appears. If D 1 s D 2 s 1, the conservation principle holds; that is, the -limit sets of solutions of Ž1.6. lie in the plane ⌺ : S q x q y s 1. In this case, one can easily show that Ec Žif it exists. is locally asymptotically stable if and only if f 1X Ž S*. x q Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. ) 0. When this in-
DISTINCT REMOVAL RATES
91
FIG. 6. m1 s 4, a1 s 0.6, D 1 s 1.1, m 2 s 5, a2 s 0.5, D 2 s 1.2. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. In this case, it can be shown that Ec exists and is unstable. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. The solution oscillates and seems to approach a periodic solution.
equality is reversed, Ec will be a repeller in ⌺, and there will be at least one periodic orbit Žby an application of the Poincare᎐Bendixson theorem.. ´ Determining the number of periodic solutions is a deep mathematical problem. Kuang w5x has shown in the case of Michaelis᎐Menten-type response functions, if yf 1X Ž S*. x y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. is small and positive, then the limit cycle is unique and asymptotically stable. We guess this is true for general response functions. In w1x, it was shown that in the case of Michaelis᎐Menten-type response functions and D 1 s D 2 s 1, Ec is globally asymptotically stable if it is locally asymptotically stable. It remains open if this is true in the case of general response functions and different removal rates.
92
LI AND KUANG
REFERENCES 1. G. J. Butler, S. B. Hsu, and P. Waltman, Coexistence of competing predators in a chemostat, J. Math. Biol. 17 Ž1983., 133᎐151. 2. J. F. Drake and H. M. Tsuchiya, Predation of escherichia coli by copoda steuii, Appl. En¨ iron. Microbiol. 31 Ž1976., 870᎐874. 3. J. K. Hale, ‘‘Ordinary Differential Equations,’’ Krieger, Malabar, FL, 1980. 4. J. L. Jost, S. F. Drake, A. G. Fredrickson, and M. Tsuchiya, Interaction of tetrahymena pyriformis, escherichia coli, azotobacter vinelandii and glucose in a minimal medium, J. Bacteriol. 113 Ž1976., 834᎐840. 5. Y. Kuang, Limit cycles in a chemostat-related models, SIAM J. Appl. Math. 49 Ž1989., 1759᎐1767. 6. J. LaSalle, Some extensions of Lyapunov’s second method, IRE Trans. Circuit CT-7 Ž1960., 520᎐527. 7. B. Li, Global asymptotic behavior of the chemostat: General response functions and different removal rates, SIAM J. Appl. Math. 59 Ž1999., 411᎐422. 8. B. Li, ‘‘Analysis of Chemostat-Related Models with Distinct Removal Rates,’’ Ph.D. thesis, Arizona State University, 1998. 9. G. Sell, What is a dynamical system? in ‘‘Studies in Ordinary Differential Equations’’ ŽJ. Hale, Ed.., MAA Studies in Mathematics, Vol. 14, Math. Assoc. America, Washington, DC, 1977. 10. H. L. Smith and P. Waltman, ‘‘The Theory of the Chemostat,’’ Cambridge Univ. Press, Cambridge, UK, 1995. 11. H. R. Thieme, Persistence under relaxed point-dissipativity with an application to an epidemic model., SIAM J. Math. Anal. 24 Ž1993., 407᎐435. 12. H. M. Tsuchiya, S. F. Drake, J. L. Jost, and A. G. Fredrickson, Predator᎐prey interactions of dictyostelium discordeum and escherichia coli in continuous culture, J. Bacteriol. 110 Ž1972., 1147᎐1153. 13. G. S. K. Wolkowicz and Z. Lu, Global dynamics of a mathematical model of competition in the chemostat: General response function and differential death rates, SIAM J. Appl. Math. 52 Ž1992., 222᎐233.
Simple Food Chain in a Chemostat with Distinct Removal Rates Bingtuan Li Department of Mathematics, Uni¨ ersity of Utah, Salt Lake City, Utah 84112
and Yang Kuang 1 Department of Mathematics, Arizona State Uni¨ ersity, Tempe, Arizona 85287-1804 Submitted by V. Lakshmikantham Received March 23, 1999
In this paper, we consider a model describing predator᎐prey interactions in a chemostat that incorporates both general response functions and distinct removal rates. In this case, the conservation law fails. To overcome this difficulty, we make use of a novel way of constructing a Lyapunov function in the study of the global stability of a predator-free steady state. Local and global stability of other steady states, persistence analysis, as well as numerical simulations are also presented. Our findings are largely in line with those of an identical removal rate case. 䊚 2000 Academic Press
Key Words: chemostat; predator; prey; food chain; persistence.
1. INTRODUCTION The chemostat is a laboratory apparatus used for the continuous culture of microorganisms. It can be used to study competition between different populations of microorganisms, and has the advantage that the parameters are readily measurable. See the monograph of Smith and Waltman w10x for a detailed description of a chemostat and for various mathematical methods for analyzing chemostat models. 1
Research partially supported by NSF Grant DMS-9306239. 75 0022-247Xr00 $35.00 Copyright 䊚 2000 by Academic Press All rights of reproduction in any form reserved.
76
LI AND KUANG
Consider a food chain in a chemostat with one predator and one prey. We assume that the predator feeds exclusively on the prey, and the prey consumes the nutrient in the chemostat. This is an interesting practical problem both mathematically and biologically. In a waste treatment process, the bacteria live on the waste Žor nutrient . while other organisms such as ciliates feed on the bacteria. The equations of interest are
S⬘ Ž t . s Ž S 0 y S Ž t . . D y
1
␥1
F1 Ž S Ž t . . x Ž t . ,
x⬘ Ž t . s x Ž t . Ž F1 Ž S Ž t . . y D 1 . y
1
␥2
F2 Ž x Ž t . . y Ž t . ,
Ž 1.1.
y⬘ Ž t . s y Ž t . Ž F2 Ž x Ž t . . y D 2 . , S Ž 0 . ) 0,
x Ž 0 . ) 0,
y Ž 0 . ) 0,
where SŽ t . denotes the concentration of nutrient at time t; x Ž t . denotes the concentration of prey at time t; y Ž t . denotes the concentration of predator at time t; S 0 denotes the input concentration of the nutrient; ␥ 1 and ␥ 2 denote the yield constants; F1 and F2 denote the specific per-capita growth rates of prey and predator organisms, respectively; D is the washout rate of the chemostat; each Di s D q i , i s 1, 2, where 1 and 2 denote the specific death rates of organisms x and y, respectively. The values D 1 s D and D 2 s D result from assuming that the death rates of x and y are negligible so that the only loss of organisms is due to ‘‘washout’’ at the same rate that the nutrient is lost. If an organism’s death rate is significant, the removal rate of this organism should be the sum of D and the death rate. We make the following assumptions on the response functions f i Fi : Rqª Rq ,
Ž 1.2.
Fi is continuously differentiable ,
Ž 1.3.
Fi Ž 0 . s 0,
Ž 1.4.
F1X Ž S . ) 0 for all S G 0
and
F2X Ž x . ) 0 for x G 0.
Ž 1.5.
By measuring concentrations of nutrient in units of S 0 , time in units of 1rD, x is units of ␥ 1 S 0 , and y in units of ␥ 1 ␥ 2 S 0 , one reduces the number
77
DISTINCT REMOVAL RATES
of parameters and obtains the following differential equations SX Ž t . s 1 y S Ž t . y f 1 Ž S Ž t . . x Ž t . , x⬘ Ž t . s x Ž t . Ž f 1 Ž S Ž t . . y D 1 . y f 2 Ž x Ž t . . y Ž t . , y⬘ Ž t . s y Ž t . Ž f 2 Ž x Ž t . . y D 2 . , S Ž 0 . ) 0,
x Ž 0 . ) 0,
Ž 1.6.
y Ž 0 . ) 0,
where D i s D irD, i s 1, 2; f 1Ž S . s D y 1 F1Ž S 0 S .; and f 2 Ž x . s Dy1 F2 Ž␥ 1 ␥ 2 S 0 x .. Clearly, f i satisfy Ž1.2. ᎐ Ž1.5.. If the functional response functions are of the Michaelis᎐Menten form and D 1 s D 2 s 1, then system Ž1.6. becomes S⬘ s 1 y S y x⬘ s x
ž
y⬘ s y
ž
S Ž 0 . ) 0,
m1 S a1 q S m2 x a2 q x
m1 S a1 q S
x,
y1 y
/
m2 x a2 q x
y,
Ž 1.7.
y1 ,
/
x Ž 0 . ) 0,
y Ž 0 . ) 0.
System Ž1.7. has been studied in w1, 4, 9, and 12x and related experiments are described in w2, 4x. In w1x, Butler et al. took advantage of the conservation principle and gave a global analysis of this model. The authors showed that the interior critical point is globally asymptotically stable if it is locally asymptotically stable, and if it is unstable then there exists a periodic orbit for Ž1.7.. The analysis in w1x critically depends on the assumption that removal rates for the prey and predator organisms are both equal to the washout rate of the chemostat. Consider the quantity T s S q x q y, where S, x, and y are solutions of Ž1.7.. Then T satisfies the differential equality T⬘ s 1 y T. In this special case, the differential equality implies that all solutions approach the plane S q x q y s 1 at an exponential rate. This in turn implies that system Ž1.7. can be reduced to a planar system by dropping the S equation and making the substitution S s 1 y x y y in the remain-
78
LI AND KUANG
ing two equations. This is often referred to as the conservation principle. However, a slight change in the removal rate of x or y destroys the form of the conservation principle and the reduction to a planar system is no longer possible. This paper is organized as follows: In the next section, we present results on the positivity and boundedness of solutions. Section 3 deals with the existence and local stability of steady states. In Section 4, we shall provide global analysis, including global stability of the boundary steady states and persistence analysis. This paper ends with a discussion section which consists of comments and numerical simulations.
2. PRELIMINARIES In this section, we shall present some preliminary results, including the positivity and boundedness of solutions. We consider first the positivity. LEMMA 2.1. The solutions SŽ t ., x Ž t ., y Ž t . of Ž1.6. are positi¨ e, and for large t, SŽ t . - 1. Proof. It is easy to see that SŽ t . stays positive. Assume the lemma is false. Let t 1 s min t : t ) 0, x Ž t . y Ž t . s 04 . Assume first that x Ž t 1 . s 0. Then y Ž t . G 0 for t g w0, t 1 x. Let A s min 0 F t F t 1 f 1Ž SŽ t .. y D 1 y Ž f 2 Ž x Ž t ..rx Ž t .. y Ž t .4 . Then, for t g w0, t 1 x, x⬘Ž t . G Ax Ž t ., which implies that x Ž t 1 . G x Ž0. e A t 1 ) 0, a contradiction. A similar argument shows that y Ž t 1 . s 0 is absurd. Finally, S⬘ - 0 for all S G 1. This proves the lemma. Define Dmax s max 1, D 1 , D 2 4
and
Dmin s min 1, D 1 , D 2 4 .
Adding the three equations in Ž1.6. yields
Ž S q x q y . ⬘ s 1 y Ž S q D1 x q D 2 y . . This leads to 1 y Dmax Ž S q x q y . F Ž S q x q y . ⬘ F 1 y Dmin Ž S q x q y . . Solving this inequality yields the following lemma. LEMMA 2.2.
For ) 0, the solutions SŽ t ., x Ž t ., y Ž t . of Ž1.6. satisfy 1 Dmax
for large t.
y F SŽ t . q xŽ t . q yŽ t . F
1 Dmin
q
Ž 2.1.
DISTINCT REMOVAL RATES
79
3. STEADY STATES AND THEIR STABILITY The washout steady state for system Ž1.6. is denoted by E1 s Ž1, 0, 0.. There is only one possible steady state involving prey organisms but not predator organisms, denoted by E2 s Ž s , Ž1 y s .rD 1 , 0. where s is defined as the unique solution of f 1Ž S . y D1 s 0 Žif it exists.. The interior steady state is denoted by Ec s Ž SU , x , x Ž f 1Ž SU . y D 1 .rD 2 . where x is defined as the unique solution of f 2 Ž x . y D2 s 0 Žif it exists., and S* is defined as the unique solution of 1 y S y f 1Ž S . x s 0
Ž 3.1.
with S* g Ž0, 1.. One can see that no steady state can exist where there is a predator species but no prey species. We say that E2 or Ec does not exist if any one of its components is negative. We now discuss the existence of steady states. The washout steady state E1 s Ž1, 0, 0. always exists. Since f 1 is increasing with f 1Ž0. s 0,
s exists, satisfying 0 - s - 1 and f 1 Ž s . s D 1 m f 1 Ž 1 . ) D 1 . Ž 3.2. In this case there is a predator-free steady state E2 s Ž s , Ž1 y s .rD 1 , 0.. Otherwise, no such steady state exists. In the case where f 1Ž S . - D 1 for all S ) 0, we regard s s q⬁. Next consider the mixed-culture Žinterior. steady state Ec . Since f 2 is increasing with f 2 Ž0. s 0,
x exists, satisfying f 2 Ž x . s D 2 m lim f 2 Ž x . ) D 2 . xª⬁
Ž 3.3.
For Ec to exist, f 1Ž S*. y D 1 must be positive or S* ) s . Note that F Ž S . s 1 y S y f 1Ž S . x is decreasing in S with F Ž0. s 1 ) 0, F Ž S*. s 0, and F Ž s . s 1 y s y D 1 x . So, S* ) s if and only if Ž1 y s .rD 1 ) x . In the case where f 2 Ž S . - D 2 for all S ) 0, we regard x s q⬁. Therefore E2 exists if and only if s - 1, and Ec exists if and only if s - 1 and Ž1 y s .rD 1 ) x . Next we investigate the local stability of these steady states by finding the eigenvalues of the associated Jacobian matrices.
80
LI AND KUANG
The Jacobian matrix of Ž1.6. takes the form y1 y xf 1X Ž S . Js
xf 1X
yf 1 Ž S . f 1Ž S . y D1 y
Ž S.
yf 2X
0
0 yf 2X
yf 2 Ž x .
Ž x.
. Ž 3.4.
f 2 Ž x . y D2
Ž x.
At E1 , y1
yf 1 Ž 1 .
0 0
f 1Ž 1. y D1 0
J Ž E1 . s
0 0 . yD 2
Ž 3.5.
The eigenvalues lie on the diagonal. They are all negative if and only if f 1Ž1. y D 1 - 0 or, equivalently, s ) 1. When E2 exists, the Jacobian matrix at E2 is y1 y J Ž E2 . s
1 y s D1
1 y s D1
f 1X Ž s .
yf 1 Ž s .
f 1X Ž s .
0 yf 2
0
0
f2
0
ž
ž
1 y s D1
1 y s D1
/
/
. Ž 3.6.
y D2
The determinant of the upper left-hand 2 = 2 matrix is positive and its trace is negative, so its eigenvalues have negative real parts. The third eigenvalue of J Ž E2 . is f 2 ŽŽ1 y s .rD 1 . y D 2 , the entry in the lower right-hand corner. Therefore E2 is asymptotically stable if and only if f 2 ŽŽ1 y s .rD 1 . y D 2 - 0 or Ž1 y s .rD 1 - x . When Ec exists, the Jacobian matrix at Ec takes the form y1 y f 1X Ž S* . x J Ž Ec . s
f 1X Ž S* . x 0
yf 1 Ž S* .
ž
1y
x D2
f 2X Ž x .
/Ž
f 1 Ž S* . y D 1 .
x Ž f 1 Ž S* . y D 1 . D2
0
f 2X Ž x .
yD 2
.
0
Ž 3.7.
81
DISTINCT REMOVAL RATES
The eigenvalues of EŽ Ec . satisfy
3 q a1 2 q a2 q a3 s 0,
Ž 3.8.
where a1 s 1 q f 1X Ž S* . x q a2 s Ž 1 q f 1X Ž S* . x .
ž
ž
x D2
x D2
f 2X Ž x . y 1
f 2X Ž x . y 1
/Ž
/Ž
f 1 Ž S* . y D 1 . ,
f 1 Ž S* . y D 1 .
q x Ž f 1 Ž S* . y D 1 . f 2X Ž c . q x f 1 Ž S* . f 1X Ž S* . , and a3 s x f 2X Ž x . Ž 1 q f 1X Ž S* . x .Ž f 1 Ž S* . y D 1 . . Note that the constant term a3 is positive, so the Routh᎐Hurwitz criterion says that Ec will be asymptotically stable if and only if a1 ) 0 and a1 a2 ) a3 . We summarize the above results in the following theorem. THEOREM 3.1. If s ) 1, then only E1 exists and E1 is locally asymptotically stable. If s - 1 and Ž1 y s .rD 1 - x , then only E1 and E2 exist, E1 is unstable, and E2 is locally asymptotically stable. If s - 1 and Ž1 y s .rD 1 ) x , then E1 , E2 , and Ec exist, and E1 and E2 are unstable. Ec is locally asymptotically stable if a1 ) 0 and a1 a2 ) a3 . If D 1 s D 2 s 1, then the limit plane of Ž1.6. is ⌺ : S q x q y s 1. If we drop the S equation, then the limit system of Ž1.6. takes the form x⬘ s x Ž f 1 Ž 1 y x y y . y 1 . y f 2 Ž x . y, y⬘ s y Ž f 2 Ž x . y 1 . , x Ž 0 . ) 0,
Ž 3.9.
y Ž 0 . ) 0.
For this reduced system, if Ec exists then the Jacobian matrix at Ec is J Ž Ec . s
y x f 1X Ž S* . y Ž x f 2X Ž x . y 1 .Ž f 1 Ž S* . y 1 . y x f 1X Ž S* . y f 2 Ž x .
x Ž f 1 Ž S* . y 1 . f 2X Ž x .
0
.
Ž 3.10.
82
LI AND KUANG
The determinant of this matrix is w x f 1X Ž S*. q f 2 Ž x .xw x Ž f 1Ž S* y 1. which is positive, and the trace is y x f 1X Ž S*. y Ž x f 2X Ž x . y 1. Ž f 1Ž S*. y 1.. Therefore if y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f Ž S*. y 1. - 0 then Ec is locally asymptotically stable; if y x f 1X Ž S*. y Ž x f 2X Ž x . y 1. Ž f 1Ž S*. y 1. ) 0 then Ec is a repeller in ⌺, and in this case there is a periodic solution in ⌺ Žby an application of the Poincare᎐Bendixson ´ teorem.. We summarize these in the following theorem. f 2U Ž x .x,
THEOREM 3.2. Assume D 1 s D 2 s 1 and Ec exists Ž i.e., x q s - 1.. If y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. - 0 then Ec is locally asymptotically stable. If y x f 1X Ž S*. y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. ) 0 then Ec is unstable, and there is a periodic solution in the plane ⌺ : S q x q y s 1.
4. GLOBAL ANALYSIS In the previous section, we showed that if only E1 exists then E1 is asymptotically stable, if E1 and E2 exist then E1 is unstable and E2 is Žlocally. asymptotically stable, if E1 , E2 , and Ec exist then E1 and E2 are unstable, and in this case Ec may or may not be stable. We shall show that E1 is globally asymptotically stable if only E1 exists. The proof is very straightforward. Most importantly, we shall show that if only E1 and E2 exist, under a reasonable additional assumption E2 is globally asymptotically stable. The proof involves the construction of a Lyapunov function and the application of the Lyapunov᎐LaSalle theorem. ŽWe shall use Theorem 2.1 in Wolkowicz and Lu w13x, which is a slightly modified version of the statements given in LaSalle w6x and Hale w3x.. We shall also show that system Ž1.6. is uniformly persistent if Ec exists. The following theorem states that E1 is a global attractor if it is the only steady state Ži.e., s ) 1.. THEOREM 4.1.
If s ) 1, then all solutions of Ž1.6. satisfy lim Ž S Ž t . , x Ž t . , y Ž t . . s Ž 1, 0, 0 . .
tª⬁
Proof. Since S Ž t . - 1 for large t and f 1Ž1. y D 1 - 0 Ži.e., s ) 1., there is ␣ ) 0 such that x⬘Ž t . - y␣ x Ž t . for t sufficiently large. This shows lim t ª⬁ x Ž t . s 0. It follows from the third equation of Ž1.6. that lim t ª⬁ y Ž t . s 0. Then the first equation of Ž1.6. yields lim t ª⬁ SŽ t . s 1. The proof is complete. Note that if Ž s q x . Dmin ) 1 then s q D 1 x ) 1. That is, 1rDmin s q x implies Ž1 y s .rD 1 - x .
83
DISTINCT REMOVAL RATES
THEOREM 4.2.
If s - 1 and 1
- s q x ,
Dmin
Ž 4.1.
then all solutions of Ž1.6. satisfy lim Ž S Ž t . , x Ž t . , y Ž t . . s s ,
ž
tª⬁
1 y s D1
/
,0 .
Proof. We choose d1 ) Dmax and d 2 - Dmin such that 1 d2
- s q x
Ž 4.2.
and that, for large t, 1 d1
- SŽ t . q xŽ t . q yŽ t . -
1 d2
.
Ž 4.3.
Let
␣s1q
max
0FxF Ž1y s .rD 1
½
f 2 Ž x . Ž Ž 1 y s . rD 1 y x . x Ž D2 y f 2 Ž x . .
5
and
s1q
max
xFxF1rDmin q1
½
␣ Ž f 2 Ž x . y D2 . D2
5
.
Let C Ž u. be a continuously differentiable function and C⬘Ž u. is given by
¡0, C⬘ Ž u . s
~
if u F
 x q s y 1rd 2
¢ ,
ž
u q s y
1 d2
/
,
if
1 d2
1 d2
y s - u - x ,
u G x .
C⬘Ž u. is simply linear on w1rd 2 y s , x x. In view of Ž4.3., xqyF
1 d2
y s
if S G s .
y s ,
84
LI AND KUANG
Therefore, if S G s , then C⬘Ž x q y . s 0. The graph of C⬘Ž x q y . is shown in Fig. 1. Define the Lyapunov function V Ž S, x, y . as follows Vs
S
Ž f 1Ž . y D1 . Ž 1 y s .
s
D1Ž 1 y .
H
q x y x* y x*ln
ž
x x*
/
d
q ␣ y q CŽ x q y.
Ž 4.4.
on the set ⌿ s Ž S, x, y . : S g Ž0, 1., x, y g Ž0, q⬁., S q x q y g Ž1rd1 , 1rd 2 .4 , where x* ' Ž1 y s .rD 1. Then the time derivative of V along solutions of the differential equation is V˙ s C⬘ Ž x q y . q 1 y q
f2 Ž x . x
ž
1 y s D1
Ž 1 y s . f 1Ž S . Ž f 1Ž S . y D1 . x D1Ž 1 y S . y x q ␣ Ž f 2 Ž x . y D 2 . y D 2 C⬘ Ž x q y . y.
/
First, note that w1 y Ž1 y 2 . f 1Ž S .rD 1Ž1 y S .xŽ f 1Ž S . y D 1 . x is nonpositive for 0 - S - 1 and equals 0 for S g w0, 1. if and only if S s s or x s 0. Since C⬘Ž x q y . s 0 for S G s and C⬘Ž u. G 0 for u G 0, C⬘Ž x q y . Ž f 1Ž S . y D 1 . x is nonpositive for S g w0, 1.. Therefore the first term in V˙ is always nonpositive and equals 0 for S g w0, 1. if and only if S s s or x s 0.
FIG. 1. A graphical depiction of C⬘Ž x q y . for S g w0, 1..
85
DISTINCT REMOVAL RATES
Define h Ž S, x, y . s
f2 Ž x . x
ž
1 y s D1
y x q ␣ Ž f 2 Ž x . y D 2 . y D 2 C⬘ Ž x q y . .
/
If 0 - x F 1 y srD 1 , then f2 Ž x . x
ž
1 y s D1
yx G0
/
and
f 2 Ž x . y D2 F 0
and they are equal to 0 if and only if x s Ž1 y s .rD 1. Note that C⬘Ž x q y . is always nonnegative. By the definition of ␣ , hŽ S, x, y . - 0 for 0 F x - Ž1 y s .rD 1 and possibly hŽ S, x, y . s 0 if x s Ž1 y s .rD 1. If Ž1 y s .rD 1 - x F x , all three terms in hŽ S, x, y . are nonpositive and one can easily see hŽ S, x, y . - 0. If x ) x , then x q y ) x and therefore C⬘Ž x q y . s  . Note that, if x ) x , only the second term in hŽ S, x, y . is nonnegative. According to the definition of  , hŽ S, x, y . - 0 if x ) x . Therefore hŽ S, x, y . - 0 for x G 0 and x / Ž1 y s .rD 1 , and it is possibly 0 if x s Ž1 y s .rD 1. By Lemma 2.1 every bounded solution of Ž1.6. is contained in ⌿, and hence by Theorem 2.1 in w13x every solution of Ž1.6. approaches the set ⌳, the largest invariant subset M of ⌽ s Ž S, x, y . g ⌿ : V˙ s 04 . ⌽ is made up of points of the following forms where S g w 0, 1 x ,
Ž S, 0, 0 . ,
ž
s ,
1 y s D1
/
,y ,
Ž s , x, 0 . ,
where y g w 0, ⬁ . , where x g w 0, ⬁ . .
Since V is bounded above, any point of the form Ž S, 0, 0. cannot be in the 3 Ž -limit set ⍀ of any solution initiating in the interior of Rq . s , x, 0. g M Ž . Ž . implies that S t s , which in turn leads to 0 s S⬘ t s 1 y s y f 1Ž s . x and hence x s Ž1 y s .rD 1. Ž s , Ž1 y s .rD 1 , y . g M implies that S Ž t . s s and x Ž t . s Ž1 y s .rD 1. The second equation of Ž1.6. implies that x⬘ s 0, which yields y s 0. Therefore M s E2 4 . This completes the proof. THEOREM 4.3. If s - 1 and Ž1 y s .rD 1 ) x , then system Ž1.6. is uniformly persistent; i.e., there exists a constant ) 0, independent of initial conditions, such that lim inf S Ž t . G , tª⬁
lim inf x Ž t . G , tª⬁
lim inf y Ž t . G . tª⬁
86
LI AND KUANG
Proof. Choose
½ ½ ½
X 1 s Ž S, x, y . ; 0 F S F 1, 0 - x F Y1 s Ž S, x, 0 . ; 0 F S F 1, 0 F x F Y2 s Ž S, 0, y . ; 0 F S F
1 Dmin
1 Dmin 1 Dmin
q 1, 0 - y F
1 Dmin
5
q1 ,
5
q1 ,
q 1, 0 F y F
1 Dmin
5
q1 ,
and X 2 s Y1 j Y2 .
FIG. 2. m1 s 3.6, a1 s 0.8, D 1 s 1.4, m 2 s 3, a2 s 0.6, D 2 s 1.2. Ž S Ž0., x Ž0., y Ž0.. s Ž0.6, 0.2, 0.7.. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. In this case, s q x - 1. Clearly, the solution approaches the predator-free steady state.
DISTINCT REMOVAL RATES
87
Then X 1 and X 2 are two disjoint subsets of R 3 , X 2 is compact, X s X 1 j X 2 is also compact, and X 1 and X 2 are positively invariant for Ž1.6.. By Lemma 2.2, X 2 and X are global attractors in the union of the S y x 3 plane and S y y plane and in Rq , respectively. We prove that X 2 is a uniformly strong repeller for X 1 Žfor the definitions of a uniformly strong repeller as well as a weak repeller, see Thieme w11x.. E1 and E2 are the only steady states in X 2 . E1 is a saddle in R 3 and its stable manifold is Ž S, 0, y .; y G 04 . E2 is also a saddle in R 3 and its stable manifold is Ž S, x, 0.; x ) 04 . Therefore E1 and E2 are weak repellers for X 1. The stable manifold structures of E1 and E2 Ž E2 is a global attractor in the S y x plane. imply that they are not cyclically chained to each other on the boundary X 2 . By Proposition 1.2 of Thieme w11x, X 2 is a uniform strong repeller for X 1; that is, there are ␦ 1 ) 0 and ␦ 2 ) 0 such that lim inf t ª⬁ x Ž t . ) ␦ 1 and lim inf t ª⬁ y Ž t . ) ␦ 2 with ␦ 1 and ␦ 2 not depend-
FIG. 3. m1 s 3.6, a1 s 0.8, D 1 s 1.4, m 2 s 3, a2 s 0.6, D 2 s 1. Ž SŽ0., x Ž0., y Ž0.. s Ž0.6, 0.2, 0.7.. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. Clearly, the solution approaches a positive steady state.
88
LI AND KUANG
ing on the initial values in X 1. Applying Proposition 2.2 of Thieme w11x to the first equation of Ž1.6. yields that there is ␦ 3 ) 0 such that lim inf t ª⬁ SŽ t . ) ␦ 3 with ␦ 3 not depending on the initial values of X 1. This completes the proof.
5. DISCUSSION In this paper, we considered a food chain with one prey and one predator in the chemostat. In this model, the prey consumes the nutrient and the predator consumes the prey but the predator does not consume the nutrient. We assumed that the functional response functions are general monotone response functions and the removal rates are different. The model we considered is more general and realistic than the models in w1, 4, 9, 12x.
FIG. 4. m1 s 8.5, a1 s 0.6, D 1 s D 2 s 1, m 2 s 6, a2 s 0.6. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. The top curve depicts SŽ t ., the middle one depicts y Ž t ., and the bottom one depicts x Ž t .. The solution appears to approach a periodic solution.
DISTINCT REMOVAL RATES
89
The main difficulty we faced was the lack of a conservation principle which was lost due to the different removal rates. In the case of different removal rates, the system cannot be reduced to a two-dimensional system and we therefore must look at the full system. We found that the washout steady state E1 is the global attractor if it is the only steady state Žthis happens when s ) 1.. This confirms the intuition that both prey and predator cannot persist if the removal rate of the prey is relatively large. When E1 and the predator-free steady state E2 are the only steady states, we found that E1 is unstable and E2 is locally asymptotically stable. By constructing a Lyapunov function, we were able to show that if E1 and E2 are the only steady states, under an additional assumption Žsee Ž4.1.., E2 is a global attractor. The construction of the Lyapunov function is rather novel and nontrivial. This novel idea has been used in w7, 8x. This condition does not depend on the specific properties of the functional response functions, and it becomes necessary if Dmin is close to both D 1 and 1. The
FIG. 5. m1 s 8.5, a1 s 0.6, D 1 s 1.1, m 2 s 6, a2 s 0.6, D 2 s 1. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. The top curve depicts SŽ t ., the middle one depicts y Ž t ., and the bottom one depicts x Ž t .. The solution oscillates but eventually approaches a positive steady state.
90
LI AND KUANG
global stability of E2 implies that the predator will be washed out in the chemostat regardless of the initial density levels of prey and predator. We also showed that, when Ec exists, the prey and predator coexist in the sense that the system is uniformly persistent. In this case, a switch of the stability of the interior steady state Ec may occur. E2 is a global attractor if it is locally asymptotically stable and Ž4.1. holds. Based on our extensive simulation work where the functional response functions take the Michaelis᎐Menten form f 1Ž S . s
m1 S a1 q S
and
f2 Ž x . s
m2 x a2 q x
,
we conjecture that this steady state remains globally asymptotically stable as long as it is locally stable Žsee Fig. 2.. As E2 becomes unstable, a locally asymptotically stable interior steady state Ec bifurcates from it. Our simulation work ŽFig. 3. suggests that Ec is a global attractor if it is locally asymptotically stable. As certain parameters increase or decrease further away, Ec loses its stability and oscillatory solutions appear. These oscillatory solutions Žsee Figs. 4 and 6. appear to be the results of Hopf bifurcations. Figure 4 shows a case in which D 1 s D 2 s 1 and system Ž1.6. possesses periodic solutions. Figure 5 indicates that perturbing D 1 Žwhile keeping other parameters in Fig. 4 fixed. leads to a bifurcation. For example, changing D 1 s 1 in Fig. 4 to D 1 s 1.1 in Fig. 5 seems to destroy the periodic solutions and possibly leads to the global stability of Ec . Therefore varying the values of D 1 and D 2 may affect the dynamics of Ž1.6. in a very surprising and significant way. Next, we see how the parameters D 1 and D 2 affect the dynamics of Ž1.6. if S 0 is fixed. D 1 s 1 q 1 and D 2 s 1 q 2 , where 1 and 2 denote the scaled specific death rates of the prey and predator, respectively. ŽNote that the analysis of the model requires no assumptions on the sign of i ’s, as long as the Di ’s all remain positive. This leaves the IDi ’s open to other interpretations.. Assume that D 1 and D 2 are large enough so that s ) 1, then both prey and predator populations will be washed out Ž E1 is stable . in the chemostat. As D 1 is gradually decreased, eventually there is a bifurcation when s - 1 and Ž1 y s .rD 1 - x . In this case, E1 loses its stability and the new bifurcated steady state E2 is asymptotically stable. As D 2 is gradually decreased, the next bifurcation occurs when s - 1 and Ž1 y s .rD 1 ) x hold. In this case E2 loses its stability, and a new interior steady state Ec appears. If D 1 s D 2 s 1, the conservation principle holds; that is, the -limit sets of solutions of Ž1.6. lie in the plane ⌺ : S q x q y s 1. In this case, one can easily show that Ec Žif it exists. is locally asymptotically stable if and only if f 1X Ž S*. x q Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. ) 0. When this in-
DISTINCT REMOVAL RATES
91
FIG. 6. m1 s 4, a1 s 0.6, D 1 s 1.1, m 2 s 5, a2 s 0.5, D 2 s 1.2. Ž SŽ0., x Ž0., y Ž0.. s Ž0.1, 0.7, 0.8.. In this case, it can be shown that Ec exists and is unstable. The top curve depicts S Ž t ., the middle one depicts x Ž t ., and the bottom one depicts y Ž t .. The solution oscillates and seems to approach a periodic solution.
equality is reversed, Ec will be a repeller in ⌺, and there will be at least one periodic orbit Žby an application of the Poincare᎐Bendixson theorem.. ´ Determining the number of periodic solutions is a deep mathematical problem. Kuang w5x has shown in the case of Michaelis᎐Menten-type response functions, if yf 1X Ž S*. x y Ž x f 2X Ž x . y 1.Ž f 1Ž S*. y 1. is small and positive, then the limit cycle is unique and asymptotically stable. We guess this is true for general response functions. In w1x, it was shown that in the case of Michaelis᎐Menten-type response functions and D 1 s D 2 s 1, Ec is globally asymptotically stable if it is locally asymptotically stable. It remains open if this is true in the case of general response functions and different removal rates.
92
LI AND KUANG
REFERENCES 1. G. J. Butler, S. B. Hsu, and P. Waltman, Coexistence of competing predators in a chemostat, J. Math. Biol. 17 Ž1983., 133᎐151. 2. J. F. Drake and H. M. Tsuchiya, Predation of escherichia coli by copoda steuii, Appl. En¨ iron. Microbiol. 31 Ž1976., 870᎐874. 3. J. K. Hale, ‘‘Ordinary Differential Equations,’’ Krieger, Malabar, FL, 1980. 4. J. L. Jost, S. F. Drake, A. G. Fredrickson, and M. Tsuchiya, Interaction of tetrahymena pyriformis, escherichia coli, azotobacter vinelandii and glucose in a minimal medium, J. Bacteriol. 113 Ž1976., 834᎐840. 5. Y. Kuang, Limit cycles in a chemostat-related models, SIAM J. Appl. Math. 49 Ž1989., 1759᎐1767. 6. J. LaSalle, Some extensions of Lyapunov’s second method, IRE Trans. Circuit CT-7 Ž1960., 520᎐527. 7. B. Li, Global asymptotic behavior of the chemostat: General response functions and different removal rates, SIAM J. Appl. Math. 59 Ž1999., 411᎐422. 8. B. Li, ‘‘Analysis of Chemostat-Related Models with Distinct Removal Rates,’’ Ph.D. thesis, Arizona State University, 1998. 9. G. Sell, What is a dynamical system? in ‘‘Studies in Ordinary Differential Equations’’ ŽJ. Hale, Ed.., MAA Studies in Mathematics, Vol. 14, Math. Assoc. America, Washington, DC, 1977. 10. H. L. Smith and P. Waltman, ‘‘The Theory of the Chemostat,’’ Cambridge Univ. Press, Cambridge, UK, 1995. 11. H. R. Thieme, Persistence under relaxed point-dissipativity with an application to an epidemic model., SIAM J. Math. Anal. 24 Ž1993., 407᎐435. 12. H. M. Tsuchiya, S. F. Drake, J. L. Jost, and A. G. Fredrickson, Predator᎐prey interactions of dictyostelium discordeum and escherichia coli in continuous culture, J. Bacteriol. 110 Ž1972., 1147᎐1153. 13. G. S. K. Wolkowicz and Z. Lu, Global dynamics of a mathematical model of competition in the chemostat: General response function and differential death rates, SIAM J. Appl. Math. 52 Ž1992., 222᎐233.