Global Stability in a Synthrophic Chain Model A. S. ELKHADER Department of Mathematics and Natural Sciences, Northern State Uniuersity, Aberdeen, South Dakota 57401 Receiced I4 March 1990; recised 5 November 1990
ABSTRACT A mathematical model formulated by G. Powell in 1986 is considered. This model describes a synthrophic chain of species X,, 1~ i G n, in a continuous chemostat culture. Species X, utilizes substrate 5, and forms a product S,,,. Substrate S,+t inhibits the growth of X, and constitutes a growth-limiting substrate of X,,,. For this model, conditions for coexistence of all species are improved, and the local stability results extended to global stability. Some interesting cases are discussed.
1.
the are
INTRODUCTION
In this paper, we study the behavior of a mixed culture, that is, a mixture of different types of organisms. The interaction between species in a mixed culture can best be tested by means of a continuous culture using the chemostat technique, which has an extensive description in the literature (see, for example, Pirt [5]). Because of the complexity of mixed-culture behavior, mathematical models of various systems are particularly valuable for describing and predicting this behavior. A mathematical model describing the dynamics of a particular class of mixed culture was introduced by Powell [6-81. This class is a synthrophic association of two species or more. Following Powell [8], a synthrophic chain of n species Xi, 1~ i < n, and n + 1 substrates Si, 1~ i < n + 1, in continuous culture is of the form
In this synthrophic chain, the growth of species Xi depends on the substrates Si and Si+i. Although the growth of X, may depend on the presence of many substrates, it is the levels of Si and Sj+r that are assumed to be growth-limiting. Species Xi utilizes substrate Sj and forms a product This product inhibits the growth of species Xi but constitutes a si+l. growth-limiting substrate for species X,+i. In this manner, an incremental increase in species Xi leads to more favorable conditions for the growth of MATHEMATICAL
BIOSCIENCES
OElsevier Science Publishing 655 Avenue of the Americas,
104:203-245
(1991)
Co., Inc., 1991 New York, NY 10010
203 0025-5564/91/$03.50
204
A. S. ELKHADER
species Xi_, and Xi+,. A mixed culture of the two species Pseudomonas and Hyphomicrobium is a typical example of this kind of association. The Pseudomonas species utilizes methane and produces methanol, which inhibits its growth. The second species, Hyphomicrobium, utilizes the methanol produced but not methane, thus overcoming the inhibition of the first species and allowing both species to grow. This cooperative interaction between the species will be the basic idea of the analysis in this work. To be precise, a mathematical model of 2n ordinary differential equations, formulated by Powell [8] to describe the dynamics of n species and n substrates in a chemostat culture, can be reduced to a tridiagonal cooperative system of ordinary differential equations. Cooperative systems are defined in the sense of Hirsch [2, 31. This system will be analyzed, and the local stability results will be extended to global stability, which will be the focus of this work. This work is organized as follows: In Section 2, we introduce the notation and discuss some results essential to the understanding of later sections. In Section 3, we describe the model formulated by Powell [S] and modify the function pn. Some of the previous results are presented. Also, we give some basic properties of solutions of the system and show that the system can be reduced to an n-dimensional system of the tridiagonal cooperative type. This structure, a result due to Smillie [9], and a result due to Markus [4] are the main tools in determining the global stability of the systems of differential equations in Sections 4 and 5. Section 4 is devoted to analyzing the system when only substrate S, is supplied to the culture and the function CL,,depends on substrates S, and system starts by studying Sn+l. The analysis of the reduced n-dimensional the flow on the boundary. This leads to establishing necessary and sufficient conditions for the existence of all boundary equilibria and a unique interior equilibrium. Arguments based on the implicit function theorem in R2 are the main tools in determining these conditions. Simple arguments based on the properties of cooperative matrices are used to study local stability of all equilibria. Description of the stable manifolds of all boundary equilibria and the structure of the variational matrix at the interior equilibrium allow us to apply the result of Smillie [9], which is that the interior equilibrium point is globally asymptotically stable. To establish the global stability of the (2n + l&dimensional system, we describe all possible equilibria and analyze their local stability. Arguments based on the Markus theorem [4] are used to determine the global stability of the system. The special case when the function pL, depends on s, only is studied. The system in [81, when n is set to equal 2 and p2 depends on S, and S,, is considered in Section 5. The analysis of this five-dimensional system is organized as in Section 4. We show that the existence of an interior equilibrium point is possible in three different ways. This improves the conditions given in [8]. The local stability analysis and the global stability
GLOBAL STABILITY IN A CHAIN MODEL
205
analysis of all possible equilibria are given. This section is concluded considering the special case when the function pz depends on s2 only. The paper ends with a brief discussion. 2.
by
PRELIMINARIES
In this section we introduce some notations and nonnegative matrices and cooperative systems that later sections. Also, a result on the implicit function given. Following Hirsch [2, 31, we say that the system of i = F(x),
XER”,
F=
state some results on are essential tools in theorem in [w2 will be equations
(f,,...,fn>
(2.1)
is cooperative (competitive) if afi /Jxj > 0 (af, /axj < 0 for i # j, tridiagonal if afi /ax, = 0 for ]i - j] > 1, and strongly cooperative tridiagonal if af, /ax, is positive for Ii - j] = 1. Also, an it X n matrix A is called a stable matrix if S(A) < 0, where S(A) = maxRe A and A runs over the eigenvalues of A. We let x,):
rwT={(x,,...,
xi > 0, 1 < i < n) ;
k!: is the interior of rWF containing the The next theorem, which will be used condition for determining the stability off-diagonal elements. It is part of a matrices in Smith [lo, 111.
positive vectors. repeatedly, provides a very simple of a matrix A with nonnegative general theorem on nonnegative
THEOREM2.1 Zf A is an n x n matrix with nonnegatiue off-diagonal elements, then S(A) < 0 is equivalent to the condition that the principal minors of A alternate in sign. That is,
A direct consequence of this theorem is the following proof is given in Appendix A. Let M(n) be -l(l+K,)
whose
K,
-(l+K)
. . ..f.............2....i.......................
K,
) -Cl+
Ki > 0, 1 Q i < n.
lemma,
K-d 1
Kn-, -(l+K,)
206
A. S. ELKHADER
LEMMA 2.2
The matrix M(n)
is a stable matrix.
(See Appendix
A.) The next theorem, due to Smillie [9], is one of the important tools in determining the asymptotic behavior of the solution of the tridiagonal cooperative systems, which will appear in later sections. THEOREM 2.3 [Smilliel
Let j, = F(x) be a strongly cooperative tridiagonal system of ordinary differential equations defined on R c R”. Assume that the component functions F, are n - 1 times differentiable. Let cp(t> be a solution of d(t) = F(cp), defined on a maximal interval of the form [0, a), 0 < a 0
t-m
for all i
whenever xi(O) > 0 for all i. DEFINITION
System (2.1) is uniformly persistent if there exists 6 > 0 such that whenever xi(O) > 0 for all i, then
lim infxi(t)>6>0 t+m
for all i
Finally, the following lemma, which is a consequence of the implicit function theorem in R2, will be used repeatedly. The lemma is illustrated in Figure 1. LEMMA 2.4
Let cp and IJ be continuously differentiable functions
on 0 f x < a. As-
sume that &a> = a, x < cp(x) for 0 < x < a, dp(x)/dx < 1, 0 < $(O), and $(x) is a linear function with positive slope. Let p = p(r,, r2) be a continu-
GLOBAL
STABILITY
IN A CHAIN
207
MODEL
Y
(ON
(W
(l?O)
FIG. 1. The function
x = q(y)
as described
in Lemma
3.4.
ously differentiable function on iwl, and let
For (x, y) E D, let r, = q(x)-
F(x,Y)
x and r2 = $,(x)-
=p(r1(x),r2(x,y)) -1
y. Define onD.
_ Assume that $L /Jr, > 0 and ap /ar, < 0. If F(x, y) = 0 at (.?,O) and (x, y), where 0 < .C < X < a and j = 1,9(x>, then there exists a unique smooth function x = v(y), 0 < y < J, such that F(T$ y), y) = 0, $0) = 2, and $ji) = X. Moreouer, dq/dy < l/(d$/dx).
208 Proof.
regions
A. S. ELKHADER
Let
The function in x because
Similarly,
@L(x)= I/J(O)+ (YX, LY> 0. Divide
the region
D
into
two
D, and D,, where
F(x, y) is continuously
F(x, y) is increasing
differentiable
on D and is decreasing
in y because
C?F
-=$+)>O. 8Y
Now consider the lines in D, of the form y = 0, 0 < 8 < $(O>. The function F(x, y) evaluated along each line y = 0 is given by
F(x,B)=~(cp(x)-x,lCI(~)-e)-l =/.&(cp(x)-x,$(O)+ax-8)-l. Since p’(x) < 1, F(x,B)
is strictly decreasing
in x for each 0. Moreover,
F(~,B)=II(~(~)--,IC~(~)--)-~
>p(cp(l)-i,cL(-i))-l = F(i,O) = 0, and F(x,B)=p(cp(Z)-x,$(-3:)-0)-1 = 0. Observe that ~(0) = 2, and let 2, = 7($(O)). The implicit function theorem guarantees that 77 is continuously differentiable. Next, consider the lines y = $,(O)+ ix, 0 ( 6 < (Y. These lines lie in D,. The function F(x, y) evaluated along any of those lines, for a fixed e, is
209
GLOBAL STABILITY IN A CHAIN MODEL
given by F( x,$(O) +
ix) = /A(q(x) - x,$(O) + ax -$(O) - w-1 =&+7(x)-X&X-$)x)-l.
It is easily seen that the function
F is decreasing
and satisfies
F(~~,~(O)+B~,)=CL('P(~~)-~~,(~--)~~)-~ > /_L( cp(&J -
io, “1,) -
1 = 0,
F(i,$(O)+i&)=&(X)-.?,(a-8)3)-l in D implies that x = q(y). Moreover,
dv
dy=
-$-1)
- c3F/ay t’F/ax
where k =
(aF/c+,(dd~ aF/ar,
-1))
> 0.
n
Remarks. (1) If, in Lemma 2.4, the condition F(_f,O) = 0 is replaced by F(0, F), 0 < jj 0.
(3.1)
Here, x0 = 0, xi denotes the cell density of species X,, si denotes the concentration of substrate Si, so is the concentration of substrate Si entering the chemostat, x (g/mole) denotes the growth yield of species X, with respect to its substrate Si, y.,i_ 1 (mole/g) denotes the product yield of substrate Si by the species Xi-i, and pi denotes the specific growth rate of species X, and depends on the concentrations si and si+, of substrates Si and Si+i, respectively. D is the dilution rate. As usual, the dot denotes the time derivative. The only difference between system (3.1) and the model in [8] is in the growth function p,. Here, CL,,depends on Sri+++in the same fashion pi depends on Si+i. Thus, we include the dynamics on the substrate SE+,, which is assumed to be zero in [8]. In order to characterize the interaction sufficiently without postulating particular functional forms for the specific growth rate functions, a number of qualitative assumptions have to be made. For convenience these are considered collectively: (Al) The specific growth rate of species concentrations si and si+i, pLi= I_L~(s~,s~+~), aPi(si,si+l) asi
>
aWi(si~si+l)
0
3
$+1
Xi depends
B pyax. Typically, the maximum specific growth rate and saturating /4= is attained in the limit of zero product concentration substrate concentration. (A31 In the absence of substrates Si and Si+,, the growth function pi is zero, that is, ~~(0,0)= 0. Substrate S,+i is derived solely from the growth of species X,. Some of the specific forms of the function pi frequently used (see Powell [8]; Pirt [S]) are
s -(W,/W,)p = Pmax s + K, +(k,/k,)p
P(S,P)
’
where k,, k,, k,, k,, k,, and k, are constants. The special case of system (3.1) when II = 2, with the assumptions that S, is solely derived from the growth of species X, (i.e., si = 0) and that the specific growth function pZ depends on S, only [i.e., S, is removed from the chemostat or ignored; pz = p&l] was analyzed by Powell [8]. It was shown that when the dilution rate D is less than the critical dilution rates for each species X, and X, (i.e., D < min{~;l”,~~ax}), a locally asymptotically stable steady state exists. System (3.1) with the assumption that Sn+l = 0 was addressed by Powell [8]. It has been shown that when the dilution rate D is less than the critical dilution rates of each species Xi, i=l , . . . , n (D < (ELF;“,pLr;l=,. . . , pzau}), a locally asymptotically stable equilibrium point with all species nonzero exists. Some of the previous results will be improved, and the global stability question in both cases will be answered. The special case when p,, depends on S, only will be discussed. 3.2.
SCALING
Before we start to analyze the above cases, it will be more convenient to reduce the number of independent parameters by scaling the system (3.1). This can be done if we perform the following substitutions: Ei = xi/kj,
si = Si/li, ai-
= y,i_,Ki_,
i;ii =
/liT
Pi(Fi,si+l)
D
’
ao=o, (3.2)
A. S. ELKHADER
212
where I, =$>O, ki=xli, for each i, 1 o mii+, = xil 2 ( asi+l i +i
m,l_l = xjz
Matrix M,, is of size i x(n - i) with (LV,,)~, = elsewhere. Matrix M,, is of size (n - i) x (n - i) except (A&J, =~~+~(xj,O)-l> 0. Matrix 0 is eigenvalues of M(E,) are the eigenvalues of tridiagonal cooperative matrix M,,, let kj = and
-Gj/asj+l
ai+ias, >o
xjapi /as2 > 0 and zero with - 1 in the diagonal of size (n-i>Xi. The M,, and M,,. For the
222
A. S. ELKHADER
For 1 d j d i, let Cj be the jth principal
-(1+kJ 1 cj = lxj
minor of M,,,
k, -(l+G . .
k, . .
. . 1
-(l+
kj)
The determinant on the right-hand side is the jth principal minor of the stable matrix described by Lemma 2.2. Since aj > 0, the signs of the Cj’s are the same as the signs of the principal minors of the stable matrix of Lemma 2.2. Hence, M,, is a stable matrix. The eigenvalues of M,, are (M,,),, > 0 and the n -(i + 1) diagonal entries mkk = - 1, i + 2 < k < n. Thus, Ei is a saddle point, stable in the xj direction j f i + 1, and unstable in the xifl direction. Finally, the variational matrix evaluated at E,, M(E,), is a tridiagonal cooperative matrix with aPi
mii_, = x:-g
> 0, 1
mi,i = x:
-
-aPi dSi
i
miifl = xi nG i
-
+
‘Pi dsi+l
O.
An analysis similar to that followed for the matrix M,, shows that ME,) a stable matrix and therefore E, is locally asymptotically stable.
is
Global Stability Analysis. To determine the global stability of the reduced system (4.3), we describe the stable manifolds of all possible equilibria. Let M+(Ei) be the stable manifold of Ei. PROPOSITION 4.2
For system (4.3), assume that E, exists. Then for each i, 0 < i < n,
Proof: We prove the proposition by induction. Let Bi = {x E H,‘: xi > 0 x,(t)) be the solution of system if and only if j 0, then x,(t) + 0 for all i. But since ~~(1~0) > 1, it is clear from system (4.3) that x,(t) does not tend to 0 unless x,(t) is identically 0. In . . particular, x y = 0, a contradtction.
GLOBAL STABILITY IN A CHAIN MODEL
Next, suppose that M+(E,)
223
= B,, 1~ i. We want to show that
M+(Ei+J = 4+1. To do so, let x0 E Bi+i, and hence xz = 0 for k > i + 2. The solution x(t) of this system tends to an equilibrium by Smillie’s theorem, Theorem 2.3, as t tends to infinity. This equilibrium must be one of E,, E,, . . . , Ei+ ,. But if this equilibrium is not Ei+,, then we have a contradiction to the hypothesis, since xi”+, > 0 and x0 is not on the M+(E,) for some I ,xi”+i > 0; then the trajectory through this point must go to the equilibrium (cpi 0 . . . 0 cpi(.Ci+i), . . . , fi+i, 0 , . . .,O)as t goesto infinity. Since x,+,(t) -+ Zi+, and x,(t) + 0, r > i +2, it follows that ~,+~(Xi+~-Xi+2,Xi+2-Xi+3)~~i+2(~i+~,o)as t+w. Since /Li+2Cfi+l,0) > 1, it is clear from system (4.3) that xi+2 + 0 is impossible unless xi+2 is identically zero. In particular, xi”,, = 0, and thus xi = 0 for k > i +2 since x0 E H,‘. Hence, M+(Ei+,)= Bi+l. W A corollary
to this proposition
is the following.
COROLLARY4.3
For system (4.3), assume that E, exists. Then E,, attracts all solutions with positive initial conditions in H,‘. The final result in this section is the global stability of system (4.2). To show that we let si, 0 < i 6 n, denote the equilibria of system (4.2). The equilibria in (x; v> coordinates are given by
8,=(0
0;o
)...)
8,=(&O
)...)
0)
)...)
0;o
,...)
~2=(qo1(~2),i*,o )...)
The variational
0) 0;o
,...)
0)
matrix of system (4.2) is of the form
where A is equal to M, the 12x n variational matrix of system (4.31, B is an n x(n + 1) matrix with bii = +L~/~v~, bi,i+i =$.L~/~v,+,, and bij = 0 for
224
A.S.ELKHADER
C is an (n + 1) X (n + 1) diagonal matrix with diagonal entries - 1, and 0 is the (n + 1) X n zero matrix. The eigenvalues of the variational matrix of system (4.2) at f?; are equal to the eigenvalues of M(E,) and the negative eigenvalues of the diagonal matrix C. This implies that ii is stable if and only if Ei is stable. The persistence of solutions and the global stability of I?,, are given by the next two theorems. The main tool in the proof of these theorems is a theorem due to Markus [4]. The statement of the theorem of Markus and the proof of the next two theorems are given in Appendix B.
j # i, i + 1. Matrix
THEOREM
4.4
For system (4.2), assume that k,, exists. Then system (4.2) persists; that is, for all i, 1 0 as t -+m for every solution with initial condition (x(O); ~(0)) E Wn’, xi(O) > 0 for all i.
Finally, the global stability of 8, is given by the following theorem. THEOREM
4.5
Under the hypothesis of Theorem Special Case.
or that
its value
4.4, .&, is globally asymptotically stable.
Assume that substrate S,,+, is removed from the culture is very small so it can be neglected. Mathematically
FIG.4. The effect of substrate
S,,
, on X,.
GLOBAL
STABILITY
IN A CHAIN
MODEL
225
speaking, the specific growth function pn depends on s, only; that is, pn is a function of one variable. In this case, in the omega limit set, s, = x,-i x,. Let ji, = p,,(x,_ 1- x,). It is clear that the only term affected by the removal of s,+l is the function pn. This suggests that an interior equilibrium point exists if and only if p,(l,O) > 1,
All the results in this section hold in this case. Furthermore, since species X,, has no inhibitor, its density at the equilibrium point Z,* is larger than .?‘,* when S,, , is present. A geometric illustration is given in Figure 4, where the slope of the function CL,,is less than the slope of the function fin, and therefore the density of the X, species at equilibrium is greater than that when Sri+++is present. 5.
5.1.
ANALYSIS OF A CHAIN OF TWO SPECIES AND THREE SUBSTRATES WHEN S, AND S, ARE SUPPLIED TO THE CULTURE DESCRIPTION
In this section, our main goal is to study system (3.2) when II is set equal to 2. In this case, the chemostat culture consists of two species, X, and X,, and three substrates, S,, S,, and S,. The scaled system of ordinary differential equations to be analyzed is given by
Sj(0)
~ O
for j= 1,2,3.
a, is a positive constant. As in Section 3, let S1(XI,UI) = 1+ ut - xt, s2(x,,x2,u2) s3(x2,uJ
= 1+ c’2+ (Ytx, - x2, =
u,+
x2.
(5.1)
226 Expressing
A. S.
ELKHADER
pI and pZ in terms of x,, x2, ur, u2, and u3 gives
Thus, system (5.1) is equivalent
to
~.1=x,[~~(x,,x2,uI,u2)--l],
~2=x2[~2(x*,x2,U~,U*,U3)-l],
l,,
GLOBAL STABILITY IN A CHAIN MODEL
where it is the solution
of pt(l-
227
x1,1+ atxt)
= 1;
E,=(%f,), where _E, is the solution
of pz(l-
x2,x,> = 1; and
x1, l+ ~ytxt - x,) = 1 and pL2(1+ where is the solution of pt(lx2, x2) = 1. We proceed to find conditions under which such equilibalxlria exist. First, let
The simplest equilibrium point is E, = (O,O>. This equilibrium Next we consider one-species equilibria.
always exists.
(a) Suppose x2 = 0. The function F,(x,,O) = ~r(l- x1, l+ a,~,)-- 1 is strictly decreasing. Its maximum value is attained at x1 = 0, and its minimum at x1 = 1. Therefore, F,(l,O) f F,(x,,O) d F&0,0). Since F,(l,O) = prCO,l+ a,>- 1 < 0 because p,(O,O) < 0, then F,(x,,O) = 0 has a nonnegative solution if and only if F,(O,O) > 0. If F,(O,O) = ~t(l,l)1 = 0, then x1 = 0 is the only solution of F&x1,0) = 0. If F,(O,O) = p(l,l> - 1 > 0, then there exists a unique it, 0 < .?t < 1, such that F,(i:,,O)= 0. Let Et = (,Ct,O). Then E, exists if and only if p,(l,l) > 1. (b) Suppose x1 = 0. The function F&O, x,) = p&lx2, x,)- 1 is strictly decreasing, so it attains its maximum value at x2 = 0 and its minimum at x2 = 1. This implies that F,(O, 1) < FJO, x2) < &CO, 0). But F,(O, 1) = ~~(0,1) - 1 < 0 because pJO,O) < 0. Then FJO, x,) = 0 has a nonnegative solution if and only if F-JO, 0) > 0. If F,(O, 0) = 0, that is, p#,O) - 1 = 0, then F,(O, x,) = 0 has only the trivial solution. If F,(O,O) > 0, that is, p2(1,0) > 1, then &(0,x,) = 0 has a unique solution .C2, 0 < 1, < 1, such that F,(O,X,) = 0. Let E, =(O,f,). Hence, E, exists if and only if &(l,O) > I. Finally, we consider two-species algebraic system of equations
F,(x,,x,) = 0, has a nonnegative solution Al-A3, on the functions it
equilibria.
This case occurs when the
F,(Xl,X,) = 0
in H;. From the qualitative and CL*,the partial derivatives
(5.4) assumptions, aF, /ax1 < 0,
A. S. ELKHADER
228
aF,/ax,> 0,dF,/dx,> 0,and aF,/dx,and FJx,,x2)attain their maximum and minimum values at the boundary of H;. We start with F,(x,,xz)=O. Since dF,/ax, 0,then minF,(x,,x2)=F,(1,0)=~I(0,1+~1)-1 1, then there exists a nontrivial solution of F,= 0 in H:. Furthermore, this solution lies on a curve of the form x1 = cpl(x,) such that F,(cp,(x& x,)= 0. This will be realized by finding solutions of F,= 0 at the boundary of Hz and then applying Lemma 2.4 as follows. We evaluate Fl(x,, x,)along the boundary lines. Along x1= 1, F,(l, Along
x2) = /_~(0,1+ qx2)
- 1 < 0.
x2 = 0,
As above, F,(x,,O)= 0 at (Z:,,O) if and only if ~Jl,l)> Along x1 = 0,
1.
This is an increasing function; its maximum its F,(O, 1) = pJl,O)1 > 0. Thus F,(O, x,)= 0 has a nonnegative solution if and only if min F,(O, x2)= F,(O,O) = pl(l,l)1 < 0. If p,(l,l) = 1, then (0,O) is the only solution of F&O,x2)= 0. If ~~(1~1) < 1, then there exists XZ, 0 < X2 < 1, such that
F,(O,i,)= 0. Fmally, along x2 = 1+ (Y~x,, F,(x,,l+qx,)=/.q(l-x1,0)-1. This is a decreasing
function,
with
minF,=F,(l,l+cu,)=~(O,O)-l 1, then we need to consider the following cases: (1) 1 , ,F, < xl. Thus, by the remark following Lemma 2.4, there exists a smooth function xt = ‘pi(x,), 0 = (pi(X& that satisfies xi = P&X:), and FI(9i(x2), x2) = 0. F2(x1, x2) = 0. Since dF2 /ax, P- 0 and dF2 /ax, < 0, then
Next, consider
minF,(x,,x,) H:
= F2(0,1)
=/+(0,1)-l
This implies that FJx,, x2> = 0 has a nonnegative
~0.
solution
if
maxF,(x,,x2)=F,(1,0)=/L,(1+a,,0)-120. H: If F,(l,O) = 0, then (1,O) is the only solution of F2(x,, x2> = 0. If F,(l,O) > 0, that is, p2(1+ (Y,,0) > 1, then there exists a nonnegative solution of F2 = 0 in Hz. Furthermore, F2 = 0 along a curve of the form x2 = ~&xi). As before, this will be shown by Lemma 2.4 as follows. According to Remark 2 following Lemma 2.4, we need to show that F, = 0 along the boundary of H;. Along xi = 1,
&(I, 4 This function
= ,q(l+
is a decreasing
and a maximum
al-
~2,
function
~2) - 1,
0
0.
Therefore,
a unique
solution
x,, 0 < xz < 1 + (it, exists such that F~(x,,Q)
at (1, X,X
= 0
GLOBAL
STABILITY
IN A CHAIN
231
MODEL
Along x2 = 0,
This is an increasing
function
with
Thus Fz(xl, 0) = 0 has a solution
if and only if
minF,(x,,O)=F,(O,O)=pL(l,O)-ll. If pJl,O) > 1, then there exists an Z2, 0 < T2 < 1, such that F2 = (XI, x2>= 0 at (0, JZ~). Finally, along x2 = 1+ (~rx,, F,( x,, 1+ cqxl)
= p2(0,1 + (~rx,) - 1 < 0.
This analysis shows that if ~~(1 + (~r,0) > 1, then F, = 0 along three different boundary lines, and therefore we consider the following cases. (1) 1 < p&1,0) < pJl+ a,,O) (see Figure 6a). In this case F2(x1,x2) = 0 at (1,x,) and (X,,O). By the remark following Lemma 2.4, there exists a unique smooth function x2 = (p&x,) that satisfies
0 =4X,), and
X2=p2(1),
+CY,,
A. S. ELKHADER
232 x2
x2
4
I
5’ n (17F22)
(1,522)
I
Xl
a. 1 < pz(l,O)
x1
< p2(1 + W,O)
b. p&,0)
c. clz(l, 0) < 1 < P2(1 t FIG. 6. The isocline
-
Xl
(,g
= 1 < p2(1 t QI,O)
0170)
F2(xlr x,) = 0 as in the discussion
of solutions
of system (5.4).
4 (2) l+p2(1,0) = 0. Substituting xi = (P(Q) in F*, we obtain F,(cpd4>x,)
Since dye, /dyz < l/cri,
Therefore,
=/12(1+~1(~2)-
then F2(~i(x2),
x2>+>
x2) is a decreasing
F2((p,(x2), x,) = 0 has a solution maxk(cp,(x,),x,)
-1.
function
with
if and only if
=&(‘~i(%),%)
z=O.
If ~,(l, 1) < 1, then X2 > 0 and cpi(X,) = 0. Thus, maxF2(cp2(x2),x2)
=~~(1+~,4+(%)-
%,&)
-1
=y~(1-~-,,x,)-1 $2
x2
Cc)
(44)
(4
~~~~ 7. The isocline F, = 0 and F2 = 0 as described in Case 1 in Section 5. (a) ~JI,I)< 1 1, then there exists a unique function x2 = (p2(x1) such that
Any solution of system (5.4) must lie on this curve. Therefore, x2 = (p2(x1), 0 G x, Q 1, into F,(x,, x,1, we obtain
4(wP*(XJ) =/-+(l-
substituting
Xl,l+Qlxl-(P2(X1))-1
Again, this is a strictly decreasing function. Since its minimum is F,(l, (~~(1)) = &O, 1 + (Y,- x,1-- 1 < 0, then F,(x,, cp,(x,)) = 0 has a nonnegative solution if and only if max F,(x,, (p2(x,)) > 0, that is,
If /*i(l, 1 - (~~(01)= 1, then xi = 0 is the only solution of F,(x,, (pJx,l) = 0. If ~,(l,lcpJ0)) > 1, then there exists a unique XI*, 0 < x;” < 1, such that F,(x,*, cp,(x$>> = 0. Hence E, = 1
and
cL~(l,l~~l~cL1(l~l-(P2~O~~~
where .Lz = ‘p2(0) and XT = (pJxf>. Case 4. ~i(l,l)> 1 and &,O) > 1 (see Figure 819. In this case, j.~i(l,l) > 1 implies cpi(O)= 1, > 0, and since ~J1,0) > 1, then /+(l+ (~i(pi(O),O) > 1. Thus the conditions for the existence of an interior equilibrium E, in Case 2 are satisfied. Similarly, the conditions for the existence of E, in Case 3 are satisfied because ~~(1,0> > 1 implies (~~(0) = Z2 > 0 and p,(l,l) > 1 implies pi(l,l(pJO>)> 1. Hence, E, exists uniquely if and only if ~i(l,l) > 1 and ~J1,01> 1. Now let 2 be the set of all possible equilibria of system (5.3). Then the above analysis concludes the proof of the following theorem.
GLOBAL
STABILITY
THEOREM
IN A CHAIN
237
MODEL
5. I
For system (5.3), the following hold: (a> Z = (Es} if and only if ~~(1, 1) < 1, /.&,O) < 1. (b) C = {E,, E,} if and only if pr(l, 1) > 1, ~,(1,0> < 1. (c) Z = {E,, I&) if and only if p,(l,l) d 1, pJ1,O) > 1. (d) C =(E,,E,,E,) if and onfy if p#,l)> 1, P&~,O)G 1 1. Here, x2 = (p2(x1) satisfies Fz(xl, x2) = 0 and ‘pz(O)= i2. (f) C = (E,, E,, E,, EJ if and only if pJl,l) > 1 und PL~(~,O)> 1. Local Stability Analysis. Since we are interested in coexistence between species, we study only the cases in which an interior equilibrium exists. We have shown that coexistence between species X, and X, is possible in three different ways. We will study the local stability of all possible equilibria of system (5.3) in each case. Let (XT, x$> be the interior equilibrium point. The Jacobian (variational) matrix of system (5.3) evaluated at an equilibrium point E is of the form
where m 11
=
(
--al+ -+ar,s
aI+
dS1
x,+~.,(1-~x,,l+CY,x,-xX2)-l,
2 1
-a/-h
m12 =TX1’ m21 = m 22
=
+2 “la~,X27
(
-ap2
-
as,
+%as,
1
x
2
+/.l(l+a 2
The three cases of C to be considered Case 1.
1x 1
> 1
and
P2(LO)
where .Cr = p,,(O) > 0. The variational
2,x2>
-1.
are the following.
C = {E,, E,, E3}. This case happens
/4(l,L)
-x
0 and ~J1,0)1 < 0. Ea is a saddle point, unstable in the x1 direction and stable in the x2 direction. The equilibrium E, has M(E,) as a variational matrix,
The eigenvalues
of M(E,)
are /_LCLZ(lf(Yr~,,O)-l >o.
and Thus, Et is unstable in the x2 direction Similarly, at E,,
and stable in the x1 direction.
@l
-“?-
+“‘as, ME,) =
** alx2 as,
\ The off-diagonal are
\
aP1
elements
of this matrix are positive.
Its principal
minors
and
ah ap2 ah ap2 * *!?L%_ x*x*--+cYx*x*-->O detWE3)=~lx2 as,as2 1 2 as,as, 1 1 2 as,as, because aP2 /as, < 0. Then Theorem asymptotically stable steady state.
2.1 implies
that
Es is a locally
Case 2. 2 ={E,, E,,E3}.The conditions under which this can happen are pr(l,l) < 1 < p&1,1 - ‘p2(0)) and &,O) > 1, where P,(O) = f2 > 0. The variational matrix evaluated at E, takes the same form as in Case 1 with eigenvalues p&1,1)- 1 < 0 and p2(1,0) 2 1 > 0. Hence E, is a saddle point, stable in the x, direction and unstable in the x2 direction. At E,, M(E,) takes the form 0
/.&(1,1-&-l
WE21
=
1
G2K
ap2
-ap2
x2
(
ap2
as2 +as, il
.
GLOBAL STABILITY IN A CHAIN MODEL
239
This matrix has eigenvalues /_L,(1,1-.?2)-1>o
and
-ak as
f, i
2
+s
as,
1 and p2(1,0)> 1. The eigenvalues of IWE,) are ~~(1, 1) - 1 > 0 and c~~(l,O)- 1 > 0. E0 is a repeller. The eigenvalues of M(E,) are the same as that in Case 1. El is unstable in the x2 direction and stable in the xi direction. Similarly, M(E,) is similar to that in Case 2. E, is unstable in the xi direction and stable in the x2 direction. Finally, M(E,) has the same form as M(E,) and ME,). Thus, E, is locally asymptotically stable. In the following proposition, we describe the stable manifolds of all possible boundary equilibria when E, exists. The proof of this proposition is similar to the proof of Proposition 4.2 in Section 4 and is therefore omitted. PROPOSITION
5.2
For system (5.31, assume that (x:,x,*>
(1)
If E,,
exists.
and E, are the only boundary equilibria, then M+( E,) = {(O, x2) E R2: 0 < x2 < I}, Mf(E,)
(2)
If
={(x1,0)
ER2: 0 must tend to the unique interior equilibrium point. This implies that E, is globally W stable. As in Section 4, and in order to analyze the global stability of system (5.2), we lei k = (xi, x,;O,O,O) represent an equilibrium point of (5.2). All equilibria Ei, i = 1,2,3, of (5.2) exist under the same conditions of existence as (5.3). The possible equilibria of (5.2) in (x; U) coordinates are 8, = (0,0;0,0,0),
8, = (i,,0;0,0,0),
& = (0,
f,;O,O,O),
and &=(xF, The equilibria
x,*;O,O,O).
in (x; s) space can be described
as follows:
&=(0,0;1,1,0), ~‘,=(x,,O;l-~F,,l+a,~,,O), &=(O &=(x1*, The variational
,&;l,l-f,,ZG,), x~;1-x~,1+Lyix1*-x2*,x2*).
matrix of system (5.2) evaluated
at k is given by
where A is the 2x2 variational matrix M(E) of system matrix with entries api /au,, 0 is the 3 X 2 zero matrix, diagonal matrix with diagonal entries - 1. It is clear that this matrix are given by the eigenvalues of M(E) and the
(5.3), B is a 2x 3 and C is the 3 x 3 the eigenvalues of negative eigenval-
GLOBAL
STABILITY
IN A CHAIN
241
MODEL
ues of the matrix C. Thus, ki is stable if and only if Ei is stable. This implies that if _& exists, then kc,, E?,, and I$ are unstable and 8, is asymptotically stable. The next two theorems are the counterparts of Theorems 4.3 and 4.4 of system (4.2); their proofs are similar and therefore are omitted. The three cases to be considered in the lemma are @a, Et, k$, &, &, k,}, and @a, ii, & &1. THEOREM
5.4
For system (5.21, assume that k3 exists. Then system (5.2) persists; that is, euery solution (x( t ); u( t )) E I+‘; with xi(O) > 0 has
lim inf xi( t) > 0, t’m
i=1,2.
The following is the main result of this section, the global stability of kj. THEOREM
5.5
Let the hypothesis of Theorem 5.4 hold. Then k3 is a global attractor for system (5.2). Special Case. As in the special case in Section 4, set S, equal to zero. Let bZ(sZ) = pz(sZ,O). The analysis in the omega limit set is done by replacing pZ by fiii2(1+ (~ixt - x,>. The function GL2is a function of one variable, and afi, /ax, > 0 and $z /ax2 < 0. The conditions of existence of k.3 are the same-simply, change p&1,01 to b,(l) and pZ(l+ crt(pi(O),O) to /&(l + a,cp,(Oll. 6.
DISCUSSION
We have considered a simple chain of n species and it + 1 substrates in a chemostat culture. The species interact via intermediate metabolites in a manner characterized by assumptions Al-A3 of Section 3. We have assumed that the growth function k,, depends on S, and the product substrate S, + t. In Section 4, we assumed that only substrate St is supplied to the culture. Previous analysis of this case in [7] was restricted to a two-species, two-substrate chain with the assumption that the growth function p2 depends on S2 only. The result in this section is that the unique equilibrium point E, is globally stable. This was achieved by determining all possible equilibria and the conditions under which they exist and their stability. This led to an improvement in the conditions of existence of k’, given in [7] for the case n = 2. Next, we used the structure of the reduced system to show that solutions with initial value in the interior of rW: go to E,. Also, we showed that _!?‘, is globally asymptotically stable for any .
A. S. ELKHADER
242
solution of system (4.2) with the initial condition xi(O)> 0, 1 Q i < n. Finally, we considered the case when the substrate S,, , is set equal to zero. This does not affect the analysis, but the function pL, becomes a function of one variable. A biological interpretation of the results is that all species are able to survive and their mass density is asymptotic to an equilibrium if each species Xi is able to survive and produce enough nutrient for the species Xi+, for it to survive also. In Section 5, we studied a simple synthrophic chain of two species and three substrates. We assumed that substrates S, and S, are supplied to the culture. An analysis in [8] when II = 2 showed that an equilibrium with nonzero species exists and is locally asymptotically stable if the dilution rate D < min(pTaX, ~7~). The analysis in this section reveals that all equilibria with zero xi are unstable. Moreover, a unique equilibrium with positive xi exists in three different cases: (1) Species Xi is able to survive on its own, but X, is not. Coexistence is possible if X, produces enough nutrient for X, to survive. (2) Species X, is able to survive on its own, but X, is not. In this case coexistence is possible if X, produces enough nutrient for X, to survive. (3) Both species are able to survive on their own, and therefore they coexist. Finally, we have shown that all species persist and reach a steady state. The special case with S, set equal to zero was analyzed. This line analysis can be continued for any large n > 2. APPENDIX
A
Proof of Lemma
2.2
By Theorem 2.1, it is enough to show that the principal minors of M(n) alternate in sign. To show this we use induction. Let D(i) be the determinant of the matrix M(i), 1~ i < n.
D(1) =detM(l)
= -(l+
K,) o.
GLOBAL
STABILITY
IN A CHAIN
243
MODEL
Suppose D(n - 1) = (- l)“-‘(l+ K,_,(l+ pand D(n) with respect to the last column.
Kn-J..
. Cl+
K,)
. ))I.
Ex-
D(n)=-(l+K,)D(n-l)-K,_,D(n-2) = -(l+
K,)(-l)“-‘(l+
- K,_,(
-1)“-2(
=( -l)n{[l+
K,_,(l+ l+ Kn_2(
K,_,(l+ + K,[l+
Knp2(
K,_,(I+
-K,_,[l+K,_,(
Kn_2( K,)
. ..(l+
K,)
. ..(l+ Kn-*(
. ..(I+
..))J
+} ..))I K,)
.))]
..++K,)+]}
=(-1)“{1+K,(1+K,_,(1+K,_,(~~~(1+K,)~~~)))}. APPENDIX
K,)
. ..(l+
n
B
Here we state a theorem due to Markus [4] and give the proofs of Theorems 4.4 and 4.5. But before we give the statement of Markus’s theorem, we need the following definition. DEFINITION
Let A: ii = fi(x, t) and A,:
ii = fi(x), 1 < i < n, be a first-order system of
ordinary differential equations. The real-valued functions f,(x, t ) and fiCX>are continuous in (x, t) for x E G, where G is an open subset of R”, and for t > t, they satisfy a local Lipschitz condition in x. The system A is said to be asymptotic to A,( A --) A,) in G if, for each compact set K in G and for each E > 0, there is a T = T(K,e) > 0 such that 1f,(x, t>- fi(x>l < E for all i, all x E K, and all t > T. THEOREM [MARKUS]
Let A + A, in G, and let P be an asymptotically stable critical point of A,. Then there is a neighborhood N of P and time T such that the omega limit set of every solution x(t) of A that intersects N at a time later than T is equal to P.
In order to apply this theorem, we extend the domain of system (4.2) to an E-open set G of lR2”+’ containing Wnl, and the functions ~~ are defined in G for every i. Proof of Theorem 4.4. Let jz = (x(O); u(0)) E W,‘, Xi(O) > 0. Let (x(t); u(t)) be the solution of system (4.2) with (x(O); u(0)) = jz. We need to show that for 1 < i < n, liminf xi > 0 as t +m. To do so, we use induction. First, assume that lim inf x1 = 0. If lim sup x, = 0, then lim, +m x1 = 0. Since I+m f-m lim f ,m uj = 0 for all j, and ui - x, + xi-, > 0 in W,, it follows that for all i>l, and ~LI(u~-xl+1,u~-x2+x~~~~,~1,0~ as lim (,,xi=O
244
A. S. ELKHADER
t -+m. Since k&1,0) > 1, it is clear from system (4.2) that x1 + 0 is impossible unless x,(t) is identically zero. In particular, x,(O) = 0, a contradiction. Now suppose that limsupx,(t) > 0. Then we can choose a sequence {t,)
such that x.,(t,) = 0, ‘lim”, --lmt, =y and lim, -rm x,(t,) = 0. As above, this implies that lim,,, x;(t,,) = 0 for all i. By arguing as above, we obtain a contradiction. Thus, liminf x,(t) > 0. Now, suppose liminf x,(t) > 0 for n+m t-m l 0, it follows that lim,,,x,(t)=O for all r b m+l.Since liminfxi>Ofori~m,thereexistsapoint(y;”,...,y~,O,...,0; I’m 0 ,..., 0) E W,, yZY’> 0, 1 Q i Q m, in a(E). But the orbit through this point tends to i,,,, so l?,,, E a(%). Now we consider the following systems for l 0. We can choose a sequence {t,} such that t-m lim n-m t n =w, i.m+l(t,) = 0, and lim n_m x(t,> = 0. As above, we obtain for all r>m+l. Since lim,,,xi(t)>O for iE n(Z), x,” > 0. By arguing as above, we obtain a contradiction. Hence, liminf xm+l(t> > 0. t-m
W
Proof of Theorem 4.5. Theorem 4.4 implies that n(X) contains a point (a; u), Pi > 0, 1 < i < n. Since every point in R(K) must have ci = 0, 1 < i < n + 1, it follows that (f, 0) E R. The trajectory x(t) with initial value x(0) = f approaches E, by Proposition 4.2. Therefore, I?;, E Cl(%). But k, is asymptotically stable; then n(x) = i,,. Thus, I!?‘, is globally asymptotically stable. n
This work formed part of my doctoral thesis at Arizona State University. I wish to express my gratitude to Professor Hal L. Smith and the referees for their helpful suggestions. REFERENCES 1 2 3
H. I. Freedman
and P. Waltman,
Persistence
in models
of three
interacting
preda-
tor-prey populations, Math. Eiosci. 68:213-231 (1984). M. W. Hirsch, Systems of differential equations which are competitive or cooperative. I. Limit set, SIAMJ. Math. Anal. 13:167-189 (1982). M. W. Hirsch, Systems of differential equations which are competitive or cooperative. II. Convergence almost everywhere, SIAM J. Math. Anal. 16:432-439 (1984).
4
L. Markus, Asymptotically autoqomous differential systems, in Contributions to the Theory of Nonlinear OsciIation, Vol. 3, Princeton Univ. Press, Princeton, N.J., 1956, pp. 17-29.
5
S. J. Pirt, Principles of Microbe and Cell Cukation, Blackwell Scientific, Oxford, 1975. G. E. Powell, Equalization of specific growth rates for synthrophic association in batch culture,/. Chem. Technol. Biotechno[. 34:97-100 (1984).
6 7 8 9 10 11
G. E. Powell, Stable coexistence of synthrophic associations in continuous culture, J. Chem. Technol. Biotechnol. 35B: 46-50 (1985). G. E. Powell, Stable coexistence of synthrophic chains in continuous culture, J. Theor. Pop. Biol. 30:17-25 (1986). J. Smillie, Competitive and cooperative tridiagonal systems of differential equations, SIAM .I. Math. Anal. 15530-534 (1984). H. L. Smith, Systems of ordinary differential equations which generate an order preserving flow. A survey of results, SIAM Rec. 30:87-113 (1988). H. L. Smith, Competing subcommunities of mutualists and a generalized Kamke theorem, SZAM J. Appl. Math. 46:856-874 (1988).
ABSTRACT A mathematical model formulated by G. Powell in 1986 is considered. This model describes a synthrophic chain of species X,, 1~ i G n, in a continuous chemostat culture. Species X, utilizes substrate 5, and forms a product S,,,. Substrate S,+t inhibits the growth of X, and constitutes a growth-limiting substrate of X,,,. For this model, conditions for coexistence of all species are improved, and the local stability results extended to global stability. Some interesting cases are discussed.
1.
the are
INTRODUCTION
In this paper, we study the behavior of a mixed culture, that is, a mixture of different types of organisms. The interaction between species in a mixed culture can best be tested by means of a continuous culture using the chemostat technique, which has an extensive description in the literature (see, for example, Pirt [5]). Because of the complexity of mixed-culture behavior, mathematical models of various systems are particularly valuable for describing and predicting this behavior. A mathematical model describing the dynamics of a particular class of mixed culture was introduced by Powell [6-81. This class is a synthrophic association of two species or more. Following Powell [8], a synthrophic chain of n species Xi, 1~ i < n, and n + 1 substrates Si, 1~ i < n + 1, in continuous culture is of the form
In this synthrophic chain, the growth of species Xi depends on the substrates Si and Si+i. Although the growth of X, may depend on the presence of many substrates, it is the levels of Si and Sj+r that are assumed to be growth-limiting. Species Xi utilizes substrate Sj and forms a product This product inhibits the growth of species Xi but constitutes a si+l. growth-limiting substrate for species X,+i. In this manner, an incremental increase in species Xi leads to more favorable conditions for the growth of MATHEMATICAL
BIOSCIENCES
OElsevier Science Publishing 655 Avenue of the Americas,
104:203-245
(1991)
Co., Inc., 1991 New York, NY 10010
203 0025-5564/91/$03.50
204
A. S. ELKHADER
species Xi_, and Xi+,. A mixed culture of the two species Pseudomonas and Hyphomicrobium is a typical example of this kind of association. The Pseudomonas species utilizes methane and produces methanol, which inhibits its growth. The second species, Hyphomicrobium, utilizes the methanol produced but not methane, thus overcoming the inhibition of the first species and allowing both species to grow. This cooperative interaction between the species will be the basic idea of the analysis in this work. To be precise, a mathematical model of 2n ordinary differential equations, formulated by Powell [8] to describe the dynamics of n species and n substrates in a chemostat culture, can be reduced to a tridiagonal cooperative system of ordinary differential equations. Cooperative systems are defined in the sense of Hirsch [2, 31. This system will be analyzed, and the local stability results will be extended to global stability, which will be the focus of this work. This work is organized as follows: In Section 2, we introduce the notation and discuss some results essential to the understanding of later sections. In Section 3, we describe the model formulated by Powell [S] and modify the function pn. Some of the previous results are presented. Also, we give some basic properties of solutions of the system and show that the system can be reduced to an n-dimensional system of the tridiagonal cooperative type. This structure, a result due to Smillie [9], and a result due to Markus [4] are the main tools in determining the global stability of the systems of differential equations in Sections 4 and 5. Section 4 is devoted to analyzing the system when only substrate S, is supplied to the culture and the function CL,,depends on substrates S, and system starts by studying Sn+l. The analysis of the reduced n-dimensional the flow on the boundary. This leads to establishing necessary and sufficient conditions for the existence of all boundary equilibria and a unique interior equilibrium. Arguments based on the implicit function theorem in R2 are the main tools in determining these conditions. Simple arguments based on the properties of cooperative matrices are used to study local stability of all equilibria. Description of the stable manifolds of all boundary equilibria and the structure of the variational matrix at the interior equilibrium allow us to apply the result of Smillie [9], which is that the interior equilibrium point is globally asymptotically stable. To establish the global stability of the (2n + l&dimensional system, we describe all possible equilibria and analyze their local stability. Arguments based on the Markus theorem [4] are used to determine the global stability of the system. The special case when the function pL, depends on s, only is studied. The system in [81, when n is set to equal 2 and p2 depends on S, and S,, is considered in Section 5. The analysis of this five-dimensional system is organized as in Section 4. We show that the existence of an interior equilibrium point is possible in three different ways. This improves the conditions given in [8]. The local stability analysis and the global stability
GLOBAL STABILITY IN A CHAIN MODEL
205
analysis of all possible equilibria are given. This section is concluded considering the special case when the function pz depends on s2 only. The paper ends with a brief discussion. 2.
by
PRELIMINARIES
In this section we introduce some notations and nonnegative matrices and cooperative systems that later sections. Also, a result on the implicit function given. Following Hirsch [2, 31, we say that the system of i = F(x),
XER”,
F=
state some results on are essential tools in theorem in [w2 will be equations
(f,,...,fn>
(2.1)
is cooperative (competitive) if afi /Jxj > 0 (af, /axj < 0 for i # j, tridiagonal if afi /ax, = 0 for ]i - j] > 1, and strongly cooperative tridiagonal if af, /ax, is positive for Ii - j] = 1. Also, an it X n matrix A is called a stable matrix if S(A) < 0, where S(A) = maxRe A and A runs over the eigenvalues of A. We let x,):
rwT={(x,,...,
xi > 0, 1 < i < n) ;
k!: is the interior of rWF containing the The next theorem, which will be used condition for determining the stability off-diagonal elements. It is part of a matrices in Smith [lo, 111.
positive vectors. repeatedly, provides a very simple of a matrix A with nonnegative general theorem on nonnegative
THEOREM2.1 Zf A is an n x n matrix with nonnegatiue off-diagonal elements, then S(A) < 0 is equivalent to the condition that the principal minors of A alternate in sign. That is,
A direct consequence of this theorem is the following proof is given in Appendix A. Let M(n) be -l(l+K,)
whose
K,
-(l+K)
. . ..f.............2....i.......................
K,
) -Cl+
Ki > 0, 1 Q i < n.
lemma,
K-d 1
Kn-, -(l+K,)
206
A. S. ELKHADER
LEMMA 2.2
The matrix M(n)
is a stable matrix.
(See Appendix
A.) The next theorem, due to Smillie [9], is one of the important tools in determining the asymptotic behavior of the solution of the tridiagonal cooperative systems, which will appear in later sections. THEOREM 2.3 [Smilliel
Let j, = F(x) be a strongly cooperative tridiagonal system of ordinary differential equations defined on R c R”. Assume that the component functions F, are n - 1 times differentiable. Let cp(t> be a solution of d(t) = F(cp), defined on a maximal interval of the form [0, a), 0 < a 0
t-m
for all i
whenever xi(O) > 0 for all i. DEFINITION
System (2.1) is uniformly persistent if there exists 6 > 0 such that whenever xi(O) > 0 for all i, then
lim infxi(t)>6>0 t+m
for all i
Finally, the following lemma, which is a consequence of the implicit function theorem in R2, will be used repeatedly. The lemma is illustrated in Figure 1. LEMMA 2.4
Let cp and IJ be continuously differentiable functions
on 0 f x < a. As-
sume that &a> = a, x < cp(x) for 0 < x < a, dp(x)/dx < 1, 0 < $(O), and $(x) is a linear function with positive slope. Let p = p(r,, r2) be a continu-
GLOBAL
STABILITY
IN A CHAIN
207
MODEL
Y
(ON
(W
(l?O)
FIG. 1. The function
x = q(y)
as described
in Lemma
3.4.
ously differentiable function on iwl, and let
For (x, y) E D, let r, = q(x)-
F(x,Y)
x and r2 = $,(x)-
=p(r1(x),r2(x,y)) -1
y. Define onD.
_ Assume that $L /Jr, > 0 and ap /ar, < 0. If F(x, y) = 0 at (.?,O) and (x, y), where 0 < .C < X < a and j = 1,9(x>, then there exists a unique smooth function x = v(y), 0 < y < J, such that F(T$ y), y) = 0, $0) = 2, and $ji) = X. Moreouer, dq/dy < l/(d$/dx).
208 Proof.
regions
A. S. ELKHADER
Let
The function in x because
Similarly,
@L(x)= I/J(O)+ (YX, LY> 0. Divide
the region
D
into
two
D, and D,, where
F(x, y) is continuously
F(x, y) is increasing
differentiable
on D and is decreasing
in y because
C?F
-=$+)>O. 8Y
Now consider the lines in D, of the form y = 0, 0 < 8 < $(O>. The function F(x, y) evaluated along each line y = 0 is given by
F(x,B)=~(cp(x)-x,lCI(~)-e)-l =/.&(cp(x)-x,$(O)+ax-8)-l. Since p’(x) < 1, F(x,B)
is strictly decreasing
in x for each 0. Moreover,
F(~,B)=II(~(~)--,IC~(~)--)-~
>p(cp(l)-i,cL(-i))-l = F(i,O) = 0, and F(x,B)=p(cp(Z)-x,$(-3:)-0)-1 = 0. Observe that ~(0) = 2, and let 2, = 7($(O)). The implicit function theorem guarantees that 77 is continuously differentiable. Next, consider the lines y = $,(O)+ ix, 0 ( 6 < (Y. These lines lie in D,. The function F(x, y) evaluated along any of those lines, for a fixed e, is
209
GLOBAL STABILITY IN A CHAIN MODEL
given by F( x,$(O) +
ix) = /A(q(x) - x,$(O) + ax -$(O) - w-1 =&+7(x)-X&X-$)x)-l.
It is easily seen that the function
F is decreasing
and satisfies
F(~~,~(O)+B~,)=CL('P(~~)-~~,(~--)~~)-~ > /_L( cp(&J -
io, “1,) -
1 = 0,
F(i,$(O)+i&)=&(X)-.?,(a-8)3)-l in D implies that x = q(y). Moreover,
dv
dy=
-$-1)
- c3F/ay t’F/ax
where k =
(aF/c+,(dd~ aF/ar,
-1))
> 0.
n
Remarks. (1) If, in Lemma 2.4, the condition F(_f,O) = 0 is replaced by F(0, F), 0 < jj 0.
(3.1)
Here, x0 = 0, xi denotes the cell density of species X,, si denotes the concentration of substrate Si, so is the concentration of substrate Si entering the chemostat, x (g/mole) denotes the growth yield of species X, with respect to its substrate Si, y.,i_ 1 (mole/g) denotes the product yield of substrate Si by the species Xi-i, and pi denotes the specific growth rate of species X, and depends on the concentrations si and si+, of substrates Si and Si+i, respectively. D is the dilution rate. As usual, the dot denotes the time derivative. The only difference between system (3.1) and the model in [8] is in the growth function p,. Here, CL,,depends on Sri+++in the same fashion pi depends on Si+i. Thus, we include the dynamics on the substrate SE+,, which is assumed to be zero in [8]. In order to characterize the interaction sufficiently without postulating particular functional forms for the specific growth rate functions, a number of qualitative assumptions have to be made. For convenience these are considered collectively: (Al) The specific growth rate of species concentrations si and si+i, pLi= I_L~(s~,s~+~), aPi(si,si+l) asi
>
aWi(si~si+l)
0
3
$+1
Xi depends
B pyax. Typically, the maximum specific growth rate and saturating /4= is attained in the limit of zero product concentration substrate concentration. (A31 In the absence of substrates Si and Si+,, the growth function pi is zero, that is, ~~(0,0)= 0. Substrate S,+i is derived solely from the growth of species X,. Some of the specific forms of the function pi frequently used (see Powell [8]; Pirt [S]) are
s -(W,/W,)p = Pmax s + K, +(k,/k,)p
P(S,P)
’
where k,, k,, k,, k,, k,, and k, are constants. The special case of system (3.1) when II = 2, with the assumptions that S, is solely derived from the growth of species X, (i.e., si = 0) and that the specific growth function pZ depends on S, only [i.e., S, is removed from the chemostat or ignored; pz = p&l] was analyzed by Powell [8]. It was shown that when the dilution rate D is less than the critical dilution rates for each species X, and X, (i.e., D < min{~;l”,~~ax}), a locally asymptotically stable steady state exists. System (3.1) with the assumption that Sn+l = 0 was addressed by Powell [8]. It has been shown that when the dilution rate D is less than the critical dilution rates of each species Xi, i=l , . . . , n (D < (ELF;“,pLr;l=,. . . , pzau}), a locally asymptotically stable equilibrium point with all species nonzero exists. Some of the previous results will be improved, and the global stability question in both cases will be answered. The special case when p,, depends on S, only will be discussed. 3.2.
SCALING
Before we start to analyze the above cases, it will be more convenient to reduce the number of independent parameters by scaling the system (3.1). This can be done if we perform the following substitutions: Ei = xi/kj,
si = Si/li, ai-
= y,i_,Ki_,
i;ii =
/liT
Pi(Fi,si+l)
D
’
ao=o, (3.2)
A. S. ELKHADER
212
where I, =$>O, ki=xli, for each i, 1 o mii+, = xil 2 ( asi+l i +i
m,l_l = xjz
Matrix M,, is of size i x(n - i) with (LV,,)~, = elsewhere. Matrix M,, is of size (n - i) x (n - i) except (A&J, =~~+~(xj,O)-l> 0. Matrix 0 is eigenvalues of M(E,) are the eigenvalues of tridiagonal cooperative matrix M,,, let kj = and
-Gj/asj+l
ai+ias, >o
xjapi /as2 > 0 and zero with - 1 in the diagonal of size (n-i>Xi. The M,, and M,,. For the
222
A. S. ELKHADER
For 1 d j d i, let Cj be the jth principal
-(1+kJ 1 cj = lxj
minor of M,,,
k, -(l+G . .
k, . .
. . 1
-(l+
kj)
The determinant on the right-hand side is the jth principal minor of the stable matrix described by Lemma 2.2. Since aj > 0, the signs of the Cj’s are the same as the signs of the principal minors of the stable matrix of Lemma 2.2. Hence, M,, is a stable matrix. The eigenvalues of M,, are (M,,),, > 0 and the n -(i + 1) diagonal entries mkk = - 1, i + 2 < k < n. Thus, Ei is a saddle point, stable in the xj direction j f i + 1, and unstable in the xifl direction. Finally, the variational matrix evaluated at E,, M(E,), is a tridiagonal cooperative matrix with aPi
mii_, = x:-g
> 0, 1
mi,i = x:
-
-aPi dSi
i
miifl = xi nG i
-
+
‘Pi dsi+l
O.
An analysis similar to that followed for the matrix M,, shows that ME,) a stable matrix and therefore E, is locally asymptotically stable.
is
Global Stability Analysis. To determine the global stability of the reduced system (4.3), we describe the stable manifolds of all possible equilibria. Let M+(Ei) be the stable manifold of Ei. PROPOSITION 4.2
For system (4.3), assume that E, exists. Then for each i, 0 < i < n,
Proof: We prove the proposition by induction. Let Bi = {x E H,‘: xi > 0 x,(t)) be the solution of system if and only if j 0, then x,(t) + 0 for all i. But since ~~(1~0) > 1, it is clear from system (4.3) that x,(t) does not tend to 0 unless x,(t) is identically 0. In . . particular, x y = 0, a contradtction.
GLOBAL STABILITY IN A CHAIN MODEL
Next, suppose that M+(E,)
223
= B,, 1~ i. We want to show that
M+(Ei+J = 4+1. To do so, let x0 E Bi+i, and hence xz = 0 for k > i + 2. The solution x(t) of this system tends to an equilibrium by Smillie’s theorem, Theorem 2.3, as t tends to infinity. This equilibrium must be one of E,, E,, . . . , Ei+ ,. But if this equilibrium is not Ei+,, then we have a contradiction to the hypothesis, since xi”+, > 0 and x0 is not on the M+(E,) for some I ,xi”+i > 0; then the trajectory through this point must go to the equilibrium (cpi 0 . . . 0 cpi(.Ci+i), . . . , fi+i, 0 , . . .,O)as t goesto infinity. Since x,+,(t) -+ Zi+, and x,(t) + 0, r > i +2, it follows that ~,+~(Xi+~-Xi+2,Xi+2-Xi+3)~~i+2(~i+~,o)as t+w. Since /Li+2Cfi+l,0) > 1, it is clear from system (4.3) that xi+2 + 0 is impossible unless xi+2 is identically zero. In particular, xi”,, = 0, and thus xi = 0 for k > i +2 since x0 E H,‘. Hence, M+(Ei+,)= Bi+l. W A corollary
to this proposition
is the following.
COROLLARY4.3
For system (4.3), assume that E, exists. Then E,, attracts all solutions with positive initial conditions in H,‘. The final result in this section is the global stability of system (4.2). To show that we let si, 0 < i 6 n, denote the equilibria of system (4.2). The equilibria in (x; v> coordinates are given by
8,=(0
0;o
)...)
8,=(&O
)...)
0)
)...)
0;o
,...)
~2=(qo1(~2),i*,o )...)
The variational
0) 0;o
,...)
0)
matrix of system (4.2) is of the form
where A is equal to M, the 12x n variational matrix of system (4.31, B is an n x(n + 1) matrix with bii = +L~/~v~, bi,i+i =$.L~/~v,+,, and bij = 0 for
224
A.S.ELKHADER
C is an (n + 1) X (n + 1) diagonal matrix with diagonal entries - 1, and 0 is the (n + 1) X n zero matrix. The eigenvalues of the variational matrix of system (4.2) at f?; are equal to the eigenvalues of M(E,) and the negative eigenvalues of the diagonal matrix C. This implies that ii is stable if and only if Ei is stable. The persistence of solutions and the global stability of I?,, are given by the next two theorems. The main tool in the proof of these theorems is a theorem due to Markus [4]. The statement of the theorem of Markus and the proof of the next two theorems are given in Appendix B.
j # i, i + 1. Matrix
THEOREM
4.4
For system (4.2), assume that k,, exists. Then system (4.2) persists; that is, for all i, 1 0 as t -+m for every solution with initial condition (x(O); ~(0)) E Wn’, xi(O) > 0 for all i.
Finally, the global stability of 8, is given by the following theorem. THEOREM
4.5
Under the hypothesis of Theorem Special Case.
or that
its value
4.4, .&, is globally asymptotically stable.
Assume that substrate S,,+, is removed from the culture is very small so it can be neglected. Mathematically
FIG.4. The effect of substrate
S,,
, on X,.
GLOBAL
STABILITY
IN A CHAIN
MODEL
225
speaking, the specific growth function pn depends on s, only; that is, pn is a function of one variable. In this case, in the omega limit set, s, = x,-i x,. Let ji, = p,,(x,_ 1- x,). It is clear that the only term affected by the removal of s,+l is the function pn. This suggests that an interior equilibrium point exists if and only if p,(l,O) > 1,
All the results in this section hold in this case. Furthermore, since species X,, has no inhibitor, its density at the equilibrium point Z,* is larger than .?‘,* when S,, , is present. A geometric illustration is given in Figure 4, where the slope of the function CL,,is less than the slope of the function fin, and therefore the density of the X, species at equilibrium is greater than that when Sri+++is present. 5.
5.1.
ANALYSIS OF A CHAIN OF TWO SPECIES AND THREE SUBSTRATES WHEN S, AND S, ARE SUPPLIED TO THE CULTURE DESCRIPTION
In this section, our main goal is to study system (3.2) when II is set equal to 2. In this case, the chemostat culture consists of two species, X, and X,, and three substrates, S,, S,, and S,. The scaled system of ordinary differential equations to be analyzed is given by
Sj(0)
~ O
for j= 1,2,3.
a, is a positive constant. As in Section 3, let S1(XI,UI) = 1+ ut - xt, s2(x,,x2,u2) s3(x2,uJ
= 1+ c’2+ (Ytx, - x2, =
u,+
x2.
(5.1)
226 Expressing
A. S.
ELKHADER
pI and pZ in terms of x,, x2, ur, u2, and u3 gives
Thus, system (5.1) is equivalent
to
~.1=x,[~~(x,,x2,uI,u2)--l],
~2=x2[~2(x*,x2,U~,U*,U3)-l],
l,,
GLOBAL STABILITY IN A CHAIN MODEL
where it is the solution
of pt(l-
227
x1,1+ atxt)
= 1;
E,=(%f,), where _E, is the solution
of pz(l-
x2,x,> = 1; and
x1, l+ ~ytxt - x,) = 1 and pL2(1+ where is the solution of pt(lx2, x2) = 1. We proceed to find conditions under which such equilibalxlria exist. First, let
The simplest equilibrium point is E, = (O,O>. This equilibrium Next we consider one-species equilibria.
always exists.
(a) Suppose x2 = 0. The function F,(x,,O) = ~r(l- x1, l+ a,~,)-- 1 is strictly decreasing. Its maximum value is attained at x1 = 0, and its minimum at x1 = 1. Therefore, F,(l,O) f F,(x,,O) d F&0,0). Since F,(l,O) = prCO,l+ a,>- 1 < 0 because p,(O,O) < 0, then F,(x,,O) = 0 has a nonnegative solution if and only if F,(O,O) > 0. If F,(O,O) = ~t(l,l)1 = 0, then x1 = 0 is the only solution of F&x1,0) = 0. If F,(O,O) = p(l,l> - 1 > 0, then there exists a unique it, 0 < .?t < 1, such that F,(i:,,O)= 0. Let Et = (,Ct,O). Then E, exists if and only if p,(l,l) > 1. (b) Suppose x1 = 0. The function F&O, x,) = p&lx2, x,)- 1 is strictly decreasing, so it attains its maximum value at x2 = 0 and its minimum at x2 = 1. This implies that F,(O, 1) < FJO, x2) < &CO, 0). But F,(O, 1) = ~~(0,1) - 1 < 0 because pJO,O) < 0. Then FJO, x,) = 0 has a nonnegative solution if and only if F-JO, 0) > 0. If F,(O, 0) = 0, that is, p#,O) - 1 = 0, then F,(O, x,) = 0 has only the trivial solution. If F,(O,O) > 0, that is, p2(1,0) > 1, then &(0,x,) = 0 has a unique solution .C2, 0 < 1, < 1, such that F,(O,X,) = 0. Let E, =(O,f,). Hence, E, exists if and only if &(l,O) > I. Finally, we consider two-species algebraic system of equations
F,(x,,x,) = 0, has a nonnegative solution Al-A3, on the functions it
equilibria.
This case occurs when the
F,(Xl,X,) = 0
in H;. From the qualitative and CL*,the partial derivatives
(5.4) assumptions, aF, /ax1 < 0,
A. S. ELKHADER
228
aF,/ax,> 0,dF,/dx,> 0,and aF,/dx,and FJx,,x2)attain their maximum and minimum values at the boundary of H;. We start with F,(x,,xz)=O. Since dF,/ax, 0,then minF,(x,,x2)=F,(1,0)=~I(0,1+~1)-1 1, then there exists a nontrivial solution of F,= 0 in H:. Furthermore, this solution lies on a curve of the form x1 = cpl(x,) such that F,(cp,(x& x,)= 0. This will be realized by finding solutions of F,= 0 at the boundary of Hz and then applying Lemma 2.4 as follows. We evaluate Fl(x,, x,)along the boundary lines. Along x1= 1, F,(l, Along
x2) = /_~(0,1+ qx2)
- 1 < 0.
x2 = 0,
As above, F,(x,,O)= 0 at (Z:,,O) if and only if ~Jl,l)> Along x1 = 0,
1.
This is an increasing function; its maximum its F,(O, 1) = pJl,O)1 > 0. Thus F,(O, x,)= 0 has a nonnegative solution if and only if min F,(O, x2)= F,(O,O) = pl(l,l)1 < 0. If p,(l,l) = 1, then (0,O) is the only solution of F&O,x2)= 0. If ~~(1~1) < 1, then there exists XZ, 0 < X2 < 1, such that
F,(O,i,)= 0. Fmally, along x2 = 1+ (Y~x,, F,(x,,l+qx,)=/.q(l-x1,0)-1. This is a decreasing
function,
with
minF,=F,(l,l+cu,)=~(O,O)-l 1, then we need to consider the following cases: (1) 1 , ,F, < xl. Thus, by the remark following Lemma 2.4, there exists a smooth function xt = ‘pi(x,), 0 = (pi(X& that satisfies xi = P&X:), and FI(9i(x2), x2) = 0. F2(x1, x2) = 0. Since dF2 /ax, P- 0 and dF2 /ax, < 0, then
Next, consider
minF,(x,,x,) H:
= F2(0,1)
=/+(0,1)-l
This implies that FJx,, x2> = 0 has a nonnegative
~0.
solution
if
maxF,(x,,x2)=F,(1,0)=/L,(1+a,,0)-120. H: If F,(l,O) = 0, then (1,O) is the only solution of F2(x,, x2> = 0. If F,(l,O) > 0, that is, p2(1+ (Y,,0) > 1, then there exists a nonnegative solution of F2 = 0 in Hz. Furthermore, F2 = 0 along a curve of the form x2 = ~&xi). As before, this will be shown by Lemma 2.4 as follows. According to Remark 2 following Lemma 2.4, we need to show that F, = 0 along the boundary of H;. Along xi = 1,
&(I, 4 This function
= ,q(l+
is a decreasing
and a maximum
al-
~2,
function
~2) - 1,
0
0.
Therefore,
a unique
solution
x,, 0 < xz < 1 + (it, exists such that F~(x,,Q)
at (1, X,X
= 0
GLOBAL
STABILITY
IN A CHAIN
231
MODEL
Along x2 = 0,
This is an increasing
function
with
Thus Fz(xl, 0) = 0 has a solution
if and only if
minF,(x,,O)=F,(O,O)=pL(l,O)-ll. If pJl,O) > 1, then there exists an Z2, 0 < T2 < 1, such that F2 = (XI, x2>= 0 at (0, JZ~). Finally, along x2 = 1+ (~rx,, F,( x,, 1+ cqxl)
= p2(0,1 + (~rx,) - 1 < 0.
This analysis shows that if ~~(1 + (~r,0) > 1, then F, = 0 along three different boundary lines, and therefore we consider the following cases. (1) 1 < p&1,0) < pJl+ a,,O) (see Figure 6a). In this case F2(x1,x2) = 0 at (1,x,) and (X,,O). By the remark following Lemma 2.4, there exists a unique smooth function x2 = (p&x,) that satisfies
0 =4X,), and
X2=p2(1),
+CY,,
A. S. ELKHADER
232 x2
x2
4
I
5’ n (17F22)
(1,522)
I
Xl
a. 1 < pz(l,O)
x1
< p2(1 + W,O)
b. p&,0)
c. clz(l, 0) < 1 < P2(1 t FIG. 6. The isocline
-
Xl
(,g
= 1 < p2(1 t QI,O)
0170)
F2(xlr x,) = 0 as in the discussion
of solutions
of system (5.4).
4 (2) l+p2(1,0) = 0. Substituting xi = (P(Q) in F*, we obtain F,(cpd4>x,)
Since dye, /dyz < l/cri,
Therefore,
=/12(1+~1(~2)-
then F2(~i(x2),
x2>+>
x2) is a decreasing
F2((p,(x2), x,) = 0 has a solution maxk(cp,(x,),x,)
-1.
function
with
if and only if
=&(‘~i(%),%)
z=O.
If ~,(l, 1) < 1, then X2 > 0 and cpi(X,) = 0. Thus, maxF2(cp2(x2),x2)
=~~(1+~,4+(%)-
%,&)
-1
=y~(1-~-,,x,)-1 $2
x2
Cc)
(44)
(4
~~~~ 7. The isocline F, = 0 and F2 = 0 as described in Case 1 in Section 5. (a) ~JI,I)< 1 1, then there exists a unique function x2 = (p2(x1) such that
Any solution of system (5.4) must lie on this curve. Therefore, x2 = (p2(x1), 0 G x, Q 1, into F,(x,, x,1, we obtain
4(wP*(XJ) =/-+(l-
substituting
Xl,l+Qlxl-(P2(X1))-1
Again, this is a strictly decreasing function. Since its minimum is F,(l, (~~(1)) = &O, 1 + (Y,- x,1-- 1 < 0, then F,(x,, cp,(x,)) = 0 has a nonnegative solution if and only if max F,(x,, (p2(x,)) > 0, that is,
If /*i(l, 1 - (~~(01)= 1, then xi = 0 is the only solution of F,(x,, (pJx,l) = 0. If ~,(l,lcpJ0)) > 1, then there exists a unique XI*, 0 < x;” < 1, such that F,(x,*, cp,(x$>> = 0. Hence E, = 1
and
cL~(l,l~~l~cL1(l~l-(P2~O~~~
where .Lz = ‘p2(0) and XT = (pJxf>. Case 4. ~i(l,l)> 1 and &,O) > 1 (see Figure 819. In this case, j.~i(l,l) > 1 implies cpi(O)= 1, > 0, and since ~J1,0) > 1, then /+(l+ (~i(pi(O),O) > 1. Thus the conditions for the existence of an interior equilibrium E, in Case 2 are satisfied. Similarly, the conditions for the existence of E, in Case 3 are satisfied because ~~(1,0> > 1 implies (~~(0) = Z2 > 0 and p,(l,l) > 1 implies pi(l,l(pJO>)> 1. Hence, E, exists uniquely if and only if ~i(l,l) > 1 and ~J1,01> 1. Now let 2 be the set of all possible equilibria of system (5.3). Then the above analysis concludes the proof of the following theorem.
GLOBAL
STABILITY
THEOREM
IN A CHAIN
237
MODEL
5. I
For system (5.3), the following hold: (a> Z = (Es} if and only if ~~(1, 1) < 1, /.&,O) < 1. (b) C = {E,, E,} if and only if pr(l, 1) > 1, ~,(1,0> < 1. (c) Z = {E,, I&) if and only if p,(l,l) d 1, pJ1,O) > 1. (d) C =(E,,E,,E,) if and onfy if p#,l)> 1, P&~,O)G 1 1. Here, x2 = (p2(x1) satisfies Fz(xl, x2) = 0 and ‘pz(O)= i2. (f) C = (E,, E,, E,, EJ if and only if pJl,l) > 1 und PL~(~,O)> 1. Local Stability Analysis. Since we are interested in coexistence between species, we study only the cases in which an interior equilibrium exists. We have shown that coexistence between species X, and X, is possible in three different ways. We will study the local stability of all possible equilibria of system (5.3) in each case. Let (XT, x$> be the interior equilibrium point. The Jacobian (variational) matrix of system (5.3) evaluated at an equilibrium point E is of the form
where m 11
=
(
--al+ -+ar,s
aI+
dS1
x,+~.,(1-~x,,l+CY,x,-xX2)-l,
2 1
-a/-h
m12 =TX1’ m21 = m 22
=
+2 “la~,X27
(
-ap2
-
as,
+%as,
1
x
2
+/.l(l+a 2
The three cases of C to be considered Case 1.
1x 1
> 1
and
P2(LO)
where .Cr = p,,(O) > 0. The variational
2,x2>
-1.
are the following.
C = {E,, E,, E3}. This case happens
/4(l,L)
-x
0 and ~J1,0)1 < 0. Ea is a saddle point, unstable in the x1 direction and stable in the x2 direction. The equilibrium E, has M(E,) as a variational matrix,
The eigenvalues
of M(E,)
are /_LCLZ(lf(Yr~,,O)-l >o.
and Thus, Et is unstable in the x2 direction Similarly, at E,,
and stable in the x1 direction.
@l
-“?-
+“‘as, ME,) =
** alx2 as,
\ The off-diagonal are
\
aP1
elements
of this matrix are positive.
Its principal
minors
and
ah ap2 ah ap2 * *!?L%_ x*x*--+cYx*x*-->O detWE3)=~lx2 as,as2 1 2 as,as, 1 1 2 as,as, because aP2 /as, < 0. Then Theorem asymptotically stable steady state.
2.1 implies
that
Es is a locally
Case 2. 2 ={E,, E,,E3}.The conditions under which this can happen are pr(l,l) < 1 < p&1,1 - ‘p2(0)) and &,O) > 1, where P,(O) = f2 > 0. The variational matrix evaluated at E, takes the same form as in Case 1 with eigenvalues p&1,1)- 1 < 0 and p2(1,0) 2 1 > 0. Hence E, is a saddle point, stable in the x, direction and unstable in the x2 direction. At E,, M(E,) takes the form 0
/.&(1,1-&-l
WE21
=
1
G2K
ap2
-ap2
x2
(
ap2
as2 +as, il
.
GLOBAL STABILITY IN A CHAIN MODEL
239
This matrix has eigenvalues /_L,(1,1-.?2)-1>o
and
-ak as
f, i
2
+s
as,
1 and p2(1,0)> 1. The eigenvalues of IWE,) are ~~(1, 1) - 1 > 0 and c~~(l,O)- 1 > 0. E0 is a repeller. The eigenvalues of M(E,) are the same as that in Case 1. El is unstable in the x2 direction and stable in the xi direction. Similarly, M(E,) is similar to that in Case 2. E, is unstable in the xi direction and stable in the x2 direction. Finally, M(E,) has the same form as M(E,) and ME,). Thus, E, is locally asymptotically stable. In the following proposition, we describe the stable manifolds of all possible boundary equilibria when E, exists. The proof of this proposition is similar to the proof of Proposition 4.2 in Section 4 and is therefore omitted. PROPOSITION
5.2
For system (5.31, assume that (x:,x,*>
(1)
If E,,
exists.
and E, are the only boundary equilibria, then M+( E,) = {(O, x2) E R2: 0 < x2 < I}, Mf(E,)
(2)
If
={(x1,0)
ER2: 0 must tend to the unique interior equilibrium point. This implies that E, is globally W stable. As in Section 4, and in order to analyze the global stability of system (5.2), we lei k = (xi, x,;O,O,O) represent an equilibrium point of (5.2). All equilibria Ei, i = 1,2,3, of (5.2) exist under the same conditions of existence as (5.3). The possible equilibria of (5.2) in (x; U) coordinates are 8, = (0,0;0,0,0),
8, = (i,,0;0,0,0),
& = (0,
f,;O,O,O),
and &=(xF, The equilibria
x,*;O,O,O).
in (x; s) space can be described
as follows:
&=(0,0;1,1,0), ~‘,=(x,,O;l-~F,,l+a,~,,O), &=(O &=(x1*, The variational
,&;l,l-f,,ZG,), x~;1-x~,1+Lyix1*-x2*,x2*).
matrix of system (5.2) evaluated
at k is given by
where A is the 2x2 variational matrix M(E) of system matrix with entries api /au,, 0 is the 3 X 2 zero matrix, diagonal matrix with diagonal entries - 1. It is clear that this matrix are given by the eigenvalues of M(E) and the
(5.3), B is a 2x 3 and C is the 3 x 3 the eigenvalues of negative eigenval-
GLOBAL
STABILITY
IN A CHAIN
241
MODEL
ues of the matrix C. Thus, ki is stable if and only if Ei is stable. This implies that if _& exists, then kc,, E?,, and I$ are unstable and 8, is asymptotically stable. The next two theorems are the counterparts of Theorems 4.3 and 4.4 of system (4.2); their proofs are similar and therefore are omitted. The three cases to be considered in the lemma are @a, Et, k$, &, &, k,}, and @a, ii, & &1. THEOREM
5.4
For system (5.21, assume that k3 exists. Then system (5.2) persists; that is, euery solution (x( t ); u( t )) E I+‘; with xi(O) > 0 has
lim inf xi( t) > 0, t’m
i=1,2.
The following is the main result of this section, the global stability of kj. THEOREM
5.5
Let the hypothesis of Theorem 5.4 hold. Then k3 is a global attractor for system (5.2). Special Case. As in the special case in Section 4, set S, equal to zero. Let bZ(sZ) = pz(sZ,O). The analysis in the omega limit set is done by replacing pZ by fiii2(1+ (~ixt - x,>. The function GL2is a function of one variable, and afi, /ax, > 0 and $z /ax2 < 0. The conditions of existence of k.3 are the same-simply, change p&1,01 to b,(l) and pZ(l+ crt(pi(O),O) to /&(l + a,cp,(Oll. 6.
DISCUSSION
We have considered a simple chain of n species and it + 1 substrates in a chemostat culture. The species interact via intermediate metabolites in a manner characterized by assumptions Al-A3 of Section 3. We have assumed that the growth function k,, depends on S, and the product substrate S, + t. In Section 4, we assumed that only substrate St is supplied to the culture. Previous analysis of this case in [7] was restricted to a two-species, two-substrate chain with the assumption that the growth function p2 depends on S2 only. The result in this section is that the unique equilibrium point E, is globally stable. This was achieved by determining all possible equilibria and the conditions under which they exist and their stability. This led to an improvement in the conditions of existence of k’, given in [7] for the case n = 2. Next, we used the structure of the reduced system to show that solutions with initial value in the interior of rW: go to E,. Also, we showed that _!?‘, is globally asymptotically stable for any .
A. S. ELKHADER
242
solution of system (4.2) with the initial condition xi(O)> 0, 1 Q i < n. Finally, we considered the case when the substrate S,, , is set equal to zero. This does not affect the analysis, but the function pL, becomes a function of one variable. A biological interpretation of the results is that all species are able to survive and their mass density is asymptotic to an equilibrium if each species Xi is able to survive and produce enough nutrient for the species Xi+, for it to survive also. In Section 5, we studied a simple synthrophic chain of two species and three substrates. We assumed that substrates S, and S, are supplied to the culture. An analysis in [8] when II = 2 showed that an equilibrium with nonzero species exists and is locally asymptotically stable if the dilution rate D < min(pTaX, ~7~). The analysis in this section reveals that all equilibria with zero xi are unstable. Moreover, a unique equilibrium with positive xi exists in three different cases: (1) Species Xi is able to survive on its own, but X, is not. Coexistence is possible if X, produces enough nutrient for X, to survive. (2) Species X, is able to survive on its own, but X, is not. In this case coexistence is possible if X, produces enough nutrient for X, to survive. (3) Both species are able to survive on their own, and therefore they coexist. Finally, we have shown that all species persist and reach a steady state. The special case with S, set equal to zero was analyzed. This line analysis can be continued for any large n > 2. APPENDIX
A
Proof of Lemma
2.2
By Theorem 2.1, it is enough to show that the principal minors of M(n) alternate in sign. To show this we use induction. Let D(i) be the determinant of the matrix M(i), 1~ i < n.
D(1) =detM(l)
= -(l+
K,) o.
GLOBAL
STABILITY
IN A CHAIN
243
MODEL
Suppose D(n - 1) = (- l)“-‘(l+ K,_,(l+ pand D(n) with respect to the last column.
Kn-J..
. Cl+
K,)
. ))I.
Ex-
D(n)=-(l+K,)D(n-l)-K,_,D(n-2) = -(l+
K,)(-l)“-‘(l+
- K,_,(
-1)“-2(
=( -l)n{[l+
K,_,(l+ l+ Kn_2(
K,_,(l+ + K,[l+
Knp2(
K,_,(I+
-K,_,[l+K,_,(
Kn_2( K,)
. ..(l+
K,)
. ..(l+ Kn-*(
. ..(I+
..))J
+} ..))I K,)
.))]
..++K,)+]}
=(-1)“{1+K,(1+K,_,(1+K,_,(~~~(1+K,)~~~)))}. APPENDIX
K,)
. ..(l+
n
B
Here we state a theorem due to Markus [4] and give the proofs of Theorems 4.4 and 4.5. But before we give the statement of Markus’s theorem, we need the following definition. DEFINITION
Let A: ii = fi(x, t) and A,:
ii = fi(x), 1 < i < n, be a first-order system of
ordinary differential equations. The real-valued functions f,(x, t ) and fiCX>are continuous in (x, t) for x E G, where G is an open subset of R”, and for t > t, they satisfy a local Lipschitz condition in x. The system A is said to be asymptotic to A,( A --) A,) in G if, for each compact set K in G and for each E > 0, there is a T = T(K,e) > 0 such that 1f,(x, t>- fi(x>l < E for all i, all x E K, and all t > T. THEOREM [MARKUS]
Let A + A, in G, and let P be an asymptotically stable critical point of A,. Then there is a neighborhood N of P and time T such that the omega limit set of every solution x(t) of A that intersects N at a time later than T is equal to P.
In order to apply this theorem, we extend the domain of system (4.2) to an E-open set G of lR2”+’ containing Wnl, and the functions ~~ are defined in G for every i. Proof of Theorem 4.4. Let jz = (x(O); u(0)) E W,‘, Xi(O) > 0. Let (x(t); u(t)) be the solution of system (4.2) with (x(O); u(0)) = jz. We need to show that for 1 < i < n, liminf xi > 0 as t +m. To do so, we use induction. First, assume that lim inf x1 = 0. If lim sup x, = 0, then lim, +m x1 = 0. Since I+m f-m lim f ,m uj = 0 for all j, and ui - x, + xi-, > 0 in W,, it follows that for all i>l, and ~LI(u~-xl+1,u~-x2+x~~~~,~1,0~ as lim (,,xi=O
244
A. S. ELKHADER
t -+m. Since k&1,0) > 1, it is clear from system (4.2) that x1 + 0 is impossible unless x,(t) is identically zero. In particular, x,(O) = 0, a contradiction. Now suppose that limsupx,(t) > 0. Then we can choose a sequence {t,)
such that x.,(t,) = 0, ‘lim”, --lmt, =y and lim, -rm x,(t,) = 0. As above, this implies that lim,,, x;(t,,) = 0 for all i. By arguing as above, we obtain a contradiction. Thus, liminf x,(t) > 0. Now, suppose liminf x,(t) > 0 for n+m t-m l 0, it follows that lim,,,x,(t)=O for all r b m+l.Since liminfxi>Ofori~m,thereexistsapoint(y;”,...,y~,O,...,0; I’m 0 ,..., 0) E W,, yZY’> 0, 1 Q i Q m, in a(E). But the orbit through this point tends to i,,,, so l?,,, E a(%). Now we consider the following systems for l 0. We can choose a sequence {t,} such that t-m lim n-m t n =w, i.m+l(t,) = 0, and lim n_m x(t,> = 0. As above, we obtain for all r>m+l. Since lim,,,xi(t)>O for iE n(Z), x,” > 0. By arguing as above, we obtain a contradiction. Hence, liminf xm+l(t> > 0. t-m
W
Proof of Theorem 4.5. Theorem 4.4 implies that n(X) contains a point (a; u), Pi > 0, 1 < i < n. Since every point in R(K) must have ci = 0, 1 < i < n + 1, it follows that (f, 0) E R. The trajectory x(t) with initial value x(0) = f approaches E, by Proposition 4.2. Therefore, I?;, E Cl(%). But k, is asymptotically stable; then n(x) = i,,. Thus, I!?‘, is globally asymptotically stable. n
This work formed part of my doctoral thesis at Arizona State University. I wish to express my gratitude to Professor Hal L. Smith and the referees for their helpful suggestions. REFERENCES 1 2 3
H. I. Freedman
and P. Waltman,
Persistence
in models
of three
interacting
preda-
tor-prey populations, Math. Eiosci. 68:213-231 (1984). M. W. Hirsch, Systems of differential equations which are competitive or cooperative. I. Limit set, SIAMJ. Math. Anal. 13:167-189 (1982). M. W. Hirsch, Systems of differential equations which are competitive or cooperative. II. Convergence almost everywhere, SIAM J. Math. Anal. 16:432-439 (1984).
4
L. Markus, Asymptotically autoqomous differential systems, in Contributions to the Theory of Nonlinear OsciIation, Vol. 3, Princeton Univ. Press, Princeton, N.J., 1956, pp. 17-29.
5
S. J. Pirt, Principles of Microbe and Cell Cukation, Blackwell Scientific, Oxford, 1975. G. E. Powell, Equalization of specific growth rates for synthrophic association in batch culture,/. Chem. Technol. Biotechno[. 34:97-100 (1984).
6 7 8 9 10 11
G. E. Powell, Stable coexistence of synthrophic associations in continuous culture, J. Chem. Technol. Biotechnol. 35B: 46-50 (1985). G. E. Powell, Stable coexistence of synthrophic chains in continuous culture, J. Theor. Pop. Biol. 30:17-25 (1986). J. Smillie, Competitive and cooperative tridiagonal systems of differential equations, SIAM .I. Math. Anal. 15530-534 (1984). H. L. Smith, Systems of ordinary differential equations which generate an order preserving flow. A survey of results, SIAM Rec. 30:87-113 (1988). H. L. Smith, Competing subcommunities of mutualists and a generalized Kamke theorem, SZAM J. Appl. Math. 46:856-874 (1988).
