CellML Discussion List

Text archives Help


[cellml-discussion] Identifying time-points at which the model becomes stiff


Chronological Thread 
  • From: ak.miller at auckland.ac.nz (Andrew Miller)
  • Subject: [cellml-discussion] Identifying time-points at which the model becomes stiff
  • Date: Wed, 18 Apr 2007 14:14:17 +1200

Hi all,

A common issue a lot of CellML models face is that they have points at
which they go from being non-stiff (from the point of view of a
numerical algorithm), and then a sudden change happens. This change is
often due to a piecewise function (i.e. the model is actually
non-differentiable at the point).

As it is desirable to have stimulus currents and the like on for the
shortest time possible, so you can see what the model itself does
immediately after the stimulus, the only way around this at the moment
is to set the maximum step size of your integrator algorithm very low
(or, for example, tweak error control parameters to stop the step size
from getting so low).

The problem with this is that the stimulus actually only occurs for a
very small proportion of the total run time of the model, so this is
needlessly slow. If there was a way to identify the bound variable
values at which we definitely want to evaluate, we could then ensure
that the integrator stops on that exact value. This would allow us to
even have singularities (for example, supporting singularity reset rules
following the reset rule MathML pattern I recently proposed).

There are two potential approaches:
1) Attempt to automatically determine the values at which the applicable
piece in piecewise equations changes (analytic determination). This can
be done for certain restricted cases, but in general is not feasible
(and would probably require finding an analytic solution of the system
in some cases, such as when variables other than the bound variable are
involved). To get around this, we would have to define a standard set
for which model authors can expect this behaviour to work.
2) Provide a metadata specification which provides additional
information to help simulators find singularities or other special time
points.

Because of the problems with the general case in 1, I prefer the second
option. However, there are still difficulties:
i) How do you specify the points? Listing explicit points is not very
good, because it doesn't generalise well to PDEs, and there might be
infinitely many such points. The time points might also depend on
parameters, and you would want these to automatically update when you
re-parameterise the model. You would probably need to allow mathematical
expressions of some form to be referenced from the metadata, but we
would still want to provide a clean separation between data and metadata.
ii) What about the case where the modeller doesn't know analytically at
what bound variable value something happens, but only a condition
involving some other variables (for example, the yeast cell cycle model
has state variables which, when they reach a threshold, cause the cell
to divide)? In this case, it would still be useful to provide some sort
of advance warning to the simulator of the condition. For example, the
modeller could provide a function which becomes zero at special points,
but further away, is non-zero, and doesn't show the same sort of
discontinuity that the rates show.

Are there any suggestions or opinions on how this case should best be
handled?

Best regards,
Andrew





Archive powered by MHonArc 2.6.18.

Top of page