Just a brief update on my activities 11 days before the R package GillespieSSA is to be released. I have been hacking away on the code, debugging, and running test simulations like there is no tomorrow. The last SSA implementation I added is the optimized tau-leaping algorithm by Cao et al. (2006) (see the reference list on the GillespieSSA page). Initially it was a bit daunting trying to wrap ones head the underlying, especially after implementing all the “straightforward” methods. For a while I was asking myself “Can I build this?” and then my two year old replied in god old Bob style “Yes we can!”. Actually, once one figures it out it is a rather ingenious implementation of the SSA.

The optimized tau-leap method (or OTL for short) choses the time increment of the simulation steps in an adaptive manner depending on the state of the system at that moment in time. There are two advantages to this method, firstly negative population sizes cannot arise (in contrast to the Poisson tau-leaping method) and secondly it is more accurate than fixed-step tau-leap methods as it switches to the direct method when the leap condition is violated (a third advantage is, of course, that it can be substantially faster than the direct method). It’s a really neat concept, switching between different methods adaptively as the simulation progresses.

Below is the trajectory of Lotka’s predator-prey prey model illustrating this point. The red point is the neutrally stable-point for the deterministic ODE and also the starting point for the simulation. The blue points is the trajectory of the populations plotted every 5 time steps. There are several interesting patterns one can see in this plot. First, while the finite population starts on the neutrally-stable equilibrium it rapidly starts to oscillate with ever increasing amplitude, inevitably leading to extinction. This nicely illustrates the important concept of deterministic models not always being able to accurately predict the dynamics of finite (and hence stochastic) populations . The second interesting feature in the plot is that can see the adaptive switches between the Poisson tau-leap method and the direct method that the optimized tau-leap method relies on, these are the regions with the dense point trajectories (lower-left corner).

By the way, Gillespie uses this model in his 1976 paper to illustrate the direct method and this is his view of the trajectory,

Thus, one way of viewing the behavior of the Lotka model is to say that the microscopic random fluctuations inherent in the reactions cause the system to execute a “drunkard’s walk” over the continuum of concentric, neutrally stable solution orbits, staggering sometimes outward and sometimes inward.