At 2:26 PM -0500 1/12/01, Pedro Mendes wrote:
>Implicit solvers are not needed explicitly because of mass conservation ,
>but rather to be stiflly stable. ODEPACK does not only have
>differential-algebraic solvers, it also has purely differential solvers
>like LSODE (the precursor to LSODA). And yes, stiff differential equations
>became an issue because of (bio)chemical dynamics.
Even the non-stiff differential equation solver in LSODA uses an
implicit algorithm (variable order Adams-Moulton). Using implicit
integration rules permits easier control of truncation error.
Stiffness arises because different variables change with different
time constants. In one of my early simulations of glycolysis, Pi
changed very slowly whereas creatine phosphate changed very rapidly.
Using a non-stiff solver required so small a time step to get
accurate integration of the creatine phosphate time course that the
change in Pi was lost in the machine round-off error and therefore
never changed over the course of the simulation.
>It turns out that in a stiff system of ODEs one can make the stiffness
>disappear by using an equilibrium approximation: the fastest reactions are
>taken to be in equilibrium and reduced to algebraic equations. This was
>Garfinkel's first approach to stiff ODEs. In doing this, one introduces the
>equilibrium assumption which carries already some approximation. No matter
>how fast a reaction is, it never proceeds at infinite speed and so it would
>be best described by a differential equation. But to do so, one ends up
>with a very stiff system, which may, nevertheless, be easy to solve with
>stiff integrators. Garfinkel's group went in this direction later on
>(possibly by the time Michael joined). Personally I favour this last
>approach. Gepasi does not accept equilibrium reactions, one is forced to
>write them explicitly and set their rate constants, approprietly, to very
>large values. This, of course, makes the system very stiff and so
>impossible to tackle with Runge-Kutta or other explicit methods...
Yes, Garfinkel's group was beginning to use the LSODE integrator
about the time I came. They had just found out how to use the
quasi-steady state approximation to solve for the enzyme forms
algebraically. A few years later we learned how to use the matrix of
stoichiometries of the fast equilibria to compute corrections to the
time derivatives of chemicals (e.g. ATP) that appear in equations for
both fast equilibria and enzymatic reactions.
> > So anyway, from my quick perusal of the literature at that time, and/or
> > from reading other people's code when available, there are some
> > home-made kinetic simulators floating around in the biochemical
> > community which "work" (or seem to), but in reality they shouldn't or
> > sometimes even don't.
>>I totally support Petr on this one. ANYONE writing code for dynamics
>should use one of the public numerical methods available. Wes hould leave
>the development of such codes to those who specialize in it. Besides, these
>public routines have been used by large numbers of people and therefore are
>reasonably optimized for speed and debugged. Why reinvent the wheel?
How true. When we were writing SCoP at Duke, a user complained that a
simulation model linked with LSODA took up too much memory on his PC.
In those days you had only a few hundred Kb for all your
applications. I wrote a stripped down Gear integrator for SCoP to
avoid this problem. It works and gives approximately the same answers
as LSODA, but it is a little slower. When you look at the truncation
error at each integration step, my program gave about 1e-10 whereas
LSODA gave about 1e-16. So the specialists have been able to optimize
the error control to minimize both the number of iterations of the
corrector formula and the single-step truncation error.
Michael C. Kohn
Laboratory of Computational Biology and Risk Analysis
National Institute of Environmental Health Sciences
P.O. Box 12233, Mail Drop A3-06
Research Triangle Park, NC 27709-2233
919-541-4929 (voice)
919-541-1470 (Fax)
E-mail: kohn at valiant.niehs.nih.gov
Web site: http://valiant.niehs.nih.gov
---
