BIOTECHNOLOGY AND BIOENGINEERING, VOL. XIX, PAGES 1831-1850 (1977)
Practical Identifiability of Growth and Substrate Consumption Models MARKKU NIHTILA and JOUKO VIRKKUNEN, Helsinki University of Technology, Electrical Engineering Department, Control Engineering Laboratory, SF-02150 Espoo 15, Finland Summary The estimation of parameters in several dynamic models, which describe growth and substrate consumption, has been carried out using a modified Gauss-Newtontype method. The four models considered are Monod, Contois, linear specific growth rate, and an enzyme kinetic model. The initial values of the differential equations are included in the parameter vector which will be estimated. The efficiency of the method and the confidence limits of the parameters were studied using simulated measurement noise. The experimental results describe Tn'chodrrma viride growing on glucose as the main carbon source.
INTRODUCTION
The use of mathematical models has become increasingly popular for the description of microbiological processes. The extensive computational effort necessary in using models has become less significant owing to the advent of powerful computers and automatic data logging facilities. In recent years several computerized pilot plants have been developed'-' and there is an urgent need to develop models and computational algorithms for the utilization of the great amounts of information which can be collected during experimental runs. In developing mathematical models the first and basic aim is to give a concise and unambiguous description of the biological mechanism of the process. The different mechanisms and the pertinent models can be distinguished by comparing measured results with the models' predictions. However, the model is most often a group of first order nonlinear differential equations with unknown parameters. The parameters must be estimated by using measured results. A satisfying fit of the model response and experimental results can often be achieved, especially if the number of free parameters is high and the number of measured 1831
@ 1977 by John Wiley & Sons, Inc.
1832
NIHTILA AND VIRKKUNEN
points is comparatively low, as is often the case in microbiological experiments. This good fit can be highly misleading as the confidence limits of the estimated parameters can be so great that the actual numerical values of the parameters have no relevance at all. Unfortunately there are no general methods for the computation of confidence limits for parameters in nonlinear estimation problems. Moreover, the parameters are normally interdependent and the estimated values can be biased. This bias is practically impossible to find by mathematical methods. Simple graphical means such as the Lineweaver-Burk plot can only be used in simple cases. The parameters in more complicated models must be estimated by using rather sophisticated mathematical methods which are based on the minimization of a criterion of the goodness of fit. In nonlinear problems, however, there can be several local minima and the algorithm can give erroneous results by leading to a secondary minimum of the criteria. In the opinion of the present authors, as stated also, e.g., by S c h ~ e p p e ,the ~ most practical and feasible way of testing the reliability of the parameters in nonlinear models is to simulate the measurements using either simulated values for measurements or using real measurements and adding to them simulated noise. A study of this type is reported here. The first section briefly describes the materials and methods used. Model structure is considered in the second section by dividing the problem into three parts: the system model, the output model, and the performance index. Mathematical model specifications which are used in parameter estimation are presented in detail in the third section. The estimation method is then qualitatively explained. (Mathematical details of the method are given in Appendix B.) Results of estimation, i.e., parameter values, confidence limits, and residual sums with relevant conclusions, are presented in the next sections. Finally, a simulation study of the method, performed using Monod’s model with feasible parameter values, is described. MATERIALS AND METHODS
The microorganism used was Trichoderma viride strain QM 9414, obtained from the U.S. Army Natick Laboratories. The inoculum (100 ml) was four days old, grown on 1% Solca Floc medium. The fermentor was a New Brunswick Magnaferm, working volume 10 liter. The medium contained in dliter: (NHk),SO,, 2.8; KH2P04, 15.0; MgSO,, 0.6; CaC12,0.6; urea, 0.6; peptone, 1.5.
GROWTH AND SUBSTRATE CONSUMFTION MODELS
I833
The carbon source was a glucose syrup (n% dry weight) containing 20% glucose. Tween 80 (0.2% v/v) and trace elements (Fe", Mn++, Zn"", C o " ) were also included in the medium. Cultivation temperature and pH were 29°C and 3.0, respectively. Three N NaOH and 3N H,SO, were used for pH control. Aeration was 0.5 v/v/m and agitation 500 rpm. Substrate was measured as reducing sugars using the dinitrosalicylic acid (DNS) method. MODEL STRUCTURE System Model
Modeling of growth and substrate consumption of a microorganism and the estimation of model parameters are in many cases problems of larger systems, which may be symbiotic growth of several microorganisms, technically complex processes, for instance, tube processes, or product formation processes. Generally used growth and substrate consumption models in batch cultivation are as follows: dc - _ - p C - kDC dt ds dt
- --
- _1
Y pc
where c is the biomass concentration (or dry weight); s is the substrate concentration; p is the specific growth rate; Y is the yield coefficient, and kD is the death rate Coefficient. The yield and death rate coefficients are often assumed to be constants, although they may in reality change during cultivation. The specific growth rate is not constant but may depend on dry weight and substrate concentration and, of course, on pH and temperature. In this study pH and temperature are assumed to be constant and consequently their effect is not directly taken into account. In some cases maintenance effect has been included in the substrate model: ds
- _-
dt
1
- - p c - mc Y
The maintenance coefficient m can be included into the death rate coefficient by substitution of p' = p
+ Ym
(3)
I834
NIHTILA AND VIRKKUNEN
into eqs. ( I ) and ( 2 ’ ) . The model equations obtained are as follows:
ds = -
-
dt
1 Y
The modified death rate coefficient is now k’,
=
k,
+ Yrn
(4)
Output Model In addition to the system model the output model, which describes the measurements, must be designed. In the simple case in which all growth components are soluble, the dry weight of the biomass can be measured directly. In the process to be considered the substrate is measured as reducing sugars. The measurement method in these experiments also gave nonzero values in the declining phase of the batch cultivation. Consequently it can be assumed that the substrate measurement gives values that deviate from true values, i.e., the measurement is biased. The bias is assumed to be constant, s,, although it probably depends on several factors (e.g., on the initial concentration of the substrate). The measurement equations are as follows: c,
s,
+ w, = s + s, + w,
=c
c, and s, are the measured values of dry weight and substrate concentration. w, and w, describe measurement errors. The model consisting of eqs. ( I ) , ( 2 ) , and ( 5 ) includes parameters k,, Y , s,, and the specific growth rate parameters. If the initial concentrations, c, and so, of the batch cultivation are not sufficiently accurately known, but are measured in the same way as the concentrations during cultivation, they must also be estimated.
Performance Index The mathematical estimation of model parameters is based on minimization of some quantity that can be calculated and that is a function of the parameters to be estimated. If the model under consideration is linear, the estimation is generally an easy task. There exists, however, no general theory for nonlinear parameter estimation.
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1835
Estimation of reaction kinetic constants, including biochemical kinetics, has frequently been performed by modification of the dynamic differential equation models resulting in static equations which are linear with respect to the parameters. After modification linear regression analysis is then applied. A difficulty in the application of this method may, however, be that the disturbances of the variables should be independent of each other. If not, this in turn distorts the results of response fitting. In the case of complex models it is not in general possible to perform this modification. In addition it must be noted that modifications of the models also change the criterion of response fitting if the fitting of modified variables is used. For instance, when estimating parameters of the Monod model in the case kD = 0, the Lineweaver-Burk plot has been used almost without exception. In this case the criterion of response fitting is the minimization of the sum of squares
hiare values of the specific growth rate at the moments ti calculated in some way using dry weight measurements. pc(kl,h,ti) are values of specific growth rate calculated using the model
where smi are measured values of substrate concentration. If the growth curve (gc) is to be fitted with the model, the performance index to be minimized should be N
Jgc(k194) =
C [cc(k1,k,,ti)- cmiI2
(8)
i=l
cc(kl,&,ti)are values of the solution of the following differential equation system for fixed Y, c,, and so: dc - klsc dt h + s dsc dt
1 k,sc Yk,+s
c(0) = c,
(9)
s(0) = so
at the measurement times ti. cmiare measured dry weight values. The optimal parameter values are clearly different for these two fittings because they depend on the form of the performance index. If several variables are measured from the process, for instance,
NIHTILA AND VIRKKUNEN
I836
dry weight and substrate concentration, the suitably weighted sum of fitting errors should be minimized (multiresponse least-squares fitting). DETAILED MODEL SPECIFICATIONS
In this study four different models of growth and substrate consumption are considered. Three of them are specifications of model equations [eqs. ( I ) and (2)] and the fourth model is based on simple enzyme kinetics. 1 ) Monod model:
c. = - k1SC -
kDc
k+s
S=---
1 k;sc Yk,+s
c(0) = c,
s(0) = so
2 ) Contois model:
s=
1 k,sc Yk&+s
s(0) = so
3 ) Linear specific growth rate model (LSGR): d = k,sc
s
=
-
kDc
-(l/Y)klsc
~ ( 0=) C,
(15)
s(0)
(16)
= so
4 ) Enzyme kinetic model of growth (EKG): c(0) = c,
(17)
s(0) = so
(18)
~ ~ (=00 )
(19)
d = YK,c, - kDc S = -K,s(c - c,) C,
=
+ K2c2
K,s(c - c,) - (K2 + K3 + k,)C,
The measurement equations are the same for all the models, i.e., eq. (5). The quadratic performance index to be minimized is N
J(V) =
1 { q [ c ( t j )-
cmi12
+ %[s(ti) - (smj
-
~BII~I
(20)
i=l
where c,,s, are the concentrations at t = 0; cmj,smi are the measured concentrations at t = ti; c ( t , ) ,s ( t i ) are the concentrations calculated using the model; s, is the bias existing in substrate
GROWTH AND SUBSTRATE CONSUMF'TION MODELS
1837
measurements; a,,% are the given weighting coefficients; cz is the fictive cell-substrate complex, and v is the extended parameter vector. The problem is to find the optimal value $ for the parameter vector v that minimizes the performance index J(v). The parameter vector v for the different models is as follows: Monod vM ,
= (CO? SO> 4
9 k Z 7
kD,
y7
SB)
Contois: VC = (cO,
SO?
'$9
4, k D ,
y ?S B )
LSGR: VL = (CO,
so, k ?k D , y, s B )
EKG: VE = ( C o ,
so, K1, KZ,K3r
kD,
Y , ss)
Monod" and Contois" models are well known. The linear specific growth rate model is obtained from the Monod model if during cultivation the substrate concentration is at all times much less than the saturation constant (s 4 k).The kinetic model is, in fact, an extension of the Monod model and reduces into it when assuming pseudo-steady-state for the cell-substrate complex (see Appendix A). The evident advantage of the kinetic model from the point of view of estimation is that the model equations are linear with respect to the parameters, although the model includes one parameter and one differential equation more than the Monod model. ESTIMATION METHOD
There exist in the literature a large number of different iterative methods for nonlinear dynamic parameter estimations. The applicability and usefulness of these methods depend mainly on the problem to be considered. The estimation problem posed in this paper has been solved using the Gauss-Newton-type method modified by Marquardt. The method was initially developed for static problems. It has also been very successfully applied by BardI3 to dynamic estimation problems. The minimization problem is solved iteratively by linearization of the output function with respect to the initial state and the parameter vector, and by solving approximately the obtained
NIHTILA AND VIRKKUNEN
I838
quadratic minimization problem. The solution of the quadratic problem is the new estimate for the initial state and for the parameter vector. The linearization is then performed again in the neighborhood of this new estimate. The quadratic problem is solved approximately in such a way that the new estimate gives a smaller value for the performance index than the old estimate gave (Appendix B). This estimation method is purely deterministic and it does not take into account any statistical properties of the measurement or modeling errors. The method can, however, be considered in stochastic sense asssuming that the process model is correct and that measurement errors, i.e., measurement noise w ( t ) (see Appendix B), can be modeled as a zero mean-valued normally distributed random process with unknown covariance. The linearized covariance matrix of the estimate according to H i m m e l b l a ~ , 'p. ~ 197, is then N
-1
R, = (2 HjTWiHi) S,2 i=l
where s,2
= J(+>/f
(22)
is the variance of the residual error of estimation. J(+) is the minimal value of the performance index, and
f = mN
- (n
+ n,)
(23)
is the degree of freedom of estimation. Other notations have been explained in Appendix B. The individual confidence limits of the estimates at probability level ( I - a ) are14*15 ui = Gi ? tCfi 1 -
a/2).(R,.ij)1'2
i
=
1, 2 ,
. . . , n + n,
(24)
where tCf, 1 - 4 2 ) is the value of the t distribution for fdegrees of freedom and for probability ( 1 - a/2).The values of (Z&)'" are standard deviations of the individual parameters u i , i.e., square roots of the diagonal elements of the covariance matrix [eq. (21)l. The (n + n,)-dimensional joint confidence region at probability level I - a is inside the hyperellipsoid described by the equation (v - f)'Rv-'(v
+
- +) = (n + n,) F(n
+ n,,f,
1 - a)
(25)
where F(n n,, f, 1 - a)is the value of F distribution for (n + np, f ) degrees of freedom and for probability 1 - a.
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1839
EXPERIMENTAL RESULTS
The experimental measurement data used for process identification are shown in Table I. The weighting coefficients of the performance index [eq. (20)] were a, = a, = 1. The estimations for the three models (Monod, Contois, and LSGR) were also successfully performed mathematically because the optima were obtained. The parameter values converged into the optimal ones in four to ten steps. Difficulties were encountered with the enzyme kinetic model because it included the greatest number of parameters and initial values to be estimated. The main difficulty was very slow convergence of some model parameters. l6 In addition, only suitable values, not optimal, were found for this model. The reasons for these difficulties were the large numbers of parameters and especially the stiffness of the differential equation model (stiffness means that the model also includes, besides slow modes, very fast ones). Consequently the integration method used, a fourth order Runge-Kutta-Gill algorithm, did not work satisfactorily in this case. TABLE I Measured Dry Weight and Substrate Concentration as a Function of Time for a Batch Cultivation of 7'. viride Grown on Glucose as the Main Carbon Source Time (daydhr)
0.0 0.09 0.12 0.15 0.17 0.20 1.o
1.08 1.14 2.0 2.09,a 2.14 2.17 3.17 4.18, 5.0 a
Dry weight (dliter) 0.4
0.7 0.9 1.1 1.5 1.8 2.5 3.9 7.7 8.2 8.7 8.5 8.5 7.9 6.4 6.0
Substrate concentration (dliter) 24.5 27.5 25.0 23.0 23.5 22.0 22.5 19.0
16.0 7.5 3.0 2.0 2.0 2.5 1.6 1.2
Two days and nine and one-half hours.
1840
NIHTILA AND VIRKKUNEN
The least-squares estimates and individual 95% confidence limits (about twice the standard deviation) are given in Table 11. The joint confidence regions were too large to yield any significance for reliability tests. The large limits are because too many parameters had to be estimated from the available measurements. Figure I presents the responses of the models calculated using least-squares estimates compared with measured dry weight and substrate concentrations. Conclusions
The identification results for the Monod and Contois models are almost equal for the respective parameter values and for the performance indices. Consequently without any statistical tests it can be said that both models describe equally well the experimental cultivation. The generally used Monod model will, however, be used in the simulation analysis of the presented estimation method. It should be noted that although the linear specific growth rate model gave clearly poorer fitting in a numerical sense, and although for the enzyme kinetic model the authors were unable to find the optimal values for the estimates, the responses of all four models were on visual inspection very similar. The individual 95% confidence limits are rather large; for all the
TABLE I1 Least-Squares Estimates and Individual 95% Confidence Limits of Initial Values and Parameters, and the Performance Index for the Four Batch Cultivation Models Linear Parameters co (&liter)
so (dliter)
kl (day-')
k2 (&liter) kl (g.day-') k, Y kD (day-') K , (g.day-3 K , (day-? K , (day-') sE (giliter) J((g/lWz) a
Monod
Contois
0.468 2 0.49 24.3' ? 1.7 2.51 If: 2.5 9.34 2 19.2
0.467 ? 0.49 24.3' ? 1.7 1.84 2 0.89
0.468 ? 0.11 0.203 2 0.12
-
1.81 ? 1.2 22.67
(LSGR) 0.135 ? 0.14 24.5' 2 2.0
-
0.668 ? 0.99
0.123 2 0.03
0.468 ? 0.21 0.202 If: 0.12
0.524 2 0.12 0.243 2 0.14
1.83 2 1.2 22.73
Measurement bias sE not included in the value.
-
1.17 ? 1.4 27.87
Kinetic (EKG) 0.300 2 0.54 24.7a 2 4.3
-
-
0.518 2 0.81 0.242 ? 0.67 0.376 ? 2.64 1.51 ? 86.8 12.1 2 5.9 1.52 k 3.8 28.03
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1841
Monod --- Contois -.-._.
LSGR _-_____ Kinetic
S
oncen tra tion
m
complex
10
20
1.0 /
/
1
2
3
5
L
Time/days
c 10
-
Dry m i g h t (9/1)
2
3
5
L Time/doys
Fig. I . Experimental measurements of T . viride cultivated on glucose as main carbon source, compared with fitted least-squares responses of the Monod, Contois, LSGR, and EKG models. Single dotted curve in the upper figure is complex concentration of the kinetic model.
1842
NIHTILA AND VIRKKUNEN
models the estimated initial concentrations of dry weight d o not deviate significantly from zero. It should be, however, evident that the calculated confidence limits are only crude approximations because of the nonlinearities of the models. The probability theory also assumes that the measurement noise w(t) has the required normalcy properties, which is not always the case. By suitable experimental design, l7 i.e., by selection of suitable measurement times and measurement density, confidence regions can be made smaller. The sensitivity functions (see, e.g., Sage,18 p. 323) of the measured variables yield information for this selection. The linear estimation theory shows that in the case of normally distributed measurement noise the estimate is unbiased, i.e., the expected value of the estimate is equal to the real parameter value. This means, in practice, that when the number of measured points is increased, the estimate becomes closer and closer to the true parameter value. For nonlinear cases the estimate is, in general, biased and only in simple special cases can the value of this bias be calculated ( S c h ~ e p p e ,p.~ 335). The biasing causes difficulties because the biased estimate may be outside the individual confidence region that was calculated with the aid of the linearized covariance matrix. The bias error and also the result of estimation are affected by the level of measurement noise (i.e., its variance), a s well as by the number and times of measurements.
SIMULATION OF THE ESTIMATION METHOD When interpreting the results of nonlinear parameter estimation care is needed, because, in general, it cannot be certain whether an obtained result is a global or local optimum and because a globally optimal result may be biased by an unknown constant whose value cannot be calculated. Also, the reliability limits obtained by linearization may be incorrect. Because of these considerations, simulation is a powerful tool for determining to what extent the results can be correct. The optimality and possible biasing can easily be revealed because the correct parameter values are known. The sample standard deviations and reliability limits that correspond to real values can also be calculated by simulating several times using independent noise every time. The Monod model is used in the simulation study because it is the best known and most frequently used model. The following parameter values
GROWTH AND SUBSTRATE CONSUMFITON MODELS
I843
were used for simulations: c,
= 0.5
S,
= 20.0 g/liter
kD = 0.2 day-’
g/liter
Y = 0.5 s, = 2.0g/liter
k, = 2.5 day-’
k, = IO.Og/liter These values are close to those given by estimation. All the initial values and parameters were estimated in this study. Modeling and measurement errors were taken into account by generating normally distributed pseudorandom noise additively into measurement as follows c,i
= C(ti)
+ w,(r,)
smi = s ( t i ) + S,
+ w&)
c ( f i )and ~ ( t are , ) the calculated dry weight and substrate concentration, respectively. w,( r ) and w,(f) are independent measurement errors. The cultivation time of the simulated batch run was five days. A total of 26 equally distributed values were “measured” from both variables. The weighting coefficients of the performance index were a, = a, = 1. Standard deviations (u)of the measurement noise vaned from 0.0 to 0.6 g/liter with steps of 0.1 g/liter, and were same for both measurements. Each standard deviation was used five times. Figure 2 presents as an example the optimal estimates of the yield Y and the maximum specific growth rate k, and the theoretical standard deviations obtained by linearization. The sample means and standard deviations were also calculated using estimation results. It is easily seen that standard deviations obtained by linearization are greater than standard deviations calculated from estimation results, i.e., from samples. The slight increasing effect of biasing can also be seen in the values of yield Y when noise level increased. The same results were also obtained starting from different initial values for iteration. DISCUSSION The presented nonlinear parameter estimation method was based on clear mathematical minimization of the cost functional, which depends on the pure unmodified measurement data and on the
I
04
7
a2
1
0
d (g/l) 0.6
14 0
0.2
0.2
d (g/L)
0.6
Fig. 2. Estimated yield Y and maximum specific growth rate k , for simulated batch runs using different levels of measurement noise with t a n d a r d deviation w. Bars represent standard deviations calculated by linearization. Dotted lines represent sample means and sample t andard deviations. ( 0 )Estimated values, (0) sample means of five estimates.
0.2
GROWTH AND SUBSTRATE CONSUMPTION MODELS
184.5
extended parameter vector to be estimated. However, nodineanties of the models made interpretation of the results impossible on the basis of any general theory of estimation. Such a theory has not yet been developed. Simulation of models and of estimation is a powerful tool when considering the possibilities of estimation based on experimental measurements and when interpreting confidence limits and reliability of estimates. The question of parameter identifiability and the effect of measurement noise on bias of estimate have been considered by simulation. Estimations using measurements from a real fermentation were successfully performed mathematically for the first three models. The convergence was obtained in four to ten steps depending on the model and on the starting values for estimation. Because of insufficient measurements, the linearized confidence limits are, however, approximate ones. Feasible parameter values were also found in the case of the enzyme kinetic model of growth, although a mathematical optimum was not obtained. The main reasons for this were the large number of parameters to be estimated compared to the number of measurements, the missing measurements of the third state variable of the model, i.e., cell-substrate complex, and especially the stiffness property of the differential equation model, which made the use of the fourth order Runge-Kutta-Gill integration method impractical. After fixing some of the model parameters, convergence was also obtained mathematically, i.e., the optimal estimate. After visual inspection the responses of all the models are, however, indistinguishable in the goodness of fit. The authors have found only few articles in the literature concerning the growth and substrate consumption of T. viride on glucose. Brown and Halstedlg considered the effect of pH on the maximum specific growth rate. The growth was exponential according to their conclusions. This means that the saturation constant of the Monod model is small compared to the values of substrate concentration, i.e., about I .O g/liter. Nagai et al.'O used an unidentified Trichoderma species growing on glucose and on cellobiose. The estimated parameters of the Monod model are compared with values found in the literature (Table 111). The values of the yield Y and the maximum specific growth rate k, are in the same range. The estimated saturation constant k, is much greater than values found in the literature. (An interesting detail is that two papers present a value of k, equal to each other.) Calculated confidence limits of the saturation constant are larger than the
1846
NIHTILA AND VIRKKUNEN
TABLE 111 Growth and Substrate Consumption Parameters of the Monod Model and Environmental Conditions from Literature Compared with Values of the Present Authors
Parameters k , (day-') k Z (&liter) Y T("C)
PH
Brown and Nagai et aLZ0 HalstedlS 3.12 0.083b 0.656d 30 6.3 + 3.5
Mitra and Wilke" 7.06 0.083 0.43 30 -3.2
-1.85 1' 0.40 28 3.0
Mandels et aLz2
Present authors
4.8" 0.3-0.4 -
2.5 1 9.3 0.468 29 3.0
Approximate values for continuous culture. 0.46mM glucose. Assumed by the present authors. * 118 g celVM glucose. a
confidence limits of the other parameters. This means that it is not possible to estimate the saturation constant accurately using only 16 measurements of both variables. In order to improve the estimates' accuracy the measurement noise should be decreased and the number of measurements increased. Using experimental design in selecting suitable measurement times the accuracy can also be improved without increasing the number of measured points. If the accuracy of one parameter estimate is to be improved, the measurements can be weighted using the inverse values of the sensitivity functions of biomass and substrate concentration with respect to the parameter to be considered. This, however, changes the original problem statement and possibly increases the confidence regions of the other parameters. APPENDIX A If it is assumed that substrate and biomass form a reversible complex that grows, the following enzyme kinetic scheme is obtained for the growth and substrate consumption process: KI
c, + s *
K3
C, 3 (1
+ u)C,
KZ
3-
0
kD
3-
kD
0
where C, is pure biomass; S is substrate; C, is the cell-substrate
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1847
complex; Y is the yield coefficient; K , , K,, and K3 are the reaction rate constants, and kD is the death rate coefficient. The result of dry weight measurements is the sum of pure biomass and complex. Consequently for the respective concentrations the following measurement equation is satisfied: c
=
c,
+ cp
(26)
From the reaction scheme and taking into account the measurement equation [eq. (26)] the differential equation model for the concentrations is obtained C
=
YK3c2- kDc
S = -K,s(c C,
=
- c,)
+ K,c,
K,s(c - c,) - ( K ,
+ K3 + k , ) ~ ,
c(0) = c,
(27)
s(0)
= so
(28)
~ ~ (=0 0)
(29)
Pseudo-steady-state assumption for the cell-substrate complex leads to the equation c,
=
K,sc/(K,
+ K3 + kD + K , s )
(30)
Substitution of this into eqs. (27) and (28) results in the Monod model with parameters (c.f., eqs. ( 1 1) and (12)):
k,
= YK,
k 2
=
(K2
(3 1)
+ K3 + k D ) / K 1
(32)
APPENDIX B This method presents, in a general form, the parameter estimation of the process model dx - = f(x, p, t ) dt
(33)
y o ) = g(x(t),p, 0 + w(d
(34)
by minimizing the performance index N
J =
C [y(ti) - g(x(ti)$p , ti)]'
wi[y(ti) - g(x(ti),p, ti)]
(35)
i=l
where x is the state variable of the process, dim(x) = n ; y is the output variable, dim(y) = m; p is the parameter vector, dim(p) = np;5 is the initial value of state x at f = 0; W iare positive definite
NIHTILA AND VIRKKUNEN
I848
weighting matrices, dim(Wi) = m X m ; N is the number of measurement times, and [ IT denotes the transpose of a vector (or matrix). The solution of the differential equation [eq. (33)] depends only on the initial state and parameter vector. Consequently, it can be expressed in symbolic form x(f) = F(t, 5, p)
(36) When this is substituted into the performance index [eq. (35)], it becomes solely the function of the initial state and parameter vector
(37)
J = J(&, p)
When linearizing the output function with regard to the initial guess,
I::[ the following equation is obtained: g(x, p, t ) = go
ago aFO + --A& agT axT
+
Derivative notations are according to C ~ a k i , i.e., * ~ elements of the matrices are of the form
The superscript zero indicates that the functions are evaluated using the initial guess. Sensitivity functions are the solutions of the following linear matrix differential equations:
d-A- afo --A dt axT
A(0)
= Z(identity
(39) (40)
M(0) = 0
where dim(A) = n x n dim(M) = n
matrix)
X
n,
1849
GROWTH AND SUBSTRATE CONSUMPTION MODELS
The output function is now expressed at all time instants form g(x(t;), p, t i ) = go +
ti
in the (43)
HiAV
where Av
=
[i:]
dim(Av)
=
n
+ np
(44)
.
where t = ti, i = 1, 2, . - . , M. When substituting eq. (43) into the cost and minimizing it with respect to Av, the linear equation for Av is obtained: AAv
=
b
(46)
A is a matrix with dimensions (n + n,) X ( n + n,) and the dimension of the gradient vector b is n + np. The change Av of the estimate is calculated using Marquardt’s method, by solving the equation
(A
+ AI)Av = b
(47)
for the value X that agrees with the equation
fito+ &(A),
PO
+ Ap(A>l < J(e0,PO)
(48)
Proceeding as described and using a suitable stopping rule the optimal cost minimizing estimate
is obtained. The authors thank the personnel of the Biotechnical Laboratory, Technical Research Center of Finland, for carrying out the experimental cultivation according to the authors’ needs. This study is a part of the research project “Development, Control, and Optimization of Biotechnical Processes,” supported by the Academy of Finland and undertaken as a joint research of the Control Engineering Laboratory and the Biotechnical Laboratory.
References I . R. P. Jefferis, 111, Biochem., 10, 15 (1975). 2. R. P. Jefferis, in Fifth International Fermentation Symposium. Abstracts of Papers, H. Dellweg, Ed., Westkreuz-Druckerei und Verlag, Berlin, 1976, p. 34.
1850
NIHTILA AND VIRKKUNEN
3. D. D. Y . Ruy and A. E. Humphrey, J. Appl. Chem. Biotechnol., 23, 283 (1973). 4. W. B. Arminger, D. W. Zabriskie, and A. E. Humphrey, in Fifth International Fermentation Symposium, Abstracts of Papers, H. Dellweg, Ed., WestkreuzDruckerei und Verlag, Berlin, 1976, p. 33. 5. H. Metz, F. Wenzel, and E. Merck, in Fifth International Fermentation Symposium, Abstracrs of Papers, H . Dellweg, Ed., Westkreuz-Druckerei und Verlag, Berlin, 1976, p. 35. 6. A. Geranton and M. Colas, Automatisme, 20,452 (1975). 7. A. Meskanen, R. Lundell, and P. Laiho, Process Biochem.. 11, 31 (1976). 8. R. Spruytenburg, N. D. P. Dang, I. J . Dunn, J. R. Mor, A. Einsele, A. Fiechter, and J. R. Bourne, Chem. Eng., 310, 447 (1976). 9. F. C. Schweppe, Uncertain Dynamic Systems, Prentice-Hall, Englewood Cliffs, N. J . , 1973. 10. J. Monod, Ann. Rev. Microbiol.. 3, 371 (1949). 1 1 . J. Lescourret, Ph.D. thesis, Universite Paul Sabatier d e Toulouse, 1974. 12. D. W. Marquardt, J. Soc. Indust. Appl. Math., 11,431 (1%3). 13. Y. Bard, SIAM J . Numer. Anal., 7, 157 (1970). 14. D. M. Himmelblau, Process Analysis by Statistical Methods, Wiley, New York, 1970. 15. G. Emig and L . H. Hosten, Chem. Eng. Sci., 29,475 (1974). 16. M. Nihtila and J. Virkkunen, in Proceedings of 5th IFACIIFIP International Conference on Digital Computer Applications to Process Control, Hague, The Netherlands, 1977. 17. D. Johnsson and P. Berthouex, Biotechnol. Bioeng., 17, 557 (1975). 18. A. P. Sage, Optimum Systems Control, Prentice-Hall, Englewood Cliffs, N. J., 1968, p. 323. 19. D. Brown and D. Halsted, Biotechnol. Bioeng., 17, I 199 (1975). 20. S. Nagai, M. Onodera, and S. Aiba, Ear. J . Appl. Microbiol., 1, 9 (1976). 21. G. Mitra and C. Wilke, Biotechnol. Bioeng., 27, I (1975). 22. M. Mandels, D. Sternberg, and R. Andreotti, in Symposium on Enzymatic Hydrolysis of Cellulose, M. Bailey, T-M. Enari, and M. Linko, Eds., The Finnish National Fund for Research and Development (SITRA), Helsinki, 1975. 23. F. Csaki, Periodica Polytechnica El. (Technical University. Budapest), 15 (1971).
Accepted for Publication May 28, 1977
Practical Identifiability of Growth and Substrate Consumption Models MARKKU NIHTILA and JOUKO VIRKKUNEN, Helsinki University of Technology, Electrical Engineering Department, Control Engineering Laboratory, SF-02150 Espoo 15, Finland Summary The estimation of parameters in several dynamic models, which describe growth and substrate consumption, has been carried out using a modified Gauss-Newtontype method. The four models considered are Monod, Contois, linear specific growth rate, and an enzyme kinetic model. The initial values of the differential equations are included in the parameter vector which will be estimated. The efficiency of the method and the confidence limits of the parameters were studied using simulated measurement noise. The experimental results describe Tn'chodrrma viride growing on glucose as the main carbon source.
INTRODUCTION
The use of mathematical models has become increasingly popular for the description of microbiological processes. The extensive computational effort necessary in using models has become less significant owing to the advent of powerful computers and automatic data logging facilities. In recent years several computerized pilot plants have been developed'-' and there is an urgent need to develop models and computational algorithms for the utilization of the great amounts of information which can be collected during experimental runs. In developing mathematical models the first and basic aim is to give a concise and unambiguous description of the biological mechanism of the process. The different mechanisms and the pertinent models can be distinguished by comparing measured results with the models' predictions. However, the model is most often a group of first order nonlinear differential equations with unknown parameters. The parameters must be estimated by using measured results. A satisfying fit of the model response and experimental results can often be achieved, especially if the number of free parameters is high and the number of measured 1831
@ 1977 by John Wiley & Sons, Inc.
1832
NIHTILA AND VIRKKUNEN
points is comparatively low, as is often the case in microbiological experiments. This good fit can be highly misleading as the confidence limits of the estimated parameters can be so great that the actual numerical values of the parameters have no relevance at all. Unfortunately there are no general methods for the computation of confidence limits for parameters in nonlinear estimation problems. Moreover, the parameters are normally interdependent and the estimated values can be biased. This bias is practically impossible to find by mathematical methods. Simple graphical means such as the Lineweaver-Burk plot can only be used in simple cases. The parameters in more complicated models must be estimated by using rather sophisticated mathematical methods which are based on the minimization of a criterion of the goodness of fit. In nonlinear problems, however, there can be several local minima and the algorithm can give erroneous results by leading to a secondary minimum of the criteria. In the opinion of the present authors, as stated also, e.g., by S c h ~ e p p e ,the ~ most practical and feasible way of testing the reliability of the parameters in nonlinear models is to simulate the measurements using either simulated values for measurements or using real measurements and adding to them simulated noise. A study of this type is reported here. The first section briefly describes the materials and methods used. Model structure is considered in the second section by dividing the problem into three parts: the system model, the output model, and the performance index. Mathematical model specifications which are used in parameter estimation are presented in detail in the third section. The estimation method is then qualitatively explained. (Mathematical details of the method are given in Appendix B.) Results of estimation, i.e., parameter values, confidence limits, and residual sums with relevant conclusions, are presented in the next sections. Finally, a simulation study of the method, performed using Monod’s model with feasible parameter values, is described. MATERIALS AND METHODS
The microorganism used was Trichoderma viride strain QM 9414, obtained from the U.S. Army Natick Laboratories. The inoculum (100 ml) was four days old, grown on 1% Solca Floc medium. The fermentor was a New Brunswick Magnaferm, working volume 10 liter. The medium contained in dliter: (NHk),SO,, 2.8; KH2P04, 15.0; MgSO,, 0.6; CaC12,0.6; urea, 0.6; peptone, 1.5.
GROWTH AND SUBSTRATE CONSUMFTION MODELS
I833
The carbon source was a glucose syrup (n% dry weight) containing 20% glucose. Tween 80 (0.2% v/v) and trace elements (Fe", Mn++, Zn"", C o " ) were also included in the medium. Cultivation temperature and pH were 29°C and 3.0, respectively. Three N NaOH and 3N H,SO, were used for pH control. Aeration was 0.5 v/v/m and agitation 500 rpm. Substrate was measured as reducing sugars using the dinitrosalicylic acid (DNS) method. MODEL STRUCTURE System Model
Modeling of growth and substrate consumption of a microorganism and the estimation of model parameters are in many cases problems of larger systems, which may be symbiotic growth of several microorganisms, technically complex processes, for instance, tube processes, or product formation processes. Generally used growth and substrate consumption models in batch cultivation are as follows: dc - _ - p C - kDC dt ds dt
- --
- _1
Y pc
where c is the biomass concentration (or dry weight); s is the substrate concentration; p is the specific growth rate; Y is the yield coefficient, and kD is the death rate Coefficient. The yield and death rate coefficients are often assumed to be constants, although they may in reality change during cultivation. The specific growth rate is not constant but may depend on dry weight and substrate concentration and, of course, on pH and temperature. In this study pH and temperature are assumed to be constant and consequently their effect is not directly taken into account. In some cases maintenance effect has been included in the substrate model: ds
- _-
dt
1
- - p c - mc Y
The maintenance coefficient m can be included into the death rate coefficient by substitution of p' = p
+ Ym
(3)
I834
NIHTILA AND VIRKKUNEN
into eqs. ( I ) and ( 2 ’ ) . The model equations obtained are as follows:
ds = -
-
dt
1 Y
The modified death rate coefficient is now k’,
=
k,
+ Yrn
(4)
Output Model In addition to the system model the output model, which describes the measurements, must be designed. In the simple case in which all growth components are soluble, the dry weight of the biomass can be measured directly. In the process to be considered the substrate is measured as reducing sugars. The measurement method in these experiments also gave nonzero values in the declining phase of the batch cultivation. Consequently it can be assumed that the substrate measurement gives values that deviate from true values, i.e., the measurement is biased. The bias is assumed to be constant, s,, although it probably depends on several factors (e.g., on the initial concentration of the substrate). The measurement equations are as follows: c,
s,
+ w, = s + s, + w,
=c
c, and s, are the measured values of dry weight and substrate concentration. w, and w, describe measurement errors. The model consisting of eqs. ( I ) , ( 2 ) , and ( 5 ) includes parameters k,, Y , s,, and the specific growth rate parameters. If the initial concentrations, c, and so, of the batch cultivation are not sufficiently accurately known, but are measured in the same way as the concentrations during cultivation, they must also be estimated.
Performance Index The mathematical estimation of model parameters is based on minimization of some quantity that can be calculated and that is a function of the parameters to be estimated. If the model under consideration is linear, the estimation is generally an easy task. There exists, however, no general theory for nonlinear parameter estimation.
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1835
Estimation of reaction kinetic constants, including biochemical kinetics, has frequently been performed by modification of the dynamic differential equation models resulting in static equations which are linear with respect to the parameters. After modification linear regression analysis is then applied. A difficulty in the application of this method may, however, be that the disturbances of the variables should be independent of each other. If not, this in turn distorts the results of response fitting. In the case of complex models it is not in general possible to perform this modification. In addition it must be noted that modifications of the models also change the criterion of response fitting if the fitting of modified variables is used. For instance, when estimating parameters of the Monod model in the case kD = 0, the Lineweaver-Burk plot has been used almost without exception. In this case the criterion of response fitting is the minimization of the sum of squares
hiare values of the specific growth rate at the moments ti calculated in some way using dry weight measurements. pc(kl,h,ti) are values of specific growth rate calculated using the model
where smi are measured values of substrate concentration. If the growth curve (gc) is to be fitted with the model, the performance index to be minimized should be N
Jgc(k194) =
C [cc(k1,k,,ti)- cmiI2
(8)
i=l
cc(kl,&,ti)are values of the solution of the following differential equation system for fixed Y, c,, and so: dc - klsc dt h + s dsc dt
1 k,sc Yk,+s
c(0) = c,
(9)
s(0) = so
at the measurement times ti. cmiare measured dry weight values. The optimal parameter values are clearly different for these two fittings because they depend on the form of the performance index. If several variables are measured from the process, for instance,
NIHTILA AND VIRKKUNEN
I836
dry weight and substrate concentration, the suitably weighted sum of fitting errors should be minimized (multiresponse least-squares fitting). DETAILED MODEL SPECIFICATIONS
In this study four different models of growth and substrate consumption are considered. Three of them are specifications of model equations [eqs. ( I ) and (2)] and the fourth model is based on simple enzyme kinetics. 1 ) Monod model:
c. = - k1SC -
kDc
k+s
S=---
1 k;sc Yk,+s
c(0) = c,
s(0) = so
2 ) Contois model:
s=
1 k,sc Yk&+s
s(0) = so
3 ) Linear specific growth rate model (LSGR): d = k,sc
s
=
-
kDc
-(l/Y)klsc
~ ( 0=) C,
(15)
s(0)
(16)
= so
4 ) Enzyme kinetic model of growth (EKG): c(0) = c,
(17)
s(0) = so
(18)
~ ~ (=00 )
(19)
d = YK,c, - kDc S = -K,s(c - c,) C,
=
+ K2c2
K,s(c - c,) - (K2 + K3 + k,)C,
The measurement equations are the same for all the models, i.e., eq. (5). The quadratic performance index to be minimized is N
J(V) =
1 { q [ c ( t j )-
cmi12
+ %[s(ti) - (smj
-
~BII~I
(20)
i=l
where c,,s, are the concentrations at t = 0; cmj,smi are the measured concentrations at t = ti; c ( t , ) ,s ( t i ) are the concentrations calculated using the model; s, is the bias existing in substrate
GROWTH AND SUBSTRATE CONSUMF'TION MODELS
1837
measurements; a,,% are the given weighting coefficients; cz is the fictive cell-substrate complex, and v is the extended parameter vector. The problem is to find the optimal value $ for the parameter vector v that minimizes the performance index J(v). The parameter vector v for the different models is as follows: Monod vM ,
= (CO? SO> 4
9 k Z 7
kD,
y7
SB)
Contois: VC = (cO,
SO?
'$9
4, k D ,
y ?S B )
LSGR: VL = (CO,
so, k ?k D , y, s B )
EKG: VE = ( C o ,
so, K1, KZ,K3r
kD,
Y , ss)
Monod" and Contois" models are well known. The linear specific growth rate model is obtained from the Monod model if during cultivation the substrate concentration is at all times much less than the saturation constant (s 4 k).The kinetic model is, in fact, an extension of the Monod model and reduces into it when assuming pseudo-steady-state for the cell-substrate complex (see Appendix A). The evident advantage of the kinetic model from the point of view of estimation is that the model equations are linear with respect to the parameters, although the model includes one parameter and one differential equation more than the Monod model. ESTIMATION METHOD
There exist in the literature a large number of different iterative methods for nonlinear dynamic parameter estimations. The applicability and usefulness of these methods depend mainly on the problem to be considered. The estimation problem posed in this paper has been solved using the Gauss-Newton-type method modified by Marquardt. The method was initially developed for static problems. It has also been very successfully applied by BardI3 to dynamic estimation problems. The minimization problem is solved iteratively by linearization of the output function with respect to the initial state and the parameter vector, and by solving approximately the obtained
NIHTILA AND VIRKKUNEN
I838
quadratic minimization problem. The solution of the quadratic problem is the new estimate for the initial state and for the parameter vector. The linearization is then performed again in the neighborhood of this new estimate. The quadratic problem is solved approximately in such a way that the new estimate gives a smaller value for the performance index than the old estimate gave (Appendix B). This estimation method is purely deterministic and it does not take into account any statistical properties of the measurement or modeling errors. The method can, however, be considered in stochastic sense asssuming that the process model is correct and that measurement errors, i.e., measurement noise w ( t ) (see Appendix B), can be modeled as a zero mean-valued normally distributed random process with unknown covariance. The linearized covariance matrix of the estimate according to H i m m e l b l a ~ , 'p. ~ 197, is then N
-1
R, = (2 HjTWiHi) S,2 i=l
where s,2
= J(+>/f
(22)
is the variance of the residual error of estimation. J(+) is the minimal value of the performance index, and
f = mN
- (n
+ n,)
(23)
is the degree of freedom of estimation. Other notations have been explained in Appendix B. The individual confidence limits of the estimates at probability level ( I - a ) are14*15 ui = Gi ? tCfi 1 -
a/2).(R,.ij)1'2
i
=
1, 2 ,
. . . , n + n,
(24)
where tCf, 1 - 4 2 ) is the value of the t distribution for fdegrees of freedom and for probability ( 1 - a/2).The values of (Z&)'" are standard deviations of the individual parameters u i , i.e., square roots of the diagonal elements of the covariance matrix [eq. (21)l. The (n + n,)-dimensional joint confidence region at probability level I - a is inside the hyperellipsoid described by the equation (v - f)'Rv-'(v
+
- +) = (n + n,) F(n
+ n,,f,
1 - a)
(25)
where F(n n,, f, 1 - a)is the value of F distribution for (n + np, f ) degrees of freedom and for probability 1 - a.
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1839
EXPERIMENTAL RESULTS
The experimental measurement data used for process identification are shown in Table I. The weighting coefficients of the performance index [eq. (20)] were a, = a, = 1. The estimations for the three models (Monod, Contois, and LSGR) were also successfully performed mathematically because the optima were obtained. The parameter values converged into the optimal ones in four to ten steps. Difficulties were encountered with the enzyme kinetic model because it included the greatest number of parameters and initial values to be estimated. The main difficulty was very slow convergence of some model parameters. l6 In addition, only suitable values, not optimal, were found for this model. The reasons for these difficulties were the large numbers of parameters and especially the stiffness of the differential equation model (stiffness means that the model also includes, besides slow modes, very fast ones). Consequently the integration method used, a fourth order Runge-Kutta-Gill algorithm, did not work satisfactorily in this case. TABLE I Measured Dry Weight and Substrate Concentration as a Function of Time for a Batch Cultivation of 7'. viride Grown on Glucose as the Main Carbon Source Time (daydhr)
0.0 0.09 0.12 0.15 0.17 0.20 1.o
1.08 1.14 2.0 2.09,a 2.14 2.17 3.17 4.18, 5.0 a
Dry weight (dliter) 0.4
0.7 0.9 1.1 1.5 1.8 2.5 3.9 7.7 8.2 8.7 8.5 8.5 7.9 6.4 6.0
Substrate concentration (dliter) 24.5 27.5 25.0 23.0 23.5 22.0 22.5 19.0
16.0 7.5 3.0 2.0 2.0 2.5 1.6 1.2
Two days and nine and one-half hours.
1840
NIHTILA AND VIRKKUNEN
The least-squares estimates and individual 95% confidence limits (about twice the standard deviation) are given in Table 11. The joint confidence regions were too large to yield any significance for reliability tests. The large limits are because too many parameters had to be estimated from the available measurements. Figure I presents the responses of the models calculated using least-squares estimates compared with measured dry weight and substrate concentrations. Conclusions
The identification results for the Monod and Contois models are almost equal for the respective parameter values and for the performance indices. Consequently without any statistical tests it can be said that both models describe equally well the experimental cultivation. The generally used Monod model will, however, be used in the simulation analysis of the presented estimation method. It should be noted that although the linear specific growth rate model gave clearly poorer fitting in a numerical sense, and although for the enzyme kinetic model the authors were unable to find the optimal values for the estimates, the responses of all four models were on visual inspection very similar. The individual 95% confidence limits are rather large; for all the
TABLE I1 Least-Squares Estimates and Individual 95% Confidence Limits of Initial Values and Parameters, and the Performance Index for the Four Batch Cultivation Models Linear Parameters co (&liter)
so (dliter)
kl (day-')
k2 (&liter) kl (g.day-') k, Y kD (day-') K , (g.day-3 K , (day-? K , (day-') sE (giliter) J((g/lWz) a
Monod
Contois
0.468 2 0.49 24.3' ? 1.7 2.51 If: 2.5 9.34 2 19.2
0.467 ? 0.49 24.3' ? 1.7 1.84 2 0.89
0.468 ? 0.11 0.203 2 0.12
-
1.81 ? 1.2 22.67
(LSGR) 0.135 ? 0.14 24.5' 2 2.0
-
0.668 ? 0.99
0.123 2 0.03
0.468 ? 0.21 0.202 If: 0.12
0.524 2 0.12 0.243 2 0.14
1.83 2 1.2 22.73
Measurement bias sE not included in the value.
-
1.17 ? 1.4 27.87
Kinetic (EKG) 0.300 2 0.54 24.7a 2 4.3
-
-
0.518 2 0.81 0.242 ? 0.67 0.376 ? 2.64 1.51 ? 86.8 12.1 2 5.9 1.52 k 3.8 28.03
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1841
Monod --- Contois -.-._.
LSGR _-_____ Kinetic
S
oncen tra tion
m
complex
10
20
1.0 /
/
1
2
3
5
L
Time/days
c 10
-
Dry m i g h t (9/1)
2
3
5
L Time/doys
Fig. I . Experimental measurements of T . viride cultivated on glucose as main carbon source, compared with fitted least-squares responses of the Monod, Contois, LSGR, and EKG models. Single dotted curve in the upper figure is complex concentration of the kinetic model.
1842
NIHTILA AND VIRKKUNEN
models the estimated initial concentrations of dry weight d o not deviate significantly from zero. It should be, however, evident that the calculated confidence limits are only crude approximations because of the nonlinearities of the models. The probability theory also assumes that the measurement noise w(t) has the required normalcy properties, which is not always the case. By suitable experimental design, l7 i.e., by selection of suitable measurement times and measurement density, confidence regions can be made smaller. The sensitivity functions (see, e.g., Sage,18 p. 323) of the measured variables yield information for this selection. The linear estimation theory shows that in the case of normally distributed measurement noise the estimate is unbiased, i.e., the expected value of the estimate is equal to the real parameter value. This means, in practice, that when the number of measured points is increased, the estimate becomes closer and closer to the true parameter value. For nonlinear cases the estimate is, in general, biased and only in simple special cases can the value of this bias be calculated ( S c h ~ e p p e ,p.~ 335). The biasing causes difficulties because the biased estimate may be outside the individual confidence region that was calculated with the aid of the linearized covariance matrix. The bias error and also the result of estimation are affected by the level of measurement noise (i.e., its variance), a s well as by the number and times of measurements.
SIMULATION OF THE ESTIMATION METHOD When interpreting the results of nonlinear parameter estimation care is needed, because, in general, it cannot be certain whether an obtained result is a global or local optimum and because a globally optimal result may be biased by an unknown constant whose value cannot be calculated. Also, the reliability limits obtained by linearization may be incorrect. Because of these considerations, simulation is a powerful tool for determining to what extent the results can be correct. The optimality and possible biasing can easily be revealed because the correct parameter values are known. The sample standard deviations and reliability limits that correspond to real values can also be calculated by simulating several times using independent noise every time. The Monod model is used in the simulation study because it is the best known and most frequently used model. The following parameter values
GROWTH AND SUBSTRATE CONSUMFITON MODELS
I843
were used for simulations: c,
= 0.5
S,
= 20.0 g/liter
kD = 0.2 day-’
g/liter
Y = 0.5 s, = 2.0g/liter
k, = 2.5 day-’
k, = IO.Og/liter These values are close to those given by estimation. All the initial values and parameters were estimated in this study. Modeling and measurement errors were taken into account by generating normally distributed pseudorandom noise additively into measurement as follows c,i
= C(ti)
+ w,(r,)
smi = s ( t i ) + S,
+ w&)
c ( f i )and ~ ( t are , ) the calculated dry weight and substrate concentration, respectively. w,( r ) and w,(f) are independent measurement errors. The cultivation time of the simulated batch run was five days. A total of 26 equally distributed values were “measured” from both variables. The weighting coefficients of the performance index were a, = a, = 1. Standard deviations (u)of the measurement noise vaned from 0.0 to 0.6 g/liter with steps of 0.1 g/liter, and were same for both measurements. Each standard deviation was used five times. Figure 2 presents as an example the optimal estimates of the yield Y and the maximum specific growth rate k, and the theoretical standard deviations obtained by linearization. The sample means and standard deviations were also calculated using estimation results. It is easily seen that standard deviations obtained by linearization are greater than standard deviations calculated from estimation results, i.e., from samples. The slight increasing effect of biasing can also be seen in the values of yield Y when noise level increased. The same results were also obtained starting from different initial values for iteration. DISCUSSION The presented nonlinear parameter estimation method was based on clear mathematical minimization of the cost functional, which depends on the pure unmodified measurement data and on the
I
04
7
a2
1
0
d (g/l) 0.6
14 0
0.2
0.2
d (g/L)
0.6
Fig. 2. Estimated yield Y and maximum specific growth rate k , for simulated batch runs using different levels of measurement noise with t a n d a r d deviation w. Bars represent standard deviations calculated by linearization. Dotted lines represent sample means and sample t andard deviations. ( 0 )Estimated values, (0) sample means of five estimates.
0.2
GROWTH AND SUBSTRATE CONSUMPTION MODELS
184.5
extended parameter vector to be estimated. However, nodineanties of the models made interpretation of the results impossible on the basis of any general theory of estimation. Such a theory has not yet been developed. Simulation of models and of estimation is a powerful tool when considering the possibilities of estimation based on experimental measurements and when interpreting confidence limits and reliability of estimates. The question of parameter identifiability and the effect of measurement noise on bias of estimate have been considered by simulation. Estimations using measurements from a real fermentation were successfully performed mathematically for the first three models. The convergence was obtained in four to ten steps depending on the model and on the starting values for estimation. Because of insufficient measurements, the linearized confidence limits are, however, approximate ones. Feasible parameter values were also found in the case of the enzyme kinetic model of growth, although a mathematical optimum was not obtained. The main reasons for this were the large number of parameters to be estimated compared to the number of measurements, the missing measurements of the third state variable of the model, i.e., cell-substrate complex, and especially the stiffness property of the differential equation model, which made the use of the fourth order Runge-Kutta-Gill integration method impractical. After fixing some of the model parameters, convergence was also obtained mathematically, i.e., the optimal estimate. After visual inspection the responses of all the models are, however, indistinguishable in the goodness of fit. The authors have found only few articles in the literature concerning the growth and substrate consumption of T. viride on glucose. Brown and Halstedlg considered the effect of pH on the maximum specific growth rate. The growth was exponential according to their conclusions. This means that the saturation constant of the Monod model is small compared to the values of substrate concentration, i.e., about I .O g/liter. Nagai et al.'O used an unidentified Trichoderma species growing on glucose and on cellobiose. The estimated parameters of the Monod model are compared with values found in the literature (Table 111). The values of the yield Y and the maximum specific growth rate k, are in the same range. The estimated saturation constant k, is much greater than values found in the literature. (An interesting detail is that two papers present a value of k, equal to each other.) Calculated confidence limits of the saturation constant are larger than the
1846
NIHTILA AND VIRKKUNEN
TABLE 111 Growth and Substrate Consumption Parameters of the Monod Model and Environmental Conditions from Literature Compared with Values of the Present Authors
Parameters k , (day-') k Z (&liter) Y T("C)
PH
Brown and Nagai et aLZ0 HalstedlS 3.12 0.083b 0.656d 30 6.3 + 3.5
Mitra and Wilke" 7.06 0.083 0.43 30 -3.2
-1.85 1' 0.40 28 3.0
Mandels et aLz2
Present authors
4.8" 0.3-0.4 -
2.5 1 9.3 0.468 29 3.0
Approximate values for continuous culture. 0.46mM glucose. Assumed by the present authors. * 118 g celVM glucose. a
confidence limits of the other parameters. This means that it is not possible to estimate the saturation constant accurately using only 16 measurements of both variables. In order to improve the estimates' accuracy the measurement noise should be decreased and the number of measurements increased. Using experimental design in selecting suitable measurement times the accuracy can also be improved without increasing the number of measured points. If the accuracy of one parameter estimate is to be improved, the measurements can be weighted using the inverse values of the sensitivity functions of biomass and substrate concentration with respect to the parameter to be considered. This, however, changes the original problem statement and possibly increases the confidence regions of the other parameters. APPENDIX A If it is assumed that substrate and biomass form a reversible complex that grows, the following enzyme kinetic scheme is obtained for the growth and substrate consumption process: KI
c, + s *
K3
C, 3 (1
+ u)C,
KZ
3-
0
kD
3-
kD
0
where C, is pure biomass; S is substrate; C, is the cell-substrate
GROWTH AND SUBSTRATE CONSUMPTION MODELS
1847
complex; Y is the yield coefficient; K , , K,, and K3 are the reaction rate constants, and kD is the death rate coefficient. The result of dry weight measurements is the sum of pure biomass and complex. Consequently for the respective concentrations the following measurement equation is satisfied: c
=
c,
+ cp
(26)
From the reaction scheme and taking into account the measurement equation [eq. (26)] the differential equation model for the concentrations is obtained C
=
YK3c2- kDc
S = -K,s(c C,
=
- c,)
+ K,c,
K,s(c - c,) - ( K ,
+ K3 + k , ) ~ ,
c(0) = c,
(27)
s(0)
= so
(28)
~ ~ (=0 0)
(29)
Pseudo-steady-state assumption for the cell-substrate complex leads to the equation c,
=
K,sc/(K,
+ K3 + kD + K , s )
(30)
Substitution of this into eqs. (27) and (28) results in the Monod model with parameters (c.f., eqs. ( 1 1) and (12)):
k,
= YK,
k 2
=
(K2
(3 1)
+ K3 + k D ) / K 1
(32)
APPENDIX B This method presents, in a general form, the parameter estimation of the process model dx - = f(x, p, t ) dt
(33)
y o ) = g(x(t),p, 0 + w(d
(34)
by minimizing the performance index N
J =
C [y(ti) - g(x(ti)$p , ti)]'
wi[y(ti) - g(x(ti),p, ti)]
(35)
i=l
where x is the state variable of the process, dim(x) = n ; y is the output variable, dim(y) = m; p is the parameter vector, dim(p) = np;5 is the initial value of state x at f = 0; W iare positive definite
NIHTILA AND VIRKKUNEN
I848
weighting matrices, dim(Wi) = m X m ; N is the number of measurement times, and [ IT denotes the transpose of a vector (or matrix). The solution of the differential equation [eq. (33)] depends only on the initial state and parameter vector. Consequently, it can be expressed in symbolic form x(f) = F(t, 5, p)
(36) When this is substituted into the performance index [eq. (35)], it becomes solely the function of the initial state and parameter vector
(37)
J = J(&, p)
When linearizing the output function with regard to the initial guess,
I::[ the following equation is obtained: g(x, p, t ) = go
ago aFO + --A& agT axT
+
Derivative notations are according to C ~ a k i , i.e., * ~ elements of the matrices are of the form
The superscript zero indicates that the functions are evaluated using the initial guess. Sensitivity functions are the solutions of the following linear matrix differential equations:
d-A- afo --A dt axT
A(0)
= Z(identity
(39) (40)
M(0) = 0
where dim(A) = n x n dim(M) = n
matrix)
X
n,
1849
GROWTH AND SUBSTRATE CONSUMPTION MODELS
The output function is now expressed at all time instants form g(x(t;), p, t i ) = go +
ti
in the (43)
HiAV
where Av
=
[i:]
dim(Av)
=
n
+ np
(44)
.
where t = ti, i = 1, 2, . - . , M. When substituting eq. (43) into the cost and minimizing it with respect to Av, the linear equation for Av is obtained: AAv
=
b
(46)
A is a matrix with dimensions (n + n,) X ( n + n,) and the dimension of the gradient vector b is n + np. The change Av of the estimate is calculated using Marquardt’s method, by solving the equation
(A
+ AI)Av = b
(47)
for the value X that agrees with the equation
fito+ &(A),
PO
+ Ap(A>l < J(e0,PO)
(48)
Proceeding as described and using a suitable stopping rule the optimal cost minimizing estimate
is obtained. The authors thank the personnel of the Biotechnical Laboratory, Technical Research Center of Finland, for carrying out the experimental cultivation according to the authors’ needs. This study is a part of the research project “Development, Control, and Optimization of Biotechnical Processes,” supported by the Academy of Finland and undertaken as a joint research of the Control Engineering Laboratory and the Biotechnical Laboratory.
References I . R. P. Jefferis, 111, Biochem., 10, 15 (1975). 2. R. P. Jefferis, in Fifth International Fermentation Symposium. Abstracts of Papers, H. Dellweg, Ed., Westkreuz-Druckerei und Verlag, Berlin, 1976, p. 34.
1850
NIHTILA AND VIRKKUNEN
3. D. D. Y . Ruy and A. E. Humphrey, J. Appl. Chem. Biotechnol., 23, 283 (1973). 4. W. B. Arminger, D. W. Zabriskie, and A. E. Humphrey, in Fifth International Fermentation Symposium, Abstracts of Papers, H. Dellweg, Ed., WestkreuzDruckerei und Verlag, Berlin, 1976, p. 33. 5. H. Metz, F. Wenzel, and E. Merck, in Fifth International Fermentation Symposium, Abstracrs of Papers, H . Dellweg, Ed., Westkreuz-Druckerei und Verlag, Berlin, 1976, p. 35. 6. A. Geranton and M. Colas, Automatisme, 20,452 (1975). 7. A. Meskanen, R. Lundell, and P. Laiho, Process Biochem.. 11, 31 (1976). 8. R. Spruytenburg, N. D. P. Dang, I. J . Dunn, J. R. Mor, A. Einsele, A. Fiechter, and J. R. Bourne, Chem. Eng., 310, 447 (1976). 9. F. C. Schweppe, Uncertain Dynamic Systems, Prentice-Hall, Englewood Cliffs, N. J . , 1973. 10. J. Monod, Ann. Rev. Microbiol.. 3, 371 (1949). 1 1 . J. Lescourret, Ph.D. thesis, Universite Paul Sabatier d e Toulouse, 1974. 12. D. W. Marquardt, J. Soc. Indust. Appl. Math., 11,431 (1%3). 13. Y. Bard, SIAM J . Numer. Anal., 7, 157 (1970). 14. D. M. Himmelblau, Process Analysis by Statistical Methods, Wiley, New York, 1970. 15. G. Emig and L . H. Hosten, Chem. Eng. Sci., 29,475 (1974). 16. M. Nihtila and J. Virkkunen, in Proceedings of 5th IFACIIFIP International Conference on Digital Computer Applications to Process Control, Hague, The Netherlands, 1977. 17. D. Johnsson and P. Berthouex, Biotechnol. Bioeng., 17, 557 (1975). 18. A. P. Sage, Optimum Systems Control, Prentice-Hall, Englewood Cliffs, N. J., 1968, p. 323. 19. D. Brown and D. Halsted, Biotechnol. Bioeng., 17, I 199 (1975). 20. S. Nagai, M. Onodera, and S. Aiba, Ear. J . Appl. Microbiol., 1, 9 (1976). 21. G. Mitra and C. Wilke, Biotechnol. Bioeng., 27, I (1975). 22. M. Mandels, D. Sternberg, and R. Andreotti, in Symposium on Enzymatic Hydrolysis of Cellulose, M. Bailey, T-M. Enari, and M. Linko, Eds., The Finnish National Fund for Research and Development (SITRA), Helsinki, 1975. 23. F. Csaki, Periodica Polytechnica El. (Technical University. Budapest), 15 (1971).
Accepted for Publication May 28, 1977
