Herbert M Sauro wrote:
>> [...] the way modern simulators work is
> that they divide the list of variable species into two groups, a dependent
> and an independent set. The idea is that you only integrate the independent
> set and compute the dependent set from the independent set. For example for
> the simple system (i.e a cycle)
>> A -> B -> A
>> there is only one independent species, A (or B but not both). One only
> integrates the ODE with respect to A and computes B from the relation, B =
> T - A, where T equals the total mass of conserved moiety in the network.
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).
@article { Bars8334,
author = "B. A. Barshop and R. W. Wrenn and C. Frieden",
title = "Analysis of numerical methods for computer simulations
of complex kinetic processes: development of {KINSIM}",
journal = "Anal. Biochem.",
volume = "130",
pages = "134--145",
year = "1983"}
In most cases the numerical results so obtained are probably just fine.
Actually I was surprised just how good the numerical accuracy was for
the home-brew KINSIM integrator in most cases, even very "stiff" ones.
However, it seems to me, a theoretician might cringe at this approach to
truncation error. For example, a generation or two or numerical
analysts seems to have worked (mostly at the Lawrence Livermore National
Lab but also at Sandia and Los Alamos) on the ODEPACK collection of ODE
solvers:
@incollection { Hind8355,
author = "A. C. Hindmarsh",
title = "{ODEPACK}: a systematized collection of {ODE} solvers",
booktitle = "Scientific Computing",
editor = "R. S. Stepleman and others",
publisher = "North Holland",
address = "Amsterdam",
year = "1983",
pages = "55--64"}
Now, I ask myself, why is it that these folks labored so long (Hindmarsh
got started in this in the sixties working with Gear: see review
below):
@article {Byrn8701,
author = "G. D. Byrne and A. C. Hindmarsh",
title = "Stiff {ODE} solvers: a review of current and coming
attractions",
journal = "J. Comput. Physics",
volume = "70",
pages = "1--62",
year = "1987"}
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:
c the following is a simple example problem, with the coding
c needed for its solution by lsodi. the problem is from chemical
c kinetics, and consists of the following three equations..
c
c dy1/dt = -.04*y1 + 1.e4*y2*y3
c dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
c 0. = y1 + y2 + y3 - 1.
c
c on the interval from t = 0.0 to t = 4.e10, with initial conditions
c y1 = 1.0, y2 = y3 = 0.
Variables y1, y2, and y3 are concentrations in a certain reaction
mechanism, and 0.04, 1.e4, 3.e7 are kinetic constants. According to the
home-brew recipe, one could just integrate y1 and y2, and in each step
simply recompute y3 from the mass balance equation y1 + y2 + y3 = 1.
But, strictly speaking, this would be wrong. Technically the mass
balance equation is an added "algebraic constraint", making the complete
system above a Differential-Algebraic system ("not ODE" system, as L.
Petzold wrote enticingly in the title of one of her articles - the first
one I cited in my previous message). This being so, the LSODI solver
("Livermore Solver of Ordinary Differential equations - Implicit") works
entirely _differently_ than say LSODA, which, as far as I know Herbert,
you are using in your simulators.
Strictly speaking then, Differential-Algebraic systems require special
techniques for their numerical solution, namely, they require the use of
so-called _implicit_ integration algorithms. Of course when the time
available for learning relatively new areas of numerical analysis is
short, and when the numerical result is not terribly important, why not
just forget the whole issue of error control and simply recompute some
species from mass balances in each integration step. One can do that,
but this would not not the correct way about it. (By "not terribly
important" I mean "less important" by comparison with what these folks
are working on at Sandia or Lawrence Livermore: something like the
kinetics of solid fuel burning in the intercontinental ballistic missile
engine ;)
---
Years ago, when deciding on a choice of numerical methods for an
simulator I was building, I've done some numerical experiments and found
that the "mass balance error control" kludge _can_ be deceiving. Not in
the simplest system like
A --> B --> A
Here, nothing can go wrong much. But already with three species in the
mechanism, (A + B + C = T), if you integrate A+B and simply "add up" the
totals to compute C from the mass balance, I found that, in some
difficult cases, the mass-balance step size control _seemed_ to work
despite the fact that both A and B were wrong: it's just that their
errors [by default] cancel out!
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.
-- Petr
_____________________________________________________________________
P e t r K u z m i c, Ph.D. mailto:pkuzmic at biokin.com
BioKin Ltd. * Software and Consulting http://www.biokin.com
