[Toss-devel] Dynamics in Toss and Symbolic Runge-Kutta

[Lukasz K. <luk...@gm...>]

Hi.

I was recently interested in some bioinformatic problems
and of course I started modelling them in Toss. This is why
we now have a cell cycle example, a simple model from Tyson
from 1991. On the way, I corrected a few bugs we had with
continuous dynamics and added animations to the JS interface.

One interesting thing about the implementation of ODEs in Toss
is that it works symbolically: you can leave free parameters, but
it is slower due to that of course. I never knew whether that was
a good decision, but now I could finally test it a bit.

The example (in examples/Cell-Cycle-Tyson-1991) results in the
following system of ODEs (using v0 ... v5 for variable names here),
assuming we leave the two significant parameters, k4 and k6, symbolic.

v0' = 0.015 + -200. * v0 * v3
v1' = -0.6 * v2 + k6 * v4
v2' = -100. * v2 + k6 * v4 + 100. * v3
v3' = -100. * v3 + 100. * v2 + -200. * v0 * v3
v4' = -k6 * v4 + 0.018 * v5 + k4 * v5 * v4 * v4
v5' = -0.018 * v5 + 200. * v0 * v3 + -k4 * v5 * v4 * v4

These two parameters, k4 and k6, range, say, from 10 to 1000 (k4)
and from 0 to 10 (k6), and normally, e.g. with k4=180, k6=1,
the system oscilates with a peak of v5 after time t=30, more or less.
One interesting question, which Francois Fages and others solve with
BIOCHAM in contraintes.inria.fr/~fages/Papers/RBFS09tcs.pdf is
whether one can find k4 and k6 such that v5 at, say, t=30, is > 0.3.

The idea of doing it symbolically is that, instead of searching through
the parameter space, one computes a conjunction of assertions over
the first-order theory of the real field and uses a solver for that. Since
Runge-Kutta equations, and we use the standard RK4 as described on
  http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
with a constant time interval (we use h=0.005 in each step) operates on
polynomials, it is easy to generate the system of polynomial equations
for a simulation of N steps. It is even linear in N and uses (N+1)*6+2
free variables: v0atK, v1atK, ..., v5atK for k=0,...,N to represent the
values of each variable vi at time step K, and the two parameters k4, k6.
I made a basic interface in TermTest to generate these constraints in
the smt2 format (used in SMT competitions), so if you first do in Toss
"make Arena/TermTest.native" and then run "./TermTest.native -steps 10"
you will get the formula for 10 steps. You can then add assertions about
k4, k6, or the final values as you wish and feed it to any QF_NRA solver.

I generated the formula for 100 steps of the above Tyson model and
added the intervals for k4, k6, and initial values (0,0,1,0,0,0). The file
in smt2 format is attached. Then I ran the best solver I know on it,
and after 15 minutes it still did nothing. This does not look promising,
because the system with only these constraints should be trivially
satisfiable. But it also shows that this kind of formulas are probably
very hard for symbolic solvers. Moreover, numercal methods and
especially evolutionary optimization techniques such as CMA-ES
 http://en.wikipedia.org/wiki/CMA-ES
 http://hal.inria.fr/docs/00/58/36/69/PDF/hansen2011impacts.pdf
give much better results here much faster. So maybe it is time to
remove some of the symbolic stuff and work more on making the
numerical part faster in Toss. At least I think I will focus more on
that in the next release. But first: 0.8 - I hope to make it soon :).

Best!

Lukasz