Benchmark ODE solver: GSL V/s Boost Odeint library

For our neural simulator, MOOSE, we use GNU Scientific Library (GSL) for random number generation, for solving system of non-linear equations, and for solving ODE system.

Recently I checked the performance of GSL ode solver V/s Boost ode solver; both using Runge-Kutta 4. The boost-odeint outperformed GSL by approximately by a factor of 4. The numerical results were same. Both implementation were compiled with -O3 switch.

Below are the numerical results. In second subplot, a molecule is changing is concentration ( calcium ) and on the top, another molecule is produced (CaMKII). This network has more than 30 reactions and approximately 20 molecules.

GSL took approximately 39 seconds to simulate the system for 1 year. While Boost-odeint took only 8.6 seconds. (This turns out to be a corner case)

Update: If I let both solvers choose the step size by themselves for given accuracy, boost usually outperforms the GSL solver by a factor of 1.2x to 2.9x. These tests were done during development of a signaling network which has many many reactions. Under no circumstances I have tested, BOOSE ODE solver was slower than GSL solver.

Advertisements

Simulating Random Walks using Langevin Equation

Random walks (Brownian motions), in addition to their theoretical potency (describes macro-scale behavior of gas starting with micro-scale description), also describes behavior of many processes in nature. A few of them; genetic network, protein expression caused by mRNA,  have been described/predicted well using stochasticity. Moreover they are ideal noise sources cause by thermal fluctuation and often found in natural processes. To get the solid foundation in this subject, see the classic “Stochatic Processes in Physics and Chemistry” by Van Kampen, and also “handbook of stochastic methods” by Gandiner.

Random walks (in fact any stochastic process) can be described by Fokker Planck equation: it describes how probability density function evolves over time. An equivalent is Master Equation which are much easier to visualize and solve (using Gillespie algorithm, a variant of Markov method). Master equation can describe “almost” all of the chemistry. In fact, Einstein built his theory of Brownian motion by writing down a variant of Fokker Planck equation. After Einstein published his work, a Frenchmen Paul Langevian discovered the same theory using a totally different approach; which is “infinitely simpler” than Einstein approach. I’d highly recommend to read the paper which does not require more than intermediate mathematics.  http://scitation.aip.org/content/aapt/journal/ajp/65/11/10.1119/1.18725

In this note, I’d present a Python recipe to solve Langevian equation to simulate random walk. The Langevian approach is computationally extremely cheap and works extremely well in practice. Of course there is a cost involve. Fokker-Planck gives you “exact” mean and variance (if you can solve it), you need to produce many Langevian trajectories to see mean and variance converging to fixed values.

But for simulating biological and physical processes, you don’t worry too much about overall mean and variance. One trajectory is good enough for introducing noise sources.

Langevian equation looks something like the following.

dx = −f(x) dt+\alpha g(x) \sqrt{dt}

where α is normally distributed with mean 0 and variance 1. Forget f(x) and g(x), the out-of-place thing about these equation is square root of dt on the right-hand side. This led to fractional calculus, and stochastic differential equations.

For the sake of “web and coding”, problem statement and python recipe which simulates this equation can be found here. 5 model trajectories of Random walk in 1D generated by this equation are attached with this note. Many others can be generated using the script solve.py. Mean is as usual; and standard deviation relates with diffusion coefficient.

PS: File solve.py implements the assignment. To produce these trajectories run the following command:

     $ python solve.py 1

To get more about “Randomness in Biology” visit http://courses.ncbs.res.in/mod/resource/view.php?id=374 (login as guest)

NOTES:

  1. Google for: Sriram Ramaswamy’s Resonance article on Einstein’s derivation. Nice readable article. Do check the Langevin original paper linked above.
  2. It is instructive for those who are interested in Control theory to read about biological processes: “How cell control its size” and other macro-molecules inside it in such. Biology is unexplored gold-mine for Control theory. Also its worth thinking how noise is minimized in cellular processes.

Rotation 1, Week 0: Understanding the problem

Its not much one can do in a 6 week rotation in lab, but when Madan Rao asked me to read some literature and if possible do something, I suggested that I will rather work on a problem and read if need arise. Reading without working is pretty stupid thing to do at my age. It might have some use for an undergraduate, at least for his/her exam. In nutshell, just reading  is another way of time-passing without feeling guilty about it. Knowledge is overrated when one wants to work with fundamental ideas.

Fortunately Madan saw it clearly that I don’t want to become Mr. Know-It-All or some sort of Pundit or Mahant with eye on a grand problem but rather a journeyman who solves problem at hand with whatever best tools available at his disposal and build his expertise by practicing his craft.

He politely suggested me to read a paper over coffee. Not a month ago, I was talking to Somya Mani about my bird-eye view of networks in biology and what I lack to deal with them: how to inject probabilistic variation into well-defined and well-behaved system. And to these ends, I want to take course offered by Mukund Thattai this semester at NCBS Bangalore, “Randomness in Biology”. After going through the paper which Madan suggested to read before further discussion what can be done during my rotation, I am happy about the problem I encountered by chance: how cells control a variable when input is mixed with random noise!

****

Cells are always trying to control certain processes. Some variables such as cell-size is kept constant in a very variable environment. One can formulate this as a control-system problem and ask what a cell can and can’t do to suppress noise or variability in a parameter under control. Bounds on the performance of a control system under noisy conditions are well-studies by many early pioneers in system science; but this is somewhat a different problem. Since molecules inside cells go through random birth and deaths (I am still not quite clear about the difference).

This paper [1] establishes bounds on networks with negative feedback. This is interesting because negative feedback is often used to minimize the influence of noise on variable under control.

In this paper, authors consider a class of system, which is generic enough to describe many biological networks, and build a mathematical framework which can tell us what a cell (or a network) can not do when a given amount of noise is present in their environment. General structure of system is captured by following three equations: x_1 \xrightarrow{u(x_2(-\infty,t))} x_1 + 1 , x_1 \xrightarrow{x_1/\tau_1} x_1 - 1 , x_2 \xrightarrow{f(x_1)} x_2 + 1

Variable x_1 is under control whose production is controlled by a control daemon u who knows past and present of x_1. Variable x_1, in turns, contol the production of x_2 since rate of its production is controlled by a function f of x_1. \tau is mean death rate of x_1. System under investigation is a system with negative feedback and signaling. Authors describe part of system using continuous differential equations but keep the controller and singaling discrete. They claim that It is extremely hard for a network to reduce noise where signal x_2 is made less frequently than the controlled component i.e. x_1. Reducing the standard deviation of x_1 tenfold can be achieved by increasing the birth rate of x_2 by a factor of 10,000.

In short, authors have developed a framework under which one can say something about what cell can not do under a noisy condition.

I am still trying to understand the theory they have used to established these bounds (see the supplementary material available in reference 1). I am also writing a simulator using SystemC/C++ to play with this idea.

REFERENCES

[1] Lestas, Ioannis, Glenn Vinnicombe, and Johan Paulsson. 2010. “Fundamental Limits on the Suppression of Molecular Fluctuations.” Nature 467 (7312) (Sep 09):Lestas, Ioannis, Glenn Vinnicombe, and Johan Paulsson. 2010. “Fundamental
Limits on the Suppression of Molecular Fluctuations.” Nature 467 (7312) (Sep 09): 174–178. doi:10.1038/nature09333. http://dx.doi.org/10.1038/nature09333.