GillesPy2

What it simulates

GillesPy2 simulates a network of chemical reactions in a well-mixed volume. You list species, rate constants and reactions. It returns the count of every species over time.

It offers three solvers. The Gillespie direct method (SSA) draws every reaction event one at a time, so it is exact. Tau-leaping fires many reactions per step, so it is faster and approximate. The ODE solver integrates the average behaviour and has no noise.

The key assumption is mass action. Each reaction fires at a rate equal to its constant times the product of the reactant counts. The volume is well mixed, so there is no space. Counts are whole molecules, so small populations can hit zero and stay there.

Input is a model, either Python objects or an SBML file. Output is one trajectory per run, as an array with a time column and one column per species.

When to use it

• You want noise at low copy numbers, for example gene expression bursts with 10 to 100 mRNAs.
• You want to know how often a population goes extinct, or how wide the spread between runs is.
• You want to compare the stochastic mean with the ODE answer for the same model.

COPASI on biosimhub covers the same model class with a desktop-style task interface, SBML import and parameter scans. Choose it when your model already exists as SBML or you need a GUI. epydemix fits better for epidemics with age groups and contact matrices. GillesPy2 fits when you want a Python script and a choice of exact or approximate solvers.

Install it yourself

Pinned to version 1.8.3 (conda-forge, checked with a micromamba dry run, 136 MB of downloads):

micromamba create -n gp -c conda-forge gillespy2=1.8.3
micromamba activate gp

The project also documents python3 -m pip install gillespy2 --user --upgrade.

The conda package has no g++ or scons. GillesPy2 logs a warning that the C++ solvers are unavailable and falls back to the NumPy solvers. Results stay valid, only slower.

Minimal script with the Python API, the same model as the example below (run in the recipe's 1.8.3 image):

import numpy as np
from gillespy2 import Model, Species, Parameter, Reaction, NumPySSASolver

model = Model(name="LV")
model.add_parameter([Parameter(name="k_birth", expression=1.0),
                     Parameter(name="k_eat", expression=0.01),
                     Parameter(name="k_death", expression=0.5)])
prey = Species(name="Prey", initial_value=100, mode="discrete")
pred = Species(name="Predator", initial_value=50, mode="discrete")
model.add_species([prey, pred])
p = model.listOfParameters
model.add_reaction([
    Reaction(name="birth", reactants={prey: 1}, products={prey: 2}, rate=p["k_birth"]),
    Reaction(name="eat", reactants={prey: 1, pred: 1}, products={pred: 2}, rate=p["k_eat"]),
    Reaction(name="death", reactants={pred: 1}, products={}, rate=p["k_death"])])
model.timespan(np.linspace(0, 20, 201))
res = model.run(solver=NumPySSASolver, number_of_trajectories=5, seed=42)
print(res[0]["time"][-1], res[0]["Prey"][-1])

res[i] is trajectory i. res[i]["time"] and res[i]["Prey"] are arrays. The official docs describe the same pattern (results[index], then a species name).

Reading the output

The example is a Lotka-Volterra predator-prey model: prey reproduce, predators eat prey and become new predators, predators die. It starts with 100 prey and 50 predators, runs to t = 20 with 201 output points and 5 SSA trajectories, seed 42.

The recipe's helper script run_gillespy2.py is not part of GillesPy2. It reads the text model (species, parameter, reaction A + B -> C ; rate lines), builds a Model with discrete species, runs the chosen solver and writes CSV files. It caps models at 50 species and 100 reactions. The seed goes to the SSA and tau-leaping solvers only.

• run_log.txt: done: 2 species, 3 reactions, 5 trajectories, solver ssa.
• trajectories.csv (13,483 bytes): columns trajectory,time,Prey,Predator. It has 5 x 201 = 1,005 data rows. Trajectory 0 starts 0,0,100,50, then reaches 146 prey and 96 predators at t = 1.
• mean.csv (6,146 bytes): per time point, the mean and standard deviation across trajectories. At t = 0 the sd is 0. At t = 1 the mean is 137.4 prey and 99.2 predators, with sd 15.3 and 5.3.

The sd uses the population formula (np.std, divide by N). With 5 trajectories the means are rough.

Your own runs

All runs below used the recipe's image, the example model and seed 42. Only the stated option changed.

Switch to the ODE solver (--solver ode). The output has one trajectory and every sd is 0. The prey peak is 133.9 at t = 1 and the predator peak 207.8 at t = 2.5. At t = 20 the ODE gives 128.8 prey and 76.9 predators. The --trajectories value is ignored.

Use tau-leaping (--solver tau_leaping). Five trajectories give a different path from the SSA run, since the algorithm draws random numbers differently. The mean prey peak is 175 at t = 12.2, against 138.8 at t = 0.9 for SSA. At t = 20 the mean is 17.2 prey and 12.0 predators. Five runs cannot say which solver is closer, so raise the count before comparing.

Raise trajectories to 30 (--trajectories 30). The mean at t = 20 is 55.7 prey (sd 53.8) and 59.5 predators (sd 68.0). Eight of 30 runs have zero prey at t = 20. The ODE gives 128.8 prey at that time, so for this model the ODE does not describe the typical stochastic run.

For your own model, edit a copy of the .rxn file. Species names are word characters. A product side of 0 means nothing is made.

Pitfalls

• One seed, one path. Seed 42 gives the same 5 trajectories every time. Report several seeds or 30 or more trajectories before drawing conclusions.
• The ODE ignores trajectories. --solver ode always writes 1 trajectory with sd 0. Do not read that zero as precision.
• Rates use counts. Prey + Predator -> 2 Predator ; k_eat fires at k_eat * Prey * Predator, with counts of molecules. Rescale k_eat by the volume ratio if you change the volume.
• Extinction is absorbing. Once prey or predators reach 0 in a run they stay there. Trajectory 0 in the example ends with 0 prey and 0 predators.
• Slow without C++ solvers. Large counts or long times take much longer in the NumPy solvers. Try tau-leaping, or install g++ and scons to use the compiled solvers.

Further reading

• Documentation: https://gillespy2.readthedocs.io/en/latest/
• Basic usage: https://gillespy2.readthedocs.io/en/latest/getting_started/basic_usage/basic_usage.html
• Installation: https://gillespy2.readthedocs.io/en/latest/getting_started/installation/installation.html
• Source: https://github.com/StochSS/GillesPy2
• Paper: Abel JH, Drawert B, Hellander A, Petzold LR (2016), IEEE Life Sci Lett 2:35-38, https://doi.org/10.1109/LLS.2017.2652448