statistics on systems of ODE's
Erix Schokker
I am doing research on kinetics of enzyme inactivation.
Therefore I want to fit my data to different inactivation models.
for example a simple model:
k1
A < = > B
k2
k3
B -> C
(in which A,B and C are types of molecules and k1,k2 and k3 are reaction constants)
To solve this model a set of ODE's is made:
dA/dt = -k1.A + k2.B
dB/dt = k1.A - k2.B - k3.B
dC/dt = k3.B
This system of ordenary differential equations can not be solved analytically but it has to be done numerically (with let's say Rung=
e-Kutta or Gear). By varying the reaction constants a good fit might be found.
To find the best fit I probably have to give initial values for the reaction constants; evaluate the model using some sort of sum of=
squares of (measured y - predicted y); adapt the reaction constants etc.
My questions are:
How can I find the best values for the reaction constants kx? What kind of statistics can I use?
Is it possible in a package such as SAS to perform and numerical and statistical analysis? Can SAS or SPSS perform Runge-Kutta or ca=
n you program it yourself?
Is there ready-to-use public-domain software for evaluating different models of chemical reactions?
Thanks in advance,
Erix Schokker
Department of Food Science
Wageningen Agricultural University
The Netherlands
Jose Iribarne
: How can I find the best values for the reaction constants kx? What kind
: perform and numerical and statistical analysis?
The NONLIN procedure of SAS can be used to estimate the best fit
parameters of any model, given a set of data. The only requirement is that
you have to supply an analytic expression of the model's results. In your
case, the ordinary differential equations can be solved analytically (for
example using Laplace transforms). You also need to input expressions of
the first derivatives (at least), but you can evaluate them numerically
(calculate the effect of a small variation in the parameters).
SAS will output the least-squares values of the parameters, and an
asymptotical estimate of their standard errors. This is very useful to
determine how many significant figures you have in each.
I've used this procedure to develop apparent kinetic models and it works
very well. When you have many parameters or they are cross-correlated, the
Marquardt method option converges better than the standard Newton method.
Regards,
Jose Iribarne
SUNY-ESF
Disclaimer: these opinions are my responsibility only.
Herman Rubin
Erix Schokker <erix.s...@zuivel.lmt.wau.nl> wrote:
>Hello,
>I am doing research on kinetics of enzyme inactivation.
>Therefore I want to fit my data to different inactivation models.
>for example a simple model:
> k1
>A < = > B
> k2
> k3
>B -> C
>(in which A,B and C are types of molecules and k1,k2 and k3 are reaction constants)
>To solve this model a set of ODE's is made:
>dA/dt = -k1.A + k2.B
>dB/dt = k1.A - k2.B - k3.B
>dC/dt = k3.B
>This system of ordenary differential equations can not be solved analytically but it has to be done numerically (with let's say Rung=
>e-Kutta or Gear).
I do not have any problem in solving this system, or any other linear
system, in closed form, given the initial conditions. One can then
get a slightly complicated non-linear fitting problem.
By varying the reaction constants a good fit might be found.
>To find the best fit I probably have to give initial values for the reaction constants; evaluate the model using some sort of sum of=
> squares of (measured y - predicted y); adapt the reaction constants etc.
>My questions are:
>How can I find the best values for the reaction constants kx? What kind of statistics can I use?
>Is it possible in a package such as SAS to perform and numerical and statistical analysis? Can SAS or SPSS perform Runge-Kutta or ca=
>n you program it yourself?
>Is there ready-to-use public-domain software for evaluating different models of chemical reactions?
One can even handle systems which can only be solved numerically. The
important thing in doing this is to compute the necessary derivatives
of the solution with respect to the parameters AT THE SAME TIME as the
the solution of the differential equations. For the first derivatives,
letting _i denote the derivative with respect to ki, this becomes
dA_1/dt = -k1.A_1 + k2.B_1 - A
dA_2/dt = -k1.A_2 + k2.B_2 + B
dA_3/dt = -k1.A_3 + k2.B_3,
etc.
It is especially important to do this if the system is not deterministic,
but has random input and is simulated.
--
Herman Rubin, Dept. of Statistics, Purdue Univ., West Lafayette IN47907-1399
hru...@stat.purdue.edu Phone: (317)494-6054 FAX: (317)494-0558
John C. Nash
>
>Erix Schokker (erix.s...@zuivel.lmt.wau.nl) wrote:
>: How can I find the best values for the reaction constants kx? What kind
>: perform and numerical and statistical analysis?
>
ostensibly can estimate kinetic models. We chose a range of styles rather
than going for numbers. Because of the variety, yet the space restraint
of the journal inviting the review, the work has been heavy. Not to
reveal all before publication (in The American Statistician), I believe
that MLAB, PCNONLIN and Scientist are all products that are closer to
estimating kinetic models than SAS, since they permit numerical
integration of the model (rather than insisting on an analytic form),
though if an analytic form is available, I recommend its use. You
can develop your own solution in Matlab or S Plus, but need to be
a programmer to do so.
For Mac users, the new Data Desk 5.0 has an interesting blend of my
own numerical methods for nonlinear models with Paul Velleman's visual
fitting. At the risk of tooting our own horn, the interface is quite
original, but we don't have integration of ODEs.
The field is quite active right now, so even interested parties like
myself should not expect to know all that is going on. I welcome news
of developments at jcn...@aix1.uottawa.ca.
John C. Nash
Professor of Management
James N. Petersen
MGA distributes a package that was developed by Dow Chemical called
Simusolv which does exactly what you want to do. Is it also included in
the review mentioned below?
I have not heard of the other packages that you describe. Could you
send additional information?
James N. Petersen
Professor
Chemical Engineering Department
Washington State University
John C. Nash (jcn...@aix1.uottawa.ca) wrote:
Michael Whitbeck
Petersen) wrote:
>MGA distributes a package that was developed by Dow Chemical called
>Simusolv which does exactly what you want to do. Is it also included in
>the review mentioned below?
>
>John C. Nash (jcn...@aix1.uottawa.ca) wrote:
>: In article <487ulb$o...@newstand.syr.edu>, jiri...@newstand.syr.edu
(Jose Iribarne) says:
>: >
>: >Erix Schokker (erix.s...@zuivel.lmt.wau.nl) wrote:
>: >: How can I find the best values for the reaction constants kx? What kind
>: >: of statistics can I use? Is it possible in a package such as SAS to
>: >: perform and numerical and statistical analysis?
>: >
I didn't see the original post so this may be off the mark.
Given a chemical reaction mechanism and experimental data and you wish to
estimate one or more rate constants in the model based on the data you
might consider React3.21 (free at infomeister.osc.edu- MSDOS, Mac, OS/2).
The program reads in a reaction mechanism file and builds the descriptive set of
differential equations which are then solved subject to integration parameters
from separate parameter and initial value files.
If experimental data is available it may optionally be used to adjust one or
more rate constants in a least-squares sense.
A reaction mechanism file might look like:
Mechanism from Edelson, J Chem Ed 52, 642 (1975)
5e-8 O2- + Cs+ --> Cs + O2
1e-12 Cs+ + e- --> Cs
3.24e-3 Cs --> Cs+ + e-
.4 O2- --> O2 + e-
1.76e-16 O2 + Cs --> CsO2
1.24e-30 O2 + O2 + e- --> O2 + O2-
1.4e-16 O2 + e- --> O2-
Lines lacking an ARROW are treated as comments.
Reactions flagged with an '*' are treated as least-squares parameters
if experimental data is also provided (in a separate file, see manual).
an initial value file might look like
e- 100 1e-30
O2- 520 1e-30
Cs+ 620 1e-30
O2 3.6e14 1e-30
Cs 1.0e12 1e-30
first column is name, then amount/concentration and last an error bound term
a parameter file looks like
0 100 1e4 100
start time, time step, stop time, max number of lines to print (for safety)
Additional inputs must be specified for fitting data- see the readme files.
Output is a table suitable for use with many common plotting programs such as
gnuplot, IGor, Cricket, etc.
This program is no substitute for a knowledge of chemistry,
kinetics or modeling. Remember the first law of computing "GIGO" Garbage
In -> Garbage Out!
use instructions are in the Mac directory at infomeister (follow
pub/chemistry/software)
MW
jcn...@aix1.uottawa.ca
>
>
>MGA distributes a package that was developed by Dow Chemical called
>Simusolv which does exactly what you want to do. Is it also included in
>the review mentioned below?
No. As far as I am aware (but welcome correction if wrong), Simusolv is
not currently officially available for the Windows platform. I believe
that there are Unix versions, but my information was all word of
mouth and received some time ago. If someone can send me an email
address or web reference, I'd like to know current status, though we
won't be adding to the set of products in this review, as it has
proven very difficult to get all the tests run. Difficulties of the
syntax and configuration variety have caused a lot of frustration
and wasted time. Also vendor's want you to read "the manual", but we
have 2 LARGE boxes of them for 5 packages!
>
>I have not heard of the other packages that you describe. Could you
>send additional information?
>
MLAB is available from Civlilized Software (http://www.civilized.com)
PCNONLIN from Scientific Consulting Inc (75450...@compuserve.com)
Scientist from Micr...@Delphi.com
MATLAB from MathWorks (I believe www.mathworks.com is home page)
SPlus from StatSci (sales rep j...@statsci.com or possibly sa...@statsci.com)
John C. Nash
University of Ottawa
