Gillespie.jl

This is an implementation of Gillespie's direct method as well as uniformization/Jensen's method for performing stochastic simulations, which are widely used in many fields, including systems biology and epidemiology. It borrows the basic interface (although none of the code) from the R library GillespieSSA by Mario Pineda-Krch, although Gillespie.jl only implements exact methods at present, whereas GillespieSSA also includes tau-leaping, etc.. It is intended to offer performance on par with hand-coded C code; please file an issue if you find an example that is significantly slower (2 to 5 times) than C.

The stable release of Gillespie.jl can be installed from the Julia REPL using the following command.

using Pkg
Pkg.add("Gillespie")

The development version from this repository can be installed as follows.

Pkg.add("https://github.com/sdwfrost/Gillespie.jl")

An example of a susceptible-infected-recovered (SIR) epidemiological model is as follows.

using Gillespie
using Gadfly
using Random

function F(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
Random.seed!(1234)

result = ssa(x0,F,nu,parms,tf)

data = ssa_data(result)

plot_theme = Theme(
    panel_fill=colorant"white",
    default_color=colorant"black"
)
p=plot(data,
    layer(x=:time,y=:x1,Geom.step,Theme(default_color=colorant"red")),
    layer(x=:time,y=:x2,Geom.step,Theme(default_color=colorant"orange")),
    layer(x=:time,y=:x3,Geom.step,Theme(default_color=colorant"blue")),
    Guide.xlabel("Time"),
    Guide.ylabel("Number"),
    Guide.title("SSA simulation"),
    Guide.manual_color_key("Subpopulation",["S","I","R"],["red","orange","blue"]),
    plot_theme
)

Julia versions of the examples used in GillespieSSA are given in the examples directory.

The development version of Gillespie.jl includes code to simulate via uniformization (a.k.a. Jensen's method); the API is the same as for the SSA, with the addition of max_rate, the maximum rate (Float64). Optionally, another argument, thin (Bool), can be set to false to return all the jumps (including the fictitious ones), and saves a bit of time by pre-allocating the time vector. This code is under development at present, and may change. Time-varying rates can be accommodated by passing a rate function with three arguments, F(x,parms,t), where x is the discrete state, parms are the parameters, and t is the simulation time.

The development version of Gillespie.jl also includes code to simulate assuming time-varying rates via the true jump method; the API is the same as for the SSA, with the exception that the rate function must have three arguments, as described above.

The speed of an SIR model in Gillespie.jl was compared to:

• A version using the R package GillespieSSA
• Handcoded versions of the SIR model in Julia, R, and Rcpp
• DifferentialEquations.jl's jump interface.

1000 simulations were performed, and the time per simulation computed (lower is better). Benchmarks were run on a Mac Pro (Late 2013), with 3 Ghz 8-core Intel Xeon E3, 64GB 1866 Mhz RAM, running OSX v 10.11.3 (El Capitan), using Julia v0.4.5 and R v.3.3. Jupyter notebooks for Julia and R with the code and benchmarks are available as gists. A plain Julia file is also provided in the benchmarks subdirectory for ease of benchmarking locally.

(smaller is better)

Julia performance for Gillespie.jl is much better than GillespieSSA, and close to a handcoded version in Julia (which is itself comparable to Rcpp); as compiler performance improves, the gap in performance should narrow.

Gillespie.jl is under development, and pull requests are welcome. Future enhancements include:

• Constrained simulations (where events are forced to occur at specific times)
• Discrete time simulation

If you use Gillespie.jl in a publication, please cite the following.

• Frost, Simon D.W. (2016) Gillespie.jl: Stochastic Simulation Algorithm in Julia. Journal of Open Source Software 1(3) doi:0.21105/joss.00042

A Bibtex entry can be found here.