|
| 1 | +# [Advice for performant jump simulations](@id jump_simulation_performance) |
| 2 | +We have previously described how to perform *jump simulations* of *chemical reaction network* (CRN) models using e.g. Gillespie's algorithm. These simulations can, however, be highly computationally intensive. Fortunately, there are several ways to increase their performance (thus reducing runtime). Here, we describe various considerations for performant jump simulations. Furthermore, all jump simulations of Catalyst models are performed using the JumpProcesses.jl package, [which documentation](https://github.com/SciML/JumpProcesses.jl) provides a more extensive description of these. |
| 3 | + |
| 4 | +#### Brief (and optional) introduction to jump simulations |
| 5 | +Jump simulations are continuous time, discrete space, simulations (where the system's state, throughout the simulation, are the discrete copy numbers of each species). The system's state changes at discrete time points (in so-called jumps). For CRNs, these jumps correspond to the occurrence of individual reactions. Here, *stochastic chemical kinetics* describes how the reactions are distributed across time. Typically, the frequency of each reaction depends on its *propensity* (which in turn depends on its *rate* and *substrates*). These depend on the state of the system, and must thus be re-computed whenever the system's state changes. Hence, jump simulations' runtimes are heavily dependent on how frequently these propensities must be recomputed. |
| 6 | + |
| 7 | +Typically, propensities are recomputed only whenever a jump occurs. This means that jump simulations' runtimes, very roughly, scale with the number of jumps. Runtimes typically become prohibitively expensive for: |
| 8 | +- Simulations of large models (either with many different species or where some species occur in large copy numbers). |
| 9 | +- Simulations over long time spans. |
| 10 | +- Simulations that are performed a large number of times. |
| 11 | + |
| 12 | +## Managing of solution saving |
| 13 | +By default, `solve` saves the value of the solution at the time of every jump. For simulations where a large number of jumps occur, this can cause memory to quickly fill up. Typically, for simulations with a large number of jumps, we want to [disable this feature](https://docs.sciml.ai/JumpProcesses/dev/tutorials/discrete_stochastic_example/#save_positions_docs) and instead set the save frequency manually. To exemplify this, let us consider a simple production/degradation model: |
| 14 | +```@example jump_simulation_performance_1 |
| 15 | +using Catalyst, JumpProcesses, Plots |
| 16 | +
|
| 17 | +rn = @reaction_network begin |
| 18 | + (p,d), 0 <--> X |
| 19 | +end |
| 20 | +u0 = [:X => 10] |
| 21 | +tspan = (0.0, 1000.0) |
| 22 | +ps = [:p => 1000.0, :d => 100.0] |
| 23 | +
|
| 24 | +dprob = DiscreteProblem(rn, u0, tspan, ps) |
| 25 | +nothing # hide |
| 26 | +``` |
| 27 | +Let us simulate this model using the default options, and plot the results. Furthermore, we use [BenchmarkTools.jl's](https://github.com/JuliaCI/BenchmarkTools.jl) `@btime` macro to measure the time it takes to plot the output solution: |
| 28 | +```@example jump_simulation_performance_1 |
| 29 | +using BenchmarkTools |
| 30 | +jprob = JumpProblem(rn, dprob, Direct()) |
| 31 | +sol = solve(jprob, SSAStepper()) |
| 32 | +@btime plot(sol) |
| 33 | +``` |
| 34 | +This simulation generates a very large number of jumps, with the solution saved and plotted at each time such a jump occurs. If the number of jumps is high enough, the memory limits may be exceeded (causing Julia to crash). Here, we do not reach that limit, however, the performance of the `plot(sol)` is affected by the number of points it must render (rendering too much in a single plot is another potential cause of crashed due to memory strain). |
| 35 | + |
| 36 | +Next, we provide the `save_positions = (false, false)` option to `JumpProblem`. This turns off the saving of the solution both before and after each jump. |
| 37 | +```@example jump_simulation_performance_1 |
| 38 | +jprob = JumpProblem(rn, dprob, Direct(); save_positions = (false, false)) |
| 39 | +nothing # hide |
| 40 | +``` |
| 41 | +However, if we were to simulate our model now, the solution would only actually be saved at its initial and final times. To remedy this, we provide the `saveat = 1.0` argument to the `solve` command, causing the solution to be saved at every `1.0`th time unit (at `t` = `0.0`, `1.0`, `2.0`, ... `999.0`, `1000.0`). |
| 42 | +```@example jump_simulation_performance_1 |
| 43 | +sol = solve(jprob, SSAStepper(),saveat=1.0) |
| 44 | +nothing # hide |
| 45 | +``` |
| 46 | +we can now plot the new solution: |
| 47 | +```@example jump_simulation_performance_1 |
| 48 | +@btime plot(sol) |
| 49 | +``` |
| 50 | +Here, we note that the time to plot the simulation was reduced (for this example, admittedly from an already low level). Furthermore, the plots look different since the trajectory is sampled at much sparser time points. |
| 51 | + |
| 52 | +## Types of jumps |
| 53 | +Each reaction in a chemical reaction network model corresponds to a possible jump of the jump simulation. These jumps can be divided into 3 categories: |
| 54 | +- `MassActionJump`s: These correspond to reactions which rates are constant throughout the simulation. They are typically generated when the rate contains parameters only. |
| 55 | +- `ConstantRateJump`s: These correspond to reactions which rates remain constant between individual jumps, but may change in response to a jump occurring. They are typically generated when the rate contains variables and parameters only. |
| 56 | +- `VariableRateJump`s: These correspond to reactions which rates may change at any time during the simulation. They are typically generated when the rate depends on the time variable ($t$). |
| 57 | + |
| 58 | +Here are some example reactions for the different types of jumps: |
| 59 | +```@example jump_simulation_performance_1 |
| 60 | +# `MassActionJump`s |
| 61 | +@reaction 1.0, X --> Y |
| 62 | +@reaction k, 2X + Y --> X2Y |
| 63 | +@reaction k1*k2+k3, 0 --> X |
| 64 | +
|
| 65 | +# `ConstantRateJump`s |
| 66 | +@reaction k*log(X), Y --> X |
| 67 | +@reaction mm(X,v,K), 0 --> X |
| 68 | +
|
| 69 | +# `VariableRateJump`s |
| 70 | +@reaction k*(1+sin(t)), 0 --> X |
| 71 | +@reaction X/t, X + Y --> XY |
| 72 | +nothing # hide |
| 73 | +``` |
| 74 | +Throughout a simulation, `VariableRateJump`s' rates require updating more frequently than `ConstantRateJump`s' rates, which in turn requires updating more frequently than `MassActionJump`s' rates. Since the performance of jump simulations is heavily affected by how frequently jump rates are computed, keeping the number of `ConstantRateJump`s and `VariableRateJump`s as small as possible is advantageous. Primarily, there exist two common cases where models are written in a way so that sub-optimal jump types are generated, both of which we describe below. |
| 75 | + |
| 76 | +#### Unnecessarily putting species in rates rather than in reactions |
| 77 | +Sometimes, a rate has been made dependent on a species, where that species instead could have been made part of the actual reaction. Consider a reaction where an enzyme ($E$) catalyses the phosphorylation of a protein ($X$) to phosphorylated form ($Xᵖ$). It can be written in two different forms: |
| 78 | +```@example jump_simulation_performance_1 |
| 79 | +r1 = @reaction k*E, X --> Xᵖ |
| 80 | +r2 = @reaction k, X + E --> Xᵖ + E |
| 81 | +nothing # hide |
| 82 | +``` |
| 83 | +These two reactions will generate identical simulations (this holds for ODE, SDE, and Jump simulations). However, while `r1` will generate a `ConstantRateJump`, `r2` will generate a `MassActionJump`. Hence, if the `r2` form is used, jump simulation performance is, at no cost, improved. Since the two forms are otherwise identical, it is always preferable to, whenever possible, put species in the reactions, rather than the rates. |
| 84 | + |
| 85 | +#### Using piecewise constant, time-dependant, function instead of callbacks. |
| 86 | +Let us consider a reaction with some input that depends on time. We assume that the input is initially $0$, but after some critical time $t=5.0$ it is increased to $1$. This can be implemented through a step function: |
| 87 | +```@example jump_simulation_performance_1 |
| 88 | +input(t) = t > 5.0 |
| 89 | +r = @reaction input(t), 0 --> X |
| 90 | +nothing # hide |
| 91 | +``` |
| 92 | +Here, the production of species $X$ is switched on at the time $t=5.0$. This reaction (which rate depends on `t`) will generate a `VariableRateJump`, which is bad for jump simulation performance. However, since the rate is piecewise constant, it can instead be implemented by setting it to a constant parameter $i$, and then use a [*callback*](@ref advanced_simulations_callbacks) to update it at the critical times. This will again generate an equivalent model, but with the reaction encoded as a `MassActionJump` (rather than a `VariableRateJump`). |
| 93 | + |
| 94 | +## Jump solver selection |
| 95 | +When creating a `JumpProblem`, a specific solver is designated using its third argument. |
| 96 | +```@example jump_simulation_performance_2 |
| 97 | +using Catalyst, JumpProcesses, Plots |
| 98 | +
|
| 99 | +rn = @reaction_network begin |
| 100 | + (p,d), 0 <--> X |
| 101 | +end |
| 102 | +u0 = [:X => 10] |
| 103 | +tspan = (0.0, 1000.0) |
| 104 | +ps = [:p => 1000.0, :d => 100.0] |
| 105 | +
|
| 106 | +dprob = DiscreteProblem(rn, u0, tspan, ps) |
| 107 | +jprob = JumpProblem(rn, dprob, Direct()) |
| 108 | +nothing # hide |
| 109 | +``` |
| 110 | +Here (as throughout most of Catalyst's documentation) we have used the `Direct()` solver (which corresponds to Gillespie's original direct method [^1][^2], also called the *stochastic simulation algorithm*). This method was originally published in 1976, and since then, many additional methods for performing jump simulations of CRN models have been developed. |
| 111 | + |
| 112 | +Gillespie's direct method will, after a jump has been performed, recompute the rates of *all* possible jumps in the system. This is typically not required. E.g. consider the following system: |
| 113 | +```@example jump_simulation_performance_2 |
| 114 | +@reaction_network begin |
| 115 | + k1, X1 --> X2 |
| 116 | + k2, X2 --> X3 |
| 117 | + k3, X3 --> 0 |
| 118 | +end |
| 119 | +``` |
| 120 | +Here, the rate of the `k1, X1 --> X2` and `k2, X2 --> X3` reactions does not depend on the amount of $X3$ in the system. Hence, their rates are unaffected by the occurrence of the `k3, X3 --> 0` reaction. Performant jump simulation methods have clever ways to determine which rates require recomputing after the occurrence of each reaction, which improves their performance. Many of these depend on so-called dependency graphs (which track which reactions' rates are affected by the occurrence of which reactions). Catalyst automatically builds such dependency graphs, which means that most jump simulators can be used without any additional input. |
| 121 | + |
| 122 | +A full list of jump simulation algorithms implemented by JumpProcesses can be found [here](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation). Generally, `RSSA()` (the rejection SSA method [^3][^4]) is recommended for small models, with `RSSACR()` (the rejection SSA with composition-rejection method [^5]) typically being more performant for larger models. |
| 123 | + |
| 124 | + |
| 125 | +## Hybrid simulations |
| 126 | +For some models, copy numbers may vary greatly between different species. E.g. consider a genetic promoter which can either be in an inactive form ($Pᵢ$) or an active form ($Pₐ$). The active promoter produces a molecule ($M$): |
| 127 | +```@example jump_simulation_performance_3 |
| 128 | +using Catalyst # hide |
| 129 | +rn = @reaction_network begin |
| 130 | + (kA,kI), Pᵢ <--> Pₐ |
| 131 | + p, Pₐ --> Pₐ + M |
| 132 | + d, M --> ∅ |
| 133 | +end |
| 134 | +``` |
| 135 | +Let us simulate this model and consider the copy numbers of each individual component: |
| 136 | +```@example jump_simulation_performance_2 |
| 137 | +using JumpProcesses, Plots # hide |
| 138 | +u0 = [:Pᵢ => 1, :Pₐ => 0, :M => 10000] |
| 139 | +tspan = (0.0, 5000.0) |
| 140 | +p = [:kA => 0.05, :kI => .01, :p => 500.0, :d => 0.0005] |
| 141 | +
|
| 142 | +dprob = DiscreteProblem(rn, u0, tspan, p) |
| 143 | +jprob = JumpProblem(rn, dprob, RSSA()) |
| 144 | +sol = solve(jprob, SSAStepper()) |
| 145 | +
|
| 146 | +plt_1 = plot(sol; idxs=2) |
| 147 | +plt_2 = plot(sol; idxs=3) |
| 148 | +plot(plt_1, plt_2, size=(1200,400)) |
| 149 | +``` |
| 150 | +We note that the copy numbers $Pₐ$ are highly stochastic. Meanwhile, $M$ exists in so large copy numbers that its trajectory can be considered approximately deterministic. For models like this one, jump simulations are required to capture the stochastic behaviour of the promoter. However, the large number of jumps generated by $M$ makes such simulations excessively expensive. Here, *hybrid simulations* can be used. Hybrid simulations employ different simulation strategies for their different components. In our example, the state of the promoter could be simulated using a jump simulation, while $M$ could be simulated using an ODE. |
| 151 | + |
| 152 | +Hybrid simulations for Catalyst models are currently not supported. However, it is a work in progress. If you require such simulations, please [raise an issue](https://github.com/SciML/Catalyst.jl/issues) and we can notify you of the current state of our implementation. Alternatively, other CRN modelling tools that support hybrid simulations exist. |
| 153 | + |
| 154 | + |
| 155 | +--- |
| 156 | +## References |
| 157 | +[^1]: [D. T. Gillespie, *A general method for numerically simulating the stochastic time evolution of coupled chemical reactions*, Journal of Computational Physics (1976).](https://www.sciencedirect.com/science/article/abs/pii/0021999176900413) |
| 158 | +[^2]: [D. T. Gillespie, *Exact Stochastic Simulation of Coupled Chemical Reactions*, The Journal of Physical Chemistry (1977).](https://pubs.acs.org/doi/10.1021/j100540a008) |
| 159 | +[^3]: [V. H. Thanh, C. Priami and R. Zunino, *Efficient rejection-based simulation of biochemical reactions with stochastic noise and delays*, Journal of Chemical Physics (2014).](https://pubmed.ncbi.nlm.nih.gov/25296793/) |
| 160 | +[^4]: [V. H. Thanh, R. Zunino and C. Priami, *On the rejection-based algorithm for simulation and analysis of large-scale reaction networks*, Journal of Chemical Physics (2015).](https://pubmed.ncbi.nlm.nih.gov/26133409/) |
| 161 | +[^5]: [V. H. Thanh, R. Zunino, and C. Priami, *Efficient constant-time complexity algorithm for stochastic simulation of large reaction networks*, IEEE/ACM Transactions on Computational Biology and Bioinformatics (2017).](https://pubmed.ncbi.nlm.nih.gov/26890923/) |
0 commit comments