Statement of need
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.
An example of a susceptible-infected-recovered (SIR) epidemiological model is as follows.
using Gillespie using Gadfly 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 srand(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 )
Jensen's method or uniformization
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,
x is the discrete state,
parms are the parameters, and
t is the simulation time.
The true jump method
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
- 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.
|Implementation||Time per simulation (ms)|
|Julia (Gillespie.jl, Static)||0.89|
|Julia (DifferentialEquations.jl, Static)||0.72|
(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.