cayenne: stochastic simulation of biochemical systems in Python
cayenne simulates biochemical reaction systems in which molecule counts are small enough that randomness changes the answer. Models are written as reactions and rate constants, then simulated with one of three algorithms to produce species trajectories across repeated runs. I wrote it with Srikiran Chandrasekaran, and the simulation loops are implemented in Cython.
At small copy numbers there is no average behaviour to describe
A rate equation describes the average behaviour of a large population. If the system starts with five bacteria in a dish, there is no such average: the population either dies out or begins dividing and growing. The outcome is itself random, and a differential equation returns a smooth curve through the middle of two outcomes that never occur. The same problem arises in gene expression, where the copy number of a transcript is a small integer rather than a concentration. Gillespie-type simulation gives the distribution of trajectories instead of the mean trajectory.
Benchmarks against three existing packages
Packages for stochastic simulation already existed, including Tellurium in Python, GillespieSSA in R and BioSimulator.jl in Julia, so the case for writing another one had to be demonstrated rather than asserted. We benchmarked all four against a subset of the stochastic component of the SBML test suite (DSMTS): 4 systems with different parameter combinations, 14 model-parameter pairs in total, and 10000 repetitions each. Exact solvers were scored with the Z and Y statistics recommended by the test suite. Approximate solvers were scored on the ratio of means and standard deviations against the analytical result.
Across these benchmarks, cayenne’s direct method was the most accurate of the four packages and at least an order of magnitude faster than GillespieSSA. Its interpolation also separated simulation accuracy from the irregular time points at which stochastic solvers return results.
Models are written in a subset of antimony
A model is a string in a subset of the antimony modelling language, with one deliberate difference: cayenne accepts rate constants rather than custom rate equations and generates the kinetics from them. The model supports zero-, first-, second- and third-order reactions.
Compartment volume enters the kinetics of second- and third-order reactions, and the model distinguishes chemical concentrations from biological copy counts.
Three algorithms are available.
algorithm |
Standard name | Reference | Character |
|---|---|---|---|
direct |
Gillespie’s Direct method | Gillespie 1976 | Exact, one reaction per step |
tau_leaping |
Standard tau leaping | Gillespie 2001 | Approximate, fixed step |
tau_adaptive |
Tau leaping with efficient step size selection | Cao et al. 2006 | Approximate, adaptive step |
All three simulation loops are implemented in Cython over shared propensity and reaction-selection machinery. These are the loops that run millions of times during a 10000-repetition benchmark.
Interpolation makes runs comparable to each other
A stochastic simulation returns states at random time points, such as t = 9.6 and then t = 10.2. Comparing runs with each other, or with an analytical solution, requires the state at t = 10. cayenne performs that interpolation in its backend.
Custom rate equations are rejected because they break the underlying assumptions
Refusing custom rate equations is the design choice users notice. A rate law such as the Monod equation can violate the assumptions stochastic simulation rests on, namely a well-stirred chamber of molecules in Brownian motion. The same restriction rules out writing Michaelis-Menten kinetics directly, even though that expression is grounded in chemical kinetic theory, because it is not a 0-3 order reaction. The elementary reactions underneath it are first and second order and can be modelled directly, at the cost of naming a few more constants.
Links
- Repository: https://github.com/Heuro-labs/cayenne
- Documentation: https://cayenne.readthedocs.io
- Benchmarks: https://github.com/Heuro-labs/cayenne-benchmarks
- Preprint: https://doi.org/10.1101/2020.10.10.334623
Citation
Kishore, D. and Chandrasekaran, S. (2020). Introducing and benchmarking the accuracy of cayenne: a Python package for stochastic simulations. bioRxiv 2020.10.10.334623. https://doi.org/10.1101/2020.10.10.334623