CC BY
62
UDC 541.14
KINFITSIM - A SOFTWARE PACKAGE TO FIT KINETIC DATA
SVIR I. B., KLYMENKO A. V., PLATZM. S.
A new software package, KinFitSim, for fitting and simulating kinetic data is presented. The KinFitSim package can be used in either chemical research or for educational purposes.
The KinFitSim package obtains the best-fit kinetic parameters - rate constants, amplitudes and others - to any selected chemical mechanism, displays the calculated and experimental absorbance versus time plots and presents the best fit data in a report to the user. The KinFitSim package has more mathematical possibilities for good kinetic and fitting simulations and a modern interface with new possibilities for users than previous packages that fit kinetic data such as KINSIM and KSIM.
The KinFitSim package consists of two general parts. The first part of the software package is a kinetic simulator (KS), which computes the time dependence of the concentrations of all species involved in a user defined reaction mechanism of any order, with defined rate constants of reactions. KS can work as a separate program. The second part of the software package is a fitting simulator (FS). FS computes the best-fit rate constants to experimental kinetic data and calculated results for a user selected reaction mechanism (in KS) by a non-linear regression procedure. FS starts it’s work after obtaining the calculated results in KS and loading the experimental data.
Kinetic Simulator (KS)
Reaction mechanism
In KS work the user first defines the reaction mechanism. Then the reaction scheme is solved to yield the corresponding ordinary differential equations (ODE) of the kinetic system.
Chemical reactions are conventionally presented by stoichiometric equations, which show the proportion of reacting species in the reaction. Thus the chemical reaction can be presented by equation [1]
Z aiYi = 0 (1)
where Yi is a chemical species, щ > 0 for reaction products and < 0 for reactants. a is a stoichiometric
vector, which can be written as а = я- p where я and p are the positive vectors which correspond to products and reactants. In general we have several chemical reactions and the corresponding matrix a . Stoichiometric vectors of reactions are the columns of matrix a.
If j -th reaction goes with rate #j , then the concentration change can be written as
y = У0 +“# , (2)
where y is the concentrations vector and yo is the initial concentrations vector.
After differentiation of equation (2), we obtained the system of ODEs for the concentrations of all species:
y' = ar, (3)
where r is a vector of reactions rates.
The rate expression can be written as
rj = kj П У?' , (4)
i
where kj is a rate constant of j -th reaction. Simulation methods
The user then specifies several parameters such as the initial concentrations of reactants and rate constants, the desired numerical method for solving the ODE system and other parameters (run time, desired accuracy, etc.) to simulate the specified reaction mechanism with the selected method. Three general numerical methods are implemented in KS for solving the ODEs [1-4]: the Euler method, the Runge-Kutta method and Gear methods, including the Adams-Moulton method. These methods can solve both stiff and non-stiffODE systems.
The Euler and Runge-Kutta methods are one-step methods used for solving non-stiff systems of ODEs. A more complicated case is a stiff system of ODEs, which can not be solved by explicit Euler and Runge-Kutta methods. The Gear method was used for such systems. This Gear method is a modification of the Adams-Moulton method.
Definition of stiff ODE system.
The initial value problem
y- Ax y)=o; y(a) = yo; x є [a, b],
is stiff in segment I c [a, b], if for each x є I when the following conditions are held
1) Re(^) < 0, i = 1,2,...,n;
2) S(x) = max Re(-li) / min Re(-li) >> 1,
i =1,...,n / i=1,...,n
where Ai are the eigenvalues of Jacobian dfj dy, which are calculated in x point.
The ratio S(x) is called “the local stiffness coefficient” ofthe initial value problem. The stiff system is a problem with large Lipshitz constants or is a problem where there is a large difference between the time constants.
The corrector iterations in the stiff-stable Gear method must be performed with the Newton method. According to Newton’s method we must iteratively solve the linear systems of algebraic equations until the required accuracy is met.
132
РИ, 2001, № 1
Fitting Simulator (FS)
The second part of the KinFitSim package is FS, which includes KS as a subprogram. The least squares method is used to obtain the best fit to experimental data. When using FS one first inputs the experimental data and a reaction mechanism is specified. The user then selects the parameters to be used in the kinetic and data simulation, and chooses the method for the minimization procedure such as the gradient-search method or the Marquardt method [5-8].
The method of least squares
The method of least squares [5] is based on the hypothesis that the optimum description of a set of data is one that minimizes the weighted sum of the squares
of the deviation of the data y from the fitting function
y(tj). FS solves the minimization problem for the goodness-of-fit factor
2 N ( І2
x = 2 wj\yj -y(tj V ^ , (5)
j=0 aeQ’ v 7
where aj indicates a unit vector in the direction of the a j coordinate axis.
In order to determine the gradient, we estimate the partial derivatives numerically:
Vr2 > aj-x4j
J j d aj Aaj ’
where Aa j is a step size by which a j is changed in order to determine the derivative.
The gradient has both magnitude and dimensions and,
if the dimensions of the various parameters a j are not
the same (which is usually the case), the components of the gradient do not have the same dimensions. We
define dimensionless parameters bj by rescaling each
of the parameters a j using the step sizes Aa j as the scaling constants, so that
Q:
t1 < a( < a{U, i = 1,,
(6)
bj =
aj
Aaj ‘
where y j are the measured variables, t j are the time values, wj is the weight of value yj , y\f j) are values
of the function calculated at tj by KS, aand a(U
J l l
are lower and upper bounds of aj variation.
The method of least squares consists of determining the values ofthe unknown parameters aj , j = 1,...,m ofthe
function y(t) that yields a minimum for the function % 2
given in equation (5). For non-linear fitting problems, there are several ways of finding this minimum value. We use either the gradient-search method [5-8] or Marquardt method [5,6] for solving the above mentioned minimization problem.
Gradient-search method
In the gradient-search method of least squares [5], all the parameters a j are incremented simultaneously,
with relative magnitudes adjusted so that the resultant direction of travel in parameter space is along the gradient (or direction of maximum variation) of ^2 . The gradient V^2 is a vector that points in the direction in which у increases most rapidly and has components
2
in parameter space equal to the rate of change of % along each axis:
V - m
j=is aj
The derivative with respect to b j then becomes
^db~Aaj * ^2 (aj+Aaj)_^2 W.
Now we define the dimensionless gradient у with unit magnitude and components
yj
d^2/Э bj
1
my 2 і \2
Т\дх2/д bj i=1V
The direction that the gradient-search method follows is the direction of steepest descent, which is the opposite of the gradient у .
The main iterative formula of the gradient-search method is
aM = Jk)+Jk) dk), k = 0,1,^,
where k is the iteration index, c№ is the step size (adjustable coefficient) in the direction ofvector dJk).
The components of the vector &/k) are
Scj = -yj Aaj .
The adjustable coefficient orA) is used at each iteration to avoid overshooting the minimum, which we search
РИ, 2001, № 1
133
for. The adjusting of or?' begins by setting its to value unity. If the inequality
T2 [ a(k+1 j<^2 [ Jk) j (7)
is met the value of c№ (and therefore the value of Jk+1) is accepted and a new iteration begins. In the
other case the value of the adjustable coefficient is divided by two and the inequality (7) is tested again. This
2
algorithm is used to ensure that function % decreases at each iteration.
The minimization procedure can drive the solution of the constraint region. A reasonable way to return to the range of parameter variation is a projection of the found
point Jk+1
onto the area Q . The new iteration formula in this case is
Jk+1 = +»*W) , k = 0,1,... (8)
where ^q(x) is a projection of the x point onto the region Q .
If q is a hyperparallelepiped (as in our case) the projection can be written as
Ыx)i
(/) ^ (/)
(l) ^ ^ (m)
Xj, aj > < Xj < aj >,
(m) ^ (m)
a} >, Xj > a}
/ 1 i
j = 1,...,m ,
(9)
Student i mine;
Date:
Solvent Temperature: r-Jcime of Precursor Laser Wavelength: Monitoring Wavelength:
Report
about KhtFitSint results
Alexey Klyiiieuko 11 05 01 15:28 13
3501
Mechanism name: Altsorbancc (350 nm)
At K 11 >£L 2. At Ki2 >C t 3. 2Bt >D
K-] K-2 K-3
Output equations: A* X\ + R*X2 + C*X?>+D*XЛ
Applied methods in KS Ruler FS: Отжііеі il -seal'd i
Initial parameters vain es
I Optimized parameters values
K+l 3R6 K+l 3233928
Ki2 1H5 К12 99380.82
K+3 1R5 K+3 93775.93
XI 1 XI 09654822
X2 0.1 X2 0 03928265
X3 1 X3 0.9735624
X4 5 X4 4.179795
Z2vi jlue before filling: 0.2991093 Xі value after fiUme: 0.0115696
Fig. 1. The report about the KinFitSim work
The Marquardt method
A more sophisticated approach to the minimization is
2
to use second partial derivatives of % to determine changes in the gradient along the search path. This
2
method is equivalent to approximating the Ж hypersurface by a parabolic surface. As a result of the
2
parabolic expansion of x one obtains a matrix equation
P = 8a-q, (10)
where p and 8a are row matrices and 7 is a symmetric matrix of order m with
Pi
5/
d aj
and Vjj
22
д X
d a j d aj
We can solve equation (10) to yield parameter increments Saj. The main iterative formula of this method is
similar to (8) with c№ = 1. If the search is already close to the minimum, this method does decrease the number of steps needed, but at the expense of more elaborate computation. This method is usually named the Gauss-Newton method.
One disadvantage of this method is that, although it converges quite rapidly to the point of minimum ^2
from points nearby, it cannot be relied on to approach the minimum with any accuracy from a point outside
the region where x hypersurface is approximately parabolic.
A convenient algorithm (Marquardt) can be obtained by increasing the diagonal terms of the matrix 7 by a factor
(l + t) . Equation (10) then becomes
, _\vjj t1 + ^) for j =j,
P = Sarj with 4jj - |^jj for j ф j _ (11)
The following algorithm of choosing 2 , was proposed by Marquardt [6]:
1. Compute xja}.
2. Start initially with 2 = 0.001.
3. Compute 8a and xja + Sa) with this choice of 2 .
4. If xja + Sa) > X2{a), increase 2 by a factor 10 and repeat step 3.
5. If X2(a + Sa) < X2( a, decrease 2 by a factor 10, consider a' = a + 8a to be the new starting point and return to step 3 substituting a' for a .
134
РИ, 2001, № 1
Input and output parameters
The KinFitSim package has a user-friendly interface for the PC, which allows convenient input-output. The user can easily obtain the calculated values for any reaction scheme and change any parameters in any stage of KinFitSim work. The KinFitSim package allows one to examine, to save and to print any graphical charts generated by KS and FS (theoretical and/or experimental). The user can correct the experimental data by hands in the special editor regime. The package produces a report, which contains the calculated and experimental graphical plots, and a list of the best-fit parameter values.
Fig. 1 shows a report of KinFitSim work: the initial and best-fit parameters. Fig. 2 demonstrates the calculated and experimental data before fitting and Fig.3 shows these curves after fitting. Fig. 2 and Fig. 3 correspond to information in the report (Fig. 1) about the chemical scheme, the applied numerical methods in KS and FS, initial and optimized parameters, and all other values.
Comparison of KinFitSim with KINSIM and KSIM
We compared the work of KinFitSim package with similar results of other packages: KIN SIM (Version 3.4 for PC, Chemical Kinetic Simulation System: B.A. Barshop, R.F. Wrenn, C. Frieden (1983) Anal. Biochem. 130: 134-145; C.T. Zimmerle, K. Patane, C. Frieden (1987) Biochemistry 26: 6545-6552; Dr. Bryce Plapp (1990), Biochemistry Department, The University of Iowa) and KSIM (Version 2.0 by Neil C. Millar, 1994).
There are four main differences between KinFitSim and other packages:
РИ, 2001, № 1
1. Limitations: KINSIM and KSIM have limitations for step count: the simulation stops when the time reaches the maximum time scale, or after 1000 steps for KSIM and 1024 steps for KINSIM, whichever happens sooner. KinFitSim can do simulations with an unlimited step count.
2. Kinetic simulating, the KINSIM package has two methods for ODE solving: Gear method and flux tolerance method (Euler method). The KSIM program uses the Runge-Kutta method and Gear method for kinetic simulation. KinFitSim has three methods for numerical simulation of ODEs: Euler, Runge-Kutta and Gear (Adams-Moulton in case of non-stiff systems) methods.
3. Fitting procedure: KS IM uses the Marquardt method for fitting to ten standard functions only (linear, exponential, hyperbolic, etc.). KINSIM uses the least squares method for fitting the calculated curves by optimization method of second order: Gauss-Newton and Marquardt methods. KinFitSim has three optimization methods for solving the non-linear regression problem: gradient-search method, Marquardt method and Gauss-Newton as a particular case of Marquardt method. The user can apply one, two or all three methods for one fitting procedure and change it easily in FS. The first order gradient-search method is often more efficient than the second order methods like Gauss-Newton or Marquardt methods. This occurs because there is less computational work for numerical approximation ofpartial derivatives of the X2 function. The gradient-search method needs m +1 calculations of
X
3 2
for gradient approximation versus — m
л— m +1 2
135
2
calculations of % for approximation of matrix of second partial derivatives in second order methods.
4. Input-output parameters: KINSIM and KSIM packages were written for the MS-DOS operating system and therefore they lack a user-friendly interface. KINSIM and KSIM require the creation a number of additional files for saving transient information in the course of the normal work of these packages. The KinFitSim package has a modern Windows interface with possibilities of easy input and output (save and print) of all parameters and all graphical plots (calculated and experimental) at any stage of work of our package.
Computing
The KinFitSim package consist of original programs written by the authors in Delphi 5 and run on a PC with Pentium III 800 MHz processor and 256 Mb of RAM.
Conclusion
The KinFitSim package offers more mathematical possibilities for good kinetic and fitting simulations and a modern interface with new possibilities for users than previous packages: KINSIM and KSIM. The KinFitSim package can be used in chemical research and for educational purposes. Students can apply and study the fitting procedure with various optimization methods and study each method using different initial, simulation parameters and various kinetic schemes.
Acknowledgements
We thank Dr. George McBane (The Ohio State University, USA) for interesting discussion of kinetic fitting problems and his kind help to the authors. Creation of KinFitSim was supported by grant from The US National Science Foundation for joint work during February of 2001 for Dr Irina Svir and Alexey 136
Klymenko in the group of Professor Matthew S. Platz in The Ohio State University.
References: 1. Современные численные методы для решения обыкновенных дифференциальных уравнений / Под ред. Дж. Холла и Дж. М. Ватта. М.: Мир, 1979. 312с. 2. БахваловН. С. Численные методы. М.: Наука, 1973. 632 с.
3. Хайрер Э, Нёрсетт С., Баннер Г. Решение обыкновенных дифференциальных уравнений. Нежесткие задачи: Пер. с англ. М.: Мир, 1990. 512 с. 4. Gear C. W. Numerical initial value problem in ordinary differential equations. Prentice-Hall, 1971. 253p. 5. Bevington P. R, Robinson D. K. Data reduction and error analysis for the physical sciences. McGraw-Hill, Inc. 1992. 6. Marquardt D.W. An algorithm for least-squares estimation of nonlinear parameters. J. Soc. Ind. Appl. Math. 1963. Vol. 11. Р. 431- 444. 7. Химмельблау Д. Прикладное нелинейное программирование. М.: Мир, 1975. 534 с. 8. Евдокимов А Г. Минимизация функций и ее применение к задачам потокораспределения в инженерных сетях. Х.: Вища шк., 1985. 280 с.
Поступила в редколлегию 15.03.2001
Irina B. Svir, Kand. Phys.-Math. Sci., Senior Scientist, doctorant of Biomedical Electronics Department Kharkov State Technical University of RadioElectronics. Scientific interests: mathematical and computer modeling of physicochemical problems. Address: KTURE, Dept. Biomedical Electronics, 14 Lenina Avenue, Kharkov 61166. E-mail: svir@kture.kharkov.ua
Alexey V. Klymenko, Master of Science in Applied Mathematics, Kharkov State Technical University of RadioElectronics. Scientific interests: numerical methods and programming. Address: KTURE, Dept. Biomedical Electronics, 14 Lenina Avenue, Kharkov — 61166.
Matthew S. Platz, Professor of Chemistry, Newman and Wolfrom Laboratory Ohio State University. Scientific interests: photochemistry, physical-organic chemistry. Address: 100 West 18th Avenue, Columbus, OH 432101185, USA.
РИ, 2001, № 1
