Estimating Kinetic Rate Constants Using Orthogonal Polynomials and Picard's Iteration Method Robert Tanner Chemical Engineering Research and Development Department, Merck and Co., Inc., Rahway, N . J . 07065
A new technique is proposed for estimating rate constants in systems of nonlinear differential equations. The method does not require numerical integration of the differential equations nor repetiiive iterations. It is implemented b y equating the first coefficients of orthogonal least-squares polynomials with those of Picard polynomials. The orthogonal polynomials are developed b y directly fitting discrete data, while the Picard polynomials are generated from the differential equations themselves. An extrapolation scheme based upon fitting the d a t a over more encompassing time domains tends to sharpen the accuracy of the orthogonal polynomial coefficients. A simulated enzyme kinetic model is used to illustrate and develop the technique. This method shows promise as a means for starting a parameter estimating gradient search algorithm, since i t easily provides the required initial starting vector.
A new technique for estimating rate constant's in differential
equations, from kinetic data, will be developed in this paper. Since our primary interest is modeling fermentation processes, the method will be illustrated mit'li a four-parameter simulated eiizyme kinetic essmple. T o underline its generality, the scheme will also be applied to the problem of fitting two parameters t o experimental data describing the pyrolytic dehydrogenation of benzene reaction. The strategy elucidating the individual reaction constants consists of first approsimating the differential equations, with polynomials in time, using Picard's it'eratioii technique. (The coefficients of these Picard polynomials are comprised of the unknowi rate constants and initial conditions.) Then, the discrete kinetic data are fitted by least-squares orthogonal polynomials. Equating the numerical coefficients of the orthogonal polynomials with their Picard counterpart's gives us the means for estimating the unknown parameters. Essentially then, what we are doing is mapping nonlinear differential equatioiis containing parameters into a set of algebraic relationships in the parameters. I n this process, we have converted the very difficult problem of estimating constants in differential equations into the much simpler one of computing them in algebraic equalit,ies, b y sacrificing some accuracy of these unknowns. These approximate algebraic relationships become more precise as the fitting domains get smaller and smaller, and as the discrete data approach a continuum. A similar technique (Darvey, el al., 1966) using Taylor's series expansions has already been proposed for distinguishing the fine struct,ure of models describing enzyme-catalyzed chemical reactions. Their method requires, however, precision fits to very dense data. In addition, only initial time data are used to evaluate two or more of the coefficients of high degree polynomials. Our goal here is somewhat different. We do not use dense data, but, rather, we attempt to obt'ain good parameter estimates from easily obtained and fairly widely spaced data poiiits. T o this end, we evaluate only the first or second terms ln low degree polynomials, repeating the process over as many time domains as we need in order to estimate all of
the unknowi constants. We differ further from the earlier method by employing an extrapolation technique which allon-s us to increase the accuracy of our estimated unknowns obtained from widely spaced data points. This extrapolation direction, moreover, can be used to further refine the estimated values with a n analog computer. Alternatively, the parameter accuracy can be improved by a digital gradient search wliicli is initialized by the estimated parameter vector. Besides being a very simple, direct method which can save computation time over other estimation techniques, this approach does not' necessarily require dat'a for all the variables in order t'o estimate all of the unknown parameters. Other ways of determining rate constants have been reviewed by I-Iimmelblau, et al. (1967), and will not be dealt with here. In their article they propose a noiiiterative method for estimating linearly appearing parameters. Their technique's applicability is limited, however, by the requirement of time profile measurements for all of the independent concentration variables. The Method
R e shall illustrate the estimation procedure b y an enzyme kinetic example. The numerical details follow in the next section. For the example, we have chosen the classical single intermediate enzyme-substrate reaction sequence which has been found t o be so useful in biochemistry. The model is expressed a s
S+E
hi
IfES ka_ P + E k?
where S = substrate; E = enzyme (catalyst); ES = enzymesubstrate complex; P = product; k , = rate constant i; i = 1, 2, 3. This four-variable, three-rate constant, isothermal model is generally described by differential equations in keeping with the law of mass action as
'' dt
=
-kl(S)(E)
+ kZ(ES)
(S)o = S*
Ind. Eng. Chem. Fundam., Vol. 11, No. 1, 1972
1
dt dt
=
-kl(S)(E)
+ (kz + k3)(ES)
(E)o= E*
kl(S)(E)
- (k2 + ka)(ES)
(ES)o = 0
=
dt
where the parentheses denote coiiceiitrations, typically in gram moles per liter. We now use Picard’s approximating functions (Kaplan, 1958) for developiiig a series solution t o the above differential equations. Picard’s method, simply stated for the system of ordinary differential equations =
x*
is just
q(t)
= 2*
+
S, t
f(k,z*)dt’
Pt
where (Birkhoff and Rota, 1962) lim xn = z ( t ) n-
m
the exact local solution. When afft/dz, is continuous, for all i and j , given t 5 T and 0 5 t 2 T , for T 5 min ( T ( K / M ) ) , where /Iz - xo// 5 K and iz1 = SUP ( 2 , t)lI. Ilere, 2: = vector of dependent variables, xn = nth iteration of x, k = vector of parameters, and /Izll = [(XI)’ . . . . ( x ~ ) ~ ] ”length ~, in R dimensional Euclidean space. Convergence in Picard’s method is uniform, meaning that convergence is independent of t at all points of t as n -.+ a ,i.e., /lxn(t) - x ( t ) < e, for n > N , and all t, where 0 i t T and e is arbitrarily small and positive. Siiice the partial derivatives of f c are continuous for kinetic equations, it is valid to apply Picard’s method t o the example problem to give
lf
+
