[Toss-devel] Dynamics in Toss and Symbolic Runge-Kutta
Status: Beta
Brought to you by:
lukaszkaiser
|
From: Lukasz K. <luk...@gm...> - 2012-05-12 15:33:47
|
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 |