Petr Kuzmic wrote:
> Yes, I noticed this is the way some simulators work. For example, the
> simulator in KINSIM runs in one of two separate modes: (a) using the
> Gear method for each species (e.g., both A and B), or (b) using
> essentially the Euler method (!) with a step-size controller based on
> mass-balance conservation. The step is made sufficiently small so as to
> prevent a significant violation of mass balances (within that step).
Petr, that is not what Herbert was mentioning. You only integrate A, not B.
In fact if you don't do what Herbert described you run real big risks of
numerical instability. See
Reder (1988) Metabolic control theory. A structural approach. J. theor.
Biol. 135, 175-201
If one does not do a reduction of the system, as described by Reder, the
resulting Jacobian matrix can be singular. So reducing the system to
eliminate the dependent variables is not only a sound thing to do from the
Algebra point of view, it is essential from the numerical aspect.
The approach you describe (used by KINSIM) is not numerically sound, as you
point out. I am surprised that it does not result in severe errors.
> Why is it that these folks labored so hard to create _implicit_ solvers
> for _differential-algebraic_ systems, when it became necessary to handle
> the mass conservation law of chemistry? Take a look at the solver LSODI
> (a part of ODEPACK) available from Netlib.org:
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). An yes, stiff differential equations
became an issue because of (bio)chemical dynamics.
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...
> 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?
--
Pedro Mendes
Virginia Bioinformatics Institute, Virginia Tech
1750 Kraft Drive, suite 1100, Blacksburg, VA 24060, USA
Phone: +1 540 2317411 Fax: +1 540 2311777
