9090
Ind. Eng. Chem. Res. 2008, 47, 9090–9098
A Rationale for Monod′s Biochemical Growth Kinetics Doraiswami Ramkrishna* and Hyun-Seob Song* School of Chemical Engineering, Purdue UniVersity, West Lafayette, Indiana 47907
A rationale is presented for Monod′s biochemical kinetics which has hitherto been viewed as purely empirical. It is shown here that this expression, which continues to be used by modelers of biochemical processes, has indeed an origin that relates to the detailed treatment of metabolism. This rationale originates from the resolution of the metabolic network into so-called elementary modes and assumes steady state for every metabolite in each elementary mode. 1. Introduction
2. Monod Kinetics
The authors are pleased to contribute to a special issue in honor of a distinguished colleague and friend whose contributions to chemical reaction engineering have enriched the subject through basic ideas and mathematical insight. It is then fitting to address in a festschrift for Arvind Varma the fundamental background of a kinetic expression that has been used in biochemical engineering and bioreactor analysis with no more than a rationale of its empirical appropriateness and a vague appeal to enzyme substrate kinetics. The objective of this paper is to dispel this notion by elucidating its background in the light of numerous developments of the past two decades in modeling metabolism in detail. Chemical reaction engineers have pursued other similar quests in providing support for mathematical abstraction of complex systems in various ways. Thus the subject of lumping of multicomponent, multireaction systems in petrochemical reaction engineering sought to find component clumps whose concentrations would serve to define their rates of change and to see the circumstances under which such simplifications would emerge. Naturally, linear reaction systems could be dealt with more elegance as, for example, in the publication of Wei and Kuo.1 Several other papers exist on the concept of lumping of nonlinear reaction systems. Metabolism involves a very complicated reaction network with multitudes of reactions enabled by the catalytic activities of a great many enzymes. It involves the catabolic breakdown of substrates into diverse intracellular metabolites that are recruited into the synthesis of biomass and many fermentation products that are released into the abiotic environment. That such a complex process can find a reasonable description with a model that can at best be compared to a single reaction catalyzed by an enzyme deserves reflection toward providing some insight into the extent of success that it has enjoyed. It is fair to warn that an explanation of this kind cannot be unique but insofar as it is based on the many insights that have become available in treating metabolism, it has all the attributes of a reasonable argument. The treatment begins with a short prologue on Monod kinetics, followed by a stoichiometric formulation of metabolism, its decomposition into so-called elementary modes, the pseudostate assumption for intracellular metabolites, and finally toward the formulation of Monod kinetics by determining the Michaelis constant and the growth rate of the organism that would fit in the expression.
The kinetic expression due to Monod applies to the straitjacketing of metabolism into the following single step
* To whom correspondence should be addressed: E-mail: ramkrish@ ecn.purdue.edu;
[email protected].
B + S f (1 + Y)B + · · · (1) Thus the foregoing reaction simply describes the yield of Y mass units of biomass from (a carbon, rate limiting) substrate of course, however, with the aid of existing biomass. The Monod kinetic expression is designed to address the rate at which transformation (1) occurs without regard to the actual mechanism by which metabolism functions to accomplish the same. The specific growth rate (i.e., the rate of synthesis of biomass per cell) is given by µs (2) K+s where µ is the maximum growth rate attained when the substrate concentration (defined per unit volume of culture, i.e., suspension of cells) is sufficiently larger than the Michaelis constant K. We use the term Michaelis constant somewhat loosely because its usage is reserved for enzyme-catalyzed reactions in which it can be expressed in terms of rate constants associated with elementary steps. In the use of Monod kinetics, it is more customary to describe K as a “half saturation” constant which represents the concentration of substrate at which the growth rate is half the maximum possible rate. For substrate concentrations of magnitude much less than that of K, the growth rate increases linearly with s. This then is the kinetic expression due to Monod which does reasonably well during the exponential phase of batch growth. Data by Herbert et al.2 are included in Figure 1 which shows some early data on growth of Aerobacter cloacae fitted to Monod growth kinetics. The fact that its capacity to describe phenomena at low substrate concentrations is greatly diminished is suppressed in addressing rG )
Figure 1. Data from Herbert et al. (1956) in support of Monod kinetics.
10.1021/ie800905d CCC: $40.75 2008 American Chemical Society Published on Web 11/08/2008
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9091
Figure 2. Comparison of growth rate curves: The solid line in black indicates the reference growth rate rG which is given by ∑k3 ) 1 µk/(Kk + s), and the solid line in gray and dotted line in black represent the approximated growth rate rG,app which is given by µ/(Kalt,j + s) (j ) 1, 2, 3) (solid in gray) or µ/(Kopt + s) (dotted in black) where µ ) ∑k3 ) 1 µk while Kalt and Kopt are determined from eq 24 and 31 with m ) 3, respectively.
batch growth but readily exposed by its failure to make predictions in continuous cultures at lower dilution rates. The maximum dilution rate (above which washout of cells from the reactor occurs) is, however, predicted by the Monod model with some reasonableness. There have been many discussions in the literature seeking to improve the prediction of growth dynamics by the Monod model with other growth laws but are not of significance to this paper as we are more concerned with deriving Monod′s model as an approximation from more detailed view of metabolism as a whole. 3. Stoichiometric Formulation of Metabolism
n
ij
j)1
j
ir r
∑r
s,i
(4)
i∈U
Letting mj, j ) 1, 2, ..., n, denote the intracellular metabolite levels (mass fractions relative to the total biomass), it is usual to make the pseudosteady state approximation for each of the metabolites which gives
∑ V r - r m ) 0, ij i
G
j
j ) 1, 2,..., n
(5)
i)1
The second term on the left-hand side of eq 5 is the dilution term due to biomass growth whose origin was pointed out by Fredrickson.3 The term is generally negligible for individual metabolites so that a satisfactory representation of eq 5 is R
p
∑ V M +∑ π P ) 0,
rs )
R
We will limit ourselves to central carbon metabolism in which the carbon source is admitted into the metabolic network in many different ways. Metabolism involves uptake of the substrate followed by numerous intracellular reactions. We assume that there are a total of R reactions numbered from 1 to R and represented by -σiS +
by U the set of indices for reactions in which substrate uptake occurs. Denoting the consumption of substrate in the ith reaction by rs,i and the total consumption by rs, we see that
i ) 1, 2, ..., R
(3)
r)1
In eq 3, Mj, j ) 1, 2, ..., n, represent the intracellular metabolites, and Pr, r ) 1, 2, ..., p, represent the fermentation products released into the abiotic environment. The stoichiometric constant σi for the substrate S is generally non-negative as we envisage only uptake of the substrate and forbid for our purposes any circumstance of it also arising as a fermentation product. For intracellular reactions, however, σi ) 0. The stoichiometric coefficients νij can assume both positive (when Mj is a product) and negative (when it is a reactant). The stoichiometric coefficients πir of the fermentation products Pr will be regarded as positive disregarding the possibility that some could reenter the cell as a substrate. That some or all of the intracellular reactions can be reversible is already accommodated by the allowable signs of the stoichiometric coefficients. We denote
∑ V r ) 0, ij i
j ) 1, 2,..., n
(6)
i)1
Thus the reaction rates, ri, i ) 1, 2, ..., R or more commonly referred to in the biological literature as “fluxes” form a row vector rT ≡ (r1, r2, ..., rn) which must be in the null space of the matrix N ≡ {νij; i ) 1, 2, ..., R; j ) 1, 2, ..., n} since eq 6 implies that rTN ) 0. Alternatively, the column vector r may also be viewed as in the null space of NT, as NTr ) 0. At this point the flux balance approach of Palsson and co-workers (see Edwards et al.4 for example) is to determine vectors in the null space of NT that would satisfy the criterion of maximizing an objective function such as the biomass yield and determine the flux vector up to a multiplicative constant. A measurement of the total substrate uptake rate rs would then provide the unique flux vector of interest. Since this would of course yield the biomass synthesis rate, the problem of providing a rationale for
9092 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
Monod kinetics boils down to showing that the sum of many enzyme-catalyzed substrate uptake rates can be represented by a single enzyme substrate kinetic expression. Although this is basically the approach of rationalization of this paper, we postpone this specific step to a subsequent stage in the interest of expounding a more detailed view of metabolism brought about by the decomposition of the metabolic network into elementary modes. 4. Decomposition of Metabolic Network into Elementary Modes A particularly interesting way to view a metabolic network appears to be in terms of elementary modes.5,6 An elementary mode is a part of the network which involves the uptake of substrate, leading it through a sequence of reactions ultimately producing one or more metabolic products and precursors for biomass synthesis. Thus the flux vector has only positive components in an elementary mode. The set of elementary modes forms a convex set and will form a convex basis of the null space of the matrix NT if all exchange fluxes are irreversible or the convex cone (NTr ) 0; ri g 0, i ) 1, 2, ..., n). The problem of identifying the elementary modes for an arbitrary metabolic network has been thoroughly investigated by Schuster et al.7 Elementary modes can be very substantial in number. Reversible reactions are accommodated in separate elementary modes so that the “backward” reaction rate is also entered as a positive component in a mode different from that which has the forward reaction rate. We assume that there are m elementary modes in all. Each mode represents a vector of non-negative fluxes (reaction rates) to be determined uniquely only when the substrate uptake rate through that mode is known. Given that the reactions have also been identified for each mode, we may represent each mode also by the reactions involved as follows. n
-σiS +
∑
p
VijMj +
j)1
∑ π P ) 0,
i ∈ Uk,
ir r
k ) 1, 2, ..., m
r)1
(7) where {Uk}km) 1 are subsets of {1, 2, ..., R} with Uk representing the selection of reactions from the network for the kth mode. Since the different modes exhaust all reactions in metabolism, m
we have ∪ Uk ) U. The steady state balance of reactions in the k)1 kth mode is given by replacing eq6 by
∑ V r ) 0, ij i
j ) 1, 2, ..., n,
k ) 1, 2, ..., m
(8)
i∈Uk
Clearly in eq8 we may set ri ) 0 for any reaction not included in the mode. What distinguishes eq 8 from eq 6 is the fact that eq 8 has only one linearly independent solution for each k, which we denote by rk ≡ [r1k, r2k, ..., rRk ]; k ) 1, 2, ..., m. It will be understood that the first component r1k in the vector represents the substrate uptake rate through that mode. Since it is only the total uptake rate of substrate that is known, the uptake rate through a given mode is known only up to a multiplicative constant. Thus, the total uptake rate of substrate is then given by
vector by r1T ≡ [r11, r12, ..., r1m], so that the entire reaction rate vector r may be represented by r ) Zr1
rG ) hTZr1
m
∑w r
k k 1
(9)
(11)
From eq 9, we can describe the total uptake rate of substrate by rs ) wTr1, where the row vector wT ≡ [w1, w2, ..., wm]. The yield coefficient Y defined earlier in representing growth from unit mass of the carbon substrate is given by wTr1 rs )Y) T rG h Zr1
(12)
Expression 12 shows substrate dependence of the yield coefficient which is acknowledged although not an original aspect of Monod′s model. 5. View of Metabolism through Elementary Modes The viewpoint of metabolism that emerges from the foregoing perspective consists in envisaging substrate uptake through each of the elementary modes. Denoting by wk the fraction (of the total substrate uptake rate) through the kth elementary mode, it is now possible to calculate all intracellular and extracellular fluxes connected with fermentation products. The weight fractions for the different modes are governed by the phenomenon of metabolic regulation. Their calculation would in fact require a detailed treatment of all the reactions in the cell with due account of regulatory processes, clearly a problem of inspiring magnitude. However, if we take refuge in the pseudosteady state hypothesis, the uptake rate through each mode could be stoichiometrically related to the intracellular fluxes through that mode. The uptake rates through the different modes have then been modeled by Michaelis-Menten kinetics with constants(µ′k,Kk),k ) 1,2, ..., m (see e.g., Provost et al.9). The implication of this approach is that uptake can be modeled by single substrate enzyme kinetics that could be varied for each of the modes to account for the regulatory effects in each mode. The total uptake rate of the substrate denoted rs can then be written as m
rs )
µks
∑ K +s k)1
k
The growth rate rG may then be written as m
rG )
µks
∑ K +s k)1
rs )
(10)
The matrix Z in the equation above arises out of expanding the reaction rate vector r in terms of the elementary mode flux vectors and representing each of the latter in terms of the uptake flux through the mode. Thus we have the entire reaction rate vector depend only on the uptake flux vector r1 through the assumption of steady state for the intracellular components. The growth rate of biomass is related to the intracellular fluxes and has been represented in the biological literature (see Stephanopoulos et al.8) by rG ) hTr, where the row vector hT is available. Thus the growth rate may be expressed in terms of the uptake vector as
k
µk ≡
∑hz µ
i ik
k
(13)
i
where zik is the (i,k)th coefficient of the matrix Z.
k)1
where wk values are positiVe weights. The total uptake of substrate is thus viewed to be distributed among various elementary modes. Furthermore, we may define the uptake (row)
6. Approximation of Growth Rate by Monod Kinetics We wish to approximate ∑m k ) 1 µks/(Kk + s) by the expression µs/(K + s) and seek the best choice of µ and K given the values
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9093
of {(µ1,K1), (µ2,K2), ..., (µm,Km)} for each of the MichaelisMenten kinetic expressions for the uptake reactions for the different modes. Alternatively, we seek to minimize the residual
constant for the Monod model. For m ) 4, by comparing coefficients of s0, s1, s2, and s3, we respectively have K)
m
µks µs µs ) Rm(s) ≡ rG (K + s) k)1 (Kk + s) (K + s)
∑
As a measure of the quality of the approximation, we consider on the space of functions defined on both the metrics | |∞ and | |2 defined by (1) |f|∞ ≡ sup |f(x)| (2)
|f|22 ≡
0exes s
∫ f(x)
2
0
(14)
m
∑µ
k
(15)
k)1
The Michaelis constant for the Monod model is then to be determined so that the residual is minimal in some sense. Toward this end, we consider below the case of m ) 3 for simplicity before stating the general case. It is convenient to begin by exploring the behavior of f3(s) ≡ (K1+s)(K2+s)(K3+s)(K+s)R3(s)/s which is given by µ1(K + s)(K2 + s)(K3 + s) + µ2(K + s)(K3 + s)(K1 + s) + µ1(K + s)(K1 + s)(K2 + s) - µ(K1 + s)(K2 + s)(K3 + s) ) [K(µ1K2K3 + µ2K3K1 + µ3K1K2) - µK1K2K3] + s[µ1(K2K3 + KK2 + KK3) + µ2(K3K1 + KK3 + KK1) + µ3(K1K2 + KK1 + KK2)-µ(K1K2 + K2K3 + K3K1)] + s2[µ1(K + K2 + K3) + µ2(K + K3 + K1) + µ3(K + K1 + K2) µ(K1 + K2 + K3)] (16) From eq 15, setting m ) 3, we have µ ) µ1 + µ2 + µ3. Inserting this into eq 16 and equating each of the square brackets equal to zero, we have three different alternates for the value of K from comparing coefficients of s0, s1, and s2. K)
µ1 + µ2 + µ3 -1 -1 µ1K-1 1 + µ2K2 + µ3K3
µ1K1 + µ2K2 + µ3K3 K) µ1 + µ2 + µ3
K)
k)1 4
l
l)1,k 4
(22)
∑µ ∑ K k
k)1
l
l)1,k
4
∑µ K
k k
K)
k)1 4
(23)
∑µ
k
k)1
The notation l ) 1, k is designed to exclude the running index l from taking on the value of k. As in the case of m ) 3, we find that if all the Michaelis constants are equal, then the formulas 20, 21, 22 and 23 all consistently yield the same equivalent Michaelis constant. We now return to the general case of m modes. Toward this end, it is of interest to consider all combinations of r elements of the set Sm ≡ {1, 2, ..., m} and denote this set by Cr(Sm). Identify the combinatory elements of this set with a discrete index, say i to associate uniquely with m a combination (i1, i2, ..., ir). Clearly there are elements in r the set Cr(Sm). This notation can be conveniently used to express certain functions of the combinations. For example, the term r ∑ i∈Cr(Sm) ∏k)1 Kik represents the sum of all possible products of r distinct Michaelis constants obtained from the original set of m constants. It is necessary to introduce one further notation. A specific element, say k may be omitted from the set Sm, which may be denoted by Sm,k before seeking r combinations of them. With this notation, it is now possible to return to the general case of m elementary modes. We note that fm(s) contains coefficients, say R0, R1, ..., Rm-1 of s0, s1, s2, ..., sm-1, respectively, each of which can be equated to zero to produce m alternates for the equivalent Michaelis constant. These can be seen to be respectively
()
m
∑µ
k
K)
k)1 m
(j ) 0),
∑
µkK-1 k
∑
µkKk
k)1 m
K)
∏
i∈Cm-j-1(Sm,k) r)1
m
m-j-1
∑µ k)1
(19)
∑
m-j-1
k)1
k
Of course eq17, 18, and 19 cannot be simultaneously true except for when K1 ) K2 ) K3. It should be intuitively clear at this stage that if K1, K2, K3 are not very disparate, one of the foregoing candidates could serve to be a reasonable Michaelis
4
∑µ K ∑ K k k
(17)
µ1K1(K2 + K3) + µ2K2(K3 + K1) + µ3K3(K1 + K2) (18) K) µ1(K2 + K3) + µ2(K3 + K1) + µ3(K1 + K2)
(20)
K ) [µ1K1(K2K3 + K2K4 + K3K4) + µ2K2(K1K3 + K1K4 + K3K4)+ µ3K3(K1K2 + K1K4 + K2K4) + µ4K4(K1K2 + K1K3 + K2K3)] ⁄ [µ1(K2K3 + K2K4 + K3K4) + µ2(K1K3 + K1K4 + K3K4) + µ3(K1K2 + K1K4 + K2K4) + µ4(K1K2 + K1K3 + K2K3)] (21) 4
dx
The residual is represented by the vector Rm ≡ {Rm(x): 0 e x e s}. The upper limit s can of course be increased indefinitely. The two norms in eq 14 for the residual represent different measures. The first is less forgiving of large deviations even in a small localized zone than the second. The second measure is perhaps sufficient but the first has the advantage of producing various candidates for the effective Michaelis constant for Monod kinetics. We shall consider both the measures in eq14. Let us address at first the first metric in (14). Since it is concerned with deviation at each point s, for suitably large values of s (considerably larger than the largest Michaelis constant in the set {Kk:k ) 1,2, ..., m}), the growth rate rG is given by the sum of the maximum rates from each of the modes so that µ)
µ1 + µ2 + µ3 + µ4 -1 -1 -1 µ1K1 + µ2K-1 2 + µ3K3 + µ4K4
∑
∏
i∈Cm-j-1(Sm,k) r)1
Kir (j ) 1, 2, ..., m - 2), Kir
m
∑µ K
k k
K)
k)1 m
∑µ
k
k)1
(j ) m - 1)
(24)
9094 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
Again, as before, if the different Michaelis constants in the original set are all the same, this common value is assumed by all the alternates in eq 24 for the equivalent Michaelis constant. Equation 24 may be viewed as providing various averaging rationales for the calculation of the half-saturation constant using the uptake kinetic constants {µk,Kk}. We now address the issue focal to this paper, namely, that we seek the circumstances under which Monod′s growth kinetics can be reasonably applied to microbial growth. That we can get an exact fit of the Monod model for the case where all the Michaelis constants are equal has already become obvious. One can then surmise that small fractional variations among the Michaelis constants would preserve the foregoing conclusion. We investigate below the various ways in which the applicability of Monod kinetics can break down. The alternates (24) for the equivalent Michaelis constant show that the constants µk and Kk are both involved in an answer to the question. The “L-infinity norm” of the residual Rm, given by |Rm |∞ ) sup |Rm(x)| 0exes
A coarse estimate of the upper bound for the foregoing metric is readily obtained as follows. We rewrite Rm(x) as m-1
Rm(x) )
∑ R g (x),
x
k)0
m
(K + x)
∏ (K + x)
(25)
2
m KKk k µjµk ln 2 + µ2k Kk(1 + 2 ln Kk) Kj (K K ) k j Kj k)1 j)1,k k)1 m
|Rm |22 )
m
∑∑
∑
m
+ µ2K(1 + 2 ln K) + 2
µµk
[ ]
∑ (K - K) ln KK k)1
k
K2
K2k k
(30)
which was obtained by demanding boundedness that eliminates two s - dependent terms, one as s and the other as ln s (other terms containing s also exist that automatically vanish at infinity). Interestingly, eq 15 is recovered again from this exercise. For the optimal choice of the equivalent Michaelis constant, we seek to minimize the residual with respect to K. We obtain on equating to zero the derivative of the residual metric in eq 30 the following algebraic equation for K m
2
k+1
gk(x) ≡
k k
from the more elaborate metabolic model. We will examine the different alternates for the equivalent Michaelis constant. We now consider the mean square deviation (or L-2 norm) of the residual |Rm|2 as defined in eq 14. The rational nature of Monod kinetics allows an exact calculation of the foregoing metric for the residual as
µµk
∑ (K - K) k)1
k
[
]
K2 ln K - K2k ln Kk + K(2 ln K + 1) + (Kk - K) m
µ2(3 + 2 ln K) ) 0, µ ≡
r
r)1
∑µ
k
(31)
k)1
From the generalization of the triangular inequality, we may write m-1
|Rm |∞ e
∑ |R ||g | k
(26)
k ∞
k)0
It is not difficult to show that |gk |∞ ) gk(x*k )
(27)
where x*k ’s must satisfy the equations m
k+1)
∑ x(K + x)
-1
r
+ x(K + x)-1,
k ) 0, 1, ..., m - 1
r)1
Note that eq 26 and 27 hold for arbitrary s. We denote the upper bound for |Rm|∞ suggested by eq 26 by m-1
|Rm |∞ ≡
∑ |R ||g | k
(28)
k ∞
k)0
Thus we have an estimate of the deviation of Monod kinetics in the sense of the above metric. The task that remains is to choose one of the alternates (i.e., an equivalent K for which we have the least value of the right-hand side of the inequality over a suitable range of substrate concentrations). If the first of eq24 is used for the value of K, the fit will be good for the lowest substrate concentration range. This is because the slope of the curve at the origin is given by m
∑µ
k
k)1
[ ]
d s ds Kk + s
m
)s)0
∑µ
k
k)1
() Kk
K2k
m
)-
∑µ K k)1
µ -1 k k )-
K (29)
The expression to the extreme right is readily seen to be the slope of the curve of the Monod growth expression at s ) 0. However, for larger values of s, there is no guarantee that this Monod growth expression will stay close to the growth rate
The solution of eq 31 serves to evaluate the different alternates for the equivalent Michaelis constant obtained from eq 24. For our demonstration, we use various “profiles” of Michaelis-Menten kinetics {µk,Kk; k ) 1, 2, ..., m}. It is of interest to characterize each profile with respect to how the µk is paired with Kk. The case of identical values of Kk is no longer of interest as it has already been seen to imply equivalence to Monod kinetics. Thus we consider pairs that are positively “correlated” by which is meant that larger values of µk are paired with larger values of Kk as well as pairs that are negatively correlated in which larger values of µk are paired with smaller values of Kk. In each case, we evaluate the optimal Michaelis constant, the various alternates in eq 24, the upper bound of the residual metric |Rm|∞. The evaluation of Monod kinetics is done as follows. We first evaluate the various alternates of eq 24 for the equivalent Michaelis constant by comparing their values with the optimal K obtained by minimizing the mean square residual. Next, we compare the mean square deviation of the “approximated growth rates using Monod kinetics” (rG,app) from the “reference growth rate using the substrate uptake rates through the different elementary modes” (rG). To assess the quality of the fit of the Monod growth rate relative to the reference growth rate, we calculate the actual fractional deviation of the former from the latter. 7. Discussion of Results Tables 1 and 2 each provide four different choices of kinetic profiles of the uptake rates for our demonstration. The four different profiles cover kinetic parameter choices for µk and Kk with positive correlation (case “b”), negative correlation (case “c”), and neutral correlation (cases “a” and “d”). Cases “a” and “d” feature different distributions of the Michaelis constants for the same maximum reaction rates. For example, in Table 1, one mode is featured with a Michaelis constant much larger than the other two in case “a”, while in case “d” the large Michaelis constant is as much larger than the other two. In Table
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9095 Table 1. Parameter Values Used for Plotting the Curves Shown in Figure 1 and the Residual Normsa figures
{µ1, µ2, µ3}
{K1, K2, K3}
Kalt () alternates of K) |R3|∞(K ) Kalt) |R3|2(K ) Kalt)
Kopt () optimal K) |R3|∞(K ) Kopt) |R3|2(K ) Kopt)
Figure 2a
{1/3, 1/3, 1/3}
{0.1, 2, 10}
{0.283, 1.752, 4.033} {0.423, 0.243, 0.315} {1.210, 0.497, 0.495}
2.788 0.291 0.333
Figure 2b
{0.1, 0.2, 0.7}
{0.1, 2, 10}
{0.855, 4.021, 7.410} {0.488, 0.154, 0.132} {1.807, 0.664, 0.258}
6.611 0.139 0.194
Figure 2c
{0.7, 0.2, 0.1}
{0.1, 2, 10}
{0.141, 0.657, 1.470} {0.199, 0.288, 0.442} {0.486, 0.266, 0.515}
0.618 0.289 0.265
Figure 2d
{1/3, 1/3, 1/3}
{1, 2.5, 5}
{1.875, 2.353, 2.833} {0.070, 0.043, 0.065} {0.279, 0.109, 0.111}
2.583 0.055 0.070
a Three elementary modes are considered. |R3|∞ and |R3|2 indicate the upper bound of the L-infinity norm of residual, and the L-2 norm, respectively, for three elementary modes.
Table 2. Parameter Values Used for Plotting the Curves Shown in Figure 2 and the Residual Normsa figures Figure 3a
Figure 3b
Figure 3c
Figure 3d
{µ1, µ2, ..., µ10} {1/10, 1/10, ..., 1/10}
{0.0010, 0.0020, 0.0039, 0.0078, 0.0156, 0.0313, 0.0626, 0.1251, 0.2502, 0.5005}
{0.5005, 0.2502, 0.1251, 0.0626, 0.0313, 0.0156, 0.0078, 0.0039, 0.0020, 0.0010}
{1/10, 1/10, ..., 1/10}
{K1, K2, ..., K10} {0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2}
{0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2}
{0.1, 0.2, 0.4, 0.8, 1.6, 3.2, 6.4, 12.8, 25.6, 51.2}
{0.1, 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 51.2}
Kalt () alternates of K) |R10|∞ (K ) Kalt) |R10|2 (K ) Kalt)
Kopt () optimalK) |R10|∞ (K ) Kopt) |R10|2 (K ) Kopt)
{0.709, 1.311, 1.5808, 1.742, 1.861, 1.962, 2.062, 2.178, 2.338, 2.610} { 0.684, 0.571, 0.476, 0.408, 0.372, 0.372, 0.408, 0.476, 0.571, 0.684} { 1.770, 1.676, 1.553, 1.398, 1.208, 0.988, 0.768, 0.641, 0.753, 1.081}
5.6406 0.486 0.640
{10.230, 11.311, 12.618, 14.213, 16.168, 18.571, 21.511, 25.069, 29.293, 34.167} {0.386, 0.342, 0.295, 0.247, 0.198, 0.153, 0.116, 0.095, 0.097, 0.129} {2.686, 2.506, 2.299, 2.059, 1.782, 1.464, 1.104, 0.711, 0.349, 0.423}
31.110
{0.150, 0.175, 0.204, 0.238, 0.276, 0.317, 0.360, 0.406, 0.453, 0.501}
0.252
{0.129, 0.097, 0.095, 0.116, 0.153, 0.198, 0.247, 0.295, 0.342, 0.386} {0.147, 0.119, 0.094, 0.077, 0.080, 0.101, 0.132, 0.167, 0.203, 0.239}
0.130
{0.713, 1.335, 1.618, 1.797, 1.939, 2.078, 2.247, 2.509, 3.097, 6.730} {0.527, 0.254, 0.186, 0.160, 0.151, 0.151, 0.161, 0.188, 0.261, 0.585} {1.087, 0.789, 0.684, 0.625, 0.583, 0.546, 0.506, 0.455, 0.400, 0.947}
3.206
0.110 0.294
0.076
0.275 0.399
a Ten elementary modes are considered. |R10|∞ and |R10|2 indicate the upper bound of the L-infinity norm of residual, and the L-2 norm, respectively, for 10 elementary modes.
2, case “d” features one mode with an even larger Michaelis constant. Each row of Tables 1 and 2 is divided into three subrows in the third and fourth columns. In the third column, among three subrows, the top shows Kalt,j (j ) 1, 2, ..., m) (the
alternates to equivalent Michaelis constant K) and the middle and bottom show |Rm|∞ (the upper bound of L-infinity norm of the residual) and |Rm|2 (the L-2 norm of the residual), respectively, when K ) Kalt, while, in the end column, the top
9096 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
Figure 3. Comparison of growth rate curves: The solid line in black indicates the reference growth rate rG which is given by ∑k10) 1 µk/(Kk + s), and the solid line in gray and dotted line in black represent the approximated growth rate rG,app which is given by µ/(Kalt,j + s) (j ) 1, 2, ...,10) (solid in gray) or µ/(Kopt + s) (dotted in black) where µ ) ∑k10) 1 µk while Kaltand Kopt are determined from eq 24 and 31 with m ) 10, respectively.
shows Kopt (the optimum choice of K) and the middle and bottom show |Rm|∞ and |Rm|2, respectively, when K ) Kopt. Figures 2 and 3 show the results of our evaluation of Kalt,j (j ) 1, 2, ..., m) and Kopt for the half-saturation constant in Monod kinetics for various types of kinetic profiles of the uptake rates. The number of elementary modes “m” was assumed to be “3” in Figure 2 and “10” in Figure 3, and the reference growth curves (solid line in black, rG) were obtained from eq 13 with the sets of µk and Kk as given in Tables 1 and 2, respectively. The various Monod kinetics approximating rG are denoted by rG,app where the equivalent Michaelis constant is given by Kalt,j (j ) 1, 2, ..., m) (solid line in gray) or Kopt (dotted line). Figure 2a shows results for a “neutral” kinetic profile displaying a spread of Michaelis constants with change in 1 order of magnitude from the smallest to the largest. Comparison of Kalt,j (j ) 1, 2, 3) with Kopt shows that the best K is chosen differently depending on the criterion, that is, Kalt,2 (i.e., 1.752) is the best in the sense of minimizing |R3|∞ while Kopt is the best regarding the minimum |R3|2 as expected (see Table 1). This observation holds good for the rest of cases considered in Tables 1 and 2. A more rigorous comparison is possible by investigating the fractional deviation of the Monod kinetics rG,app to the reference growth rate rG as shown in Figures 4 and 5. As can be seen in Figure 4a, Kalt,2 (i.e., gray line in the middle) exhibits less deviation than Kopt (dotted) in the very low range of s, but shows larger deviation on average. Figures 2a,b, and 4a,b reflect the strain in the optimal Monod approximation (dotted Kopt line) to the reference growth rate represented by the solid line. The reason of course is the disparate nature of the Michaelis constants. This observation is in striking contrast to the close proximity between the Monod (dotted) and the reference growth rates in Figure 2c,d. Clearly, the reason for Figure 2d to provide excellent fits of Monod kinetics to the reference curve lies in the limited variation in Kk’s. Curiously, Figure 2c shows a good quality fit of Monod kinetics while Figure 2b shows a poor fit with the same variation of Kk′s. This is a consequence of the dependence of the fit on the distribution
of µk′s that result from the different uptake rates. In fact, the negative correlation between µk′s and Kk′s will render one specific mode with the maximum µk and minimum Kk to be the most dominant while others are less important. Then, the growth kinetics is governed predominantly by that of one individual mode. Obviously, this is a situation where Monod kinetics fits well. On the other hand, the positive correlation tends to make the characteristics of individual uptake curves differ from one another and consequently, the effect of disparateness in the Michaelis constants may become even more pronounced. Now, we turn to Figures 3 and 5 where the number of elementary modes was assumed to be 10. The sequence of Kk′s in Figure 3a-c was generated by the relation of Kk ) 0.1 × 2k-1 (k ) 1, 2, ..., 10), and in Figure 3d the Kk′s were adjusted such that two modes show a significant variance in the Michalies constants while others are similar. Again, Figure 3 panels b and c correspond to positive and negative correlations, respectively, and Figure 3 panels a and d are neutral. Behavior quite similar to that in the previous examples is observed in Figure 3a-c. In particular, the effect of negative correlation is strikingly pronounced in this case as seen in Figures 3c and 5c. Figure 3d shows that Monod kinetics can fit the growth curve well even though Kk′s of a few elementary modes are significantly different from others if most modes have similar Kk′s as in Table 2. Our discussion has not included any assessment of the prediction of fermentation products as, for example, in a wellknown paper by Luedeking and Piret.10 This paper assumes that the rate of fermentation product (specifically lactic acid) may be expressed as the sum of a growth associated term and a nongrowth associated term. The first is assumed to be proportional to the growth rate while the second to the biomass concentration with the constants of proportionality found to be a function of pH from experiments. The concepts in this paper may be related to the foregoing model by considering as being growth associated the set of elementary modes contributing to biomass formation, and those that do not as being nongrowth associated. Mathematically, we establish the link between the
Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008 9097
Figure 4. Fractional deviation of the approximated kinetic profiles rG,appfrom the reference growth curve rG shown in Figure 2. The solid line in gray indicates {µ/(Kalt,j + s) - ∑k3 ) 1 µk/(Kk + s)}/∑k3 ) 1 µk/(Kk + s) (j ) 1, 2, 3) and the dotted line in black represents {µ/(Kopt + s) -∑k3 ) 1 µk/(Kk + s)}/ ∑k3 ) 1 µk/(Kk+s).
Figure 5. Fractional deviation of the approximated kinetic profiles rG,appfrom the reference growth curve rG shown in Figure 3. The solid line in gray indicates {µ/(Kalt,j + s) - ∑k10) 1 µk/(Kk + s)}/∑k10) 1 µk/(Kk + s) (j ) 1, 2, ...,10) and the dotted line in black represents {µ/(Kopt + s) -∑k10) 1 µk/(Kk + s)}/ ∑k10) 1 µk/(Kk+s).
treatment here and the rate of formation of fermentation product as formulated by Luedeking and Piret as follows. Let there be p fermentation products with Sp a (p × n) matrix such that Spr () Sp Zr1) represents the rates of formation of the p fermentation products. This rate of formation can be expressed as the sum of that due to growth and that not due to growth. Further, define matrices PG and PM such that the rates of formation of fermentation products by growth-associated processes and nongrowth associated processes are respectively
given by PGSpZr1 and PMSpZr1. Then the total rate of formation of the ith fermentation product is given by eTi PGSpZr1 for growth associated processes and by eTi PMSpZr1 for nongrowth associated processes. Thus any proportionality of the rate of growthassociated formation of the ith fermentation product to the growth rate depends on the time-invariance of the ratio eTi PGSpZr1/(hTZr1) in which the denominator is the growth rate rG of the cell as represented in eq 11. Similarly, the proportionality of the nongrowth associated formation of the ith
9098 Ind. Eng. Chem. Res., Vol. 47, No. 23, 2008
fermentation product depends on the time invariance of eTi PMSpZr1. The time-invariance of both quantities depends upon the time-invariance of the uptake rate vector r1. For large enough concentration of the substrate, we expect r1 to be time-invariant thus vindicating the Luedeking and Piret approach. However, when the substrate concentration drops to small values, regulatory effects associated with growth and maintenance will no longer sustain the assumptions made by Luedeking and Piret.
predominance of one specific mode with the maximum µk and minimum Kk. Even in the case that a few Kk’s are significantly different from the others, the quality of fit by Monod equation can be good if most of Kk’s have similar values. On the other hand, Monod kinetics may fail when Kk’s show significant variance and µk’s vary independently of Kk or Kk’s have positive correlation with µk’s. Literature Cited
8. Factors Against Monod Kinetics While support for Monod growth kinetics could be obtained from viewing metabolism as being the result of substrate uptake through several elementary modes, it was assumed that the relative distribution of substrate through different modes remains largely invariant. This is in fact not true as one is well aware that at very low substrate concentrations, the organisms regulate their metabolism to consume substrate more for maintenance of viability than for growth. This means that the kinetic profiles are themselves subject to change under different conditions. This situation is manifestly apparent in continuous cultures but virtually invisible in batch growth as change in biomass does not occur toward the end of batch growth. Thus the applicability of Monod kinetics to continuous cultures leads to a constant biomass concentration at steady state regardless of the dilution rate (flow rate/hold-up volume) as long at it is below the maximum dilution rate. Experimental data, on the other hand, show decrease in steady state concentration of biomass presumably due to maintenance effects. Clearly, the complexities of regulatory processes cannot be described by Monod′s growth kinetics. As long as the Michaelis constants are not widely disparate, Monod’s growth kinetics has a domain of application, however. With respect to the applicability of Monod kinetics the scenario of equal Michaelis constants is not very different from that when the Michaelis constants are not widely disparate. When regulatory processes are accounted for, however, even small differences in Michaelis constants lead to rather complex behavior in continuous cultures. Kim11 has recently shown that both steady state multiplicity and oscillatory dynamics can result with the cybernetic models of Ramkrishna and co-workers.12-15 9. Conclusions This paper has demonstrated that a rationale exists for the applicability of Monod′s growth kinetics. It also shows how its failure must be expected as it cannot accept the burden of describing regulatory processes. Monod kinetics becomes attractive when the Michaelis constants are not widely disparate but even when they are the existence of negative correlation can give rise to reasonable Monod approximations to the growth rate. The reason for the latter can be explained by the
(1) Wei, J.; Kuo, J. C. W. A lumping analysis in monomolecular reaction systems: Analysis of the exactly lumpable system. Ind. Eng. Chem. Fundam. 1969, 8, 114. (2) Herbert, D.; Elsworth, R.; Telling, R. C. The continuous culture of bacteriasA theoretical and experimental study. J. Gen. Microbiol. 1956, 14, 601. (3) Fredrickson, A. G. Formulation of structured growth models. Biotechnol. Bioeng. 1976, 18, 1481. (4) Edwards, J. S.; Ibarra, R. U.; Palsson, B. O. In silico predictions of Escherichia coli metabolic capabilities are consistent with experimental data. Nat. Biotechnol. 2001, 19, 125. (5) Schuster, S.; Dandekar, T.; Fell, D. A. Detection of elementary flux modes in biochemical networks: a promising tool for pathway analysis and metabolic engineering. Trends Biotechnol. 1999, 17, 53. (6) Schuster, S.; Fell, D. A.; Dandekar, T. A general definition of metabolic pathways useful for systematic organization and analysis of complex metabolic networks. Nat. Biotechnol. 2000, 18, 326. (7) Schuster, S.; Hilgetag, C; Woods, J. H.; Fell, D. A. Reaction routes in biochemical reaction systems: Algebraic properties, validated calculation procedure and example from nucleotide metabolism. J. Math. Biol. 2002, 45, 153. (8) Stephanopoulos, G.; Aristidou, A. A.; Nielsen, J. Metabolic Engineeringsprinciples and Methodologies; Academic Press: San Diego, CA, 1998. (9) Provost, A.; Bastin, G.; Agathos, S. N.; Schneider, Y.-J. Metabolic design of macroscopic bioreaction models: Application to Chinese hamster ovary cells. Bioprocess Biosyst. Eng. 2006, 29, 349. (10) Luedeking, R.; Piret, E. L. A kinetic study of lactic acid fermentation-batch process at controlled pH. J. Biochem. Microbiol. Technol. Eng. 1959, 1, 393. (11) Kim, J. I. , A hybrid cybernetic modeling of the growth for Escherichia coli in glucose-pyruvate mixtures. Ph.D. Thesis. Purdue University, West Lafayette, IN, 2008. (12) Kompala, D. S.; Ramkrishna, D.; Jansen, N. B.; Tsao, G. T. Investigation of bacterial-growth on mixed substrates - Experimental evaluation of cybernetic models. Biotechnol. Bioeng. 1986, 28, 1044. (13) Young, J. D.; Ramkrishna, D. On the matching and proportional laws of cybernetic models. Biotechnol. Progress 2007, 23, 83. (14) Young, J. D.; Henne, K. L.; Morgan, J. A.; Konopka, A. E.; Ramkrishna, D. Integrating cybernetic modeling with pathway analysis. A dynamic systems level description of metabolic control. Biotechnol. Bioeng. 2008, 100, 542. (15) Kim, J. I.; Varner, J. D.; Ramkrishna, D. A. hybrid model of anaerobic E. coli GJT001: Combination of elementary flux modes and cybernetic variables. Biotechnol. Progress 2008, in press.
ReceiVed for reView June 9, 2008 ReVised manuscript receiVed September 17, 2008 Accepted September 18, 2008 IE800905D
