|
1 | 1 | # [Advice for performant jump simulations](@id jump_simulation_performance)
|
2 |
| -We have previously described how to perform simulations of stochastic chemical kinetics *chemical reaction network* (CRN) jump process 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 stochastic chemical kinetics simulations, which we will subsequently refer to as jump process simulations. All jump process simulations arising from stochastic chemical kinetics representations of Catalyst models are performed using stochastic simulation algorithms (SSAs) from JumpProcesses.jl. Please see the [JumpProcesses documentation](https://github.com/SciML/JumpProcesses.jl) for a more extensive introduction to the package and the available solvers. |
| 2 | +We have previously introduced how to use *stochastic chemical kinetics* to perform jump simulations of *chemical reaction network* (CRRN) models (using e.g. Gillespie's algorithm). These simulations can, however, be very computationally intensive. Fortunately, there are several ways to increase their performance (thus reducing runtime). Here, we describe various considerations for performant stochastic chemical kinetics simulations, which we will here refer to as jump process simulations. All jump process simulations arising from stochastic chemical kinetics representations of Catalyst models are performed using stochastic simulation algorithms (SSAs) from JumpProcesses.jl. Please see the [JumpProcesses documentation](https://github.com/SciML/JumpProcesses.jl) for a more extensive introduction to the package and the available solvers. |
3 | 3 |
|
4 | 4 | #### Brief (and optional) introduction to jump simulations
|
5 |
| -Jump processes are continuous-time, discrete-space stochastic processes. Exact realizations of these processes can be generated using Stochastic Simulation Algorithms (SSAs), of which Gillespie's Direct method is one popular choice. In the chemical reaction modeling context the discrete-state variables typically correspond to the integer-valued number of each chemical species at each time. A system's state changes at discrete time points, the jump times, when the amount of one or more species are increased by integer amounts (for example, the creation of a new protein due to translation, or the removal of one protein due to degradation). For CRNs, these jumps correspond to the occurrence of individual reactions. Typically, the frequency of each reaction depends on its *propensity* (which in turn depends on its *rate* and *substrates*). The propensity of a reaction represents its rate law, i.e. probability per time that it occurs (also known as the associated jump process' intensity function). For example, the reaction `k, A + B --> C + D` has a propensity of $k*A(t)*B(t)$ at time $t$. See [Reaction rate laws used in simulations](@ref) for more details of what propensity function Catalyst generates for a given stochastic chemical kinetics reaction. |
| 5 | +Jump processes are continuous-time, discrete-space stochastic processes. Exact realizations of these processes can be generated using Stochastic Simulation Algorithms (SSAs), of which Gillespie's Direct method is the most well-known choice. In the chemical reaction modelling contexts, the discrete-state variables typically correspond to the integer-valued number of each chemical species at each time. A system's state changes at discrete time points (the jump times) when the amount of one (or more) species are increased by integer amount(s) (for example, the creation of a new protein due to translation, or the removal of one protein due to degradation). For CRNs, these jumps correspond to the occurrence of individual reactions. Typically, the frequency of each reaction depends on its *propensity* (which in turn depends on its *rate* and *substrates*). The propensity of a reaction (which is computed through as its *rate law*) represents the probability per time that it occurs (also known as the associated jump process' intensity function). For example, the reaction `k, A + B --> C + D` has a rate $K$ and a propensity $k*A*B$. See [Reaction rate laws used in simulations](@ref) for more details of what propensity function Catalyst generates for a given stochastic chemical kinetics reaction. |
6 | 6 |
|
7 |
| -The time of the next occurrence of some reaction, and the type of reaction that occurs, are sampled from distributions that depend on the current values of the propensity. As the latter depend on the state of the system, they must be recomputed whenever the system's state changes (for example due to a reaction occurring). Hence, jump simulations' run-times are heavily dependent on how frequently these propensities must be recomputed, and how many must be recomputed when a reaction occurs. |
| 7 | +During a simulation, after the occurrence of a reaction, the simulation algorithm computes both the time to the next reaction and which reaction will occur at this time point. Both these (reaction time and type) are computed by making random draws from probability distributions that depends on the system's reactions' propensities. As the latter depend on the state of the system, they must be recomputed whenever the system's state changes (typically due to the occurrence of a reaction). Hence, jump simulations' run times are heavily dependent on how frequently these propensities must be recomputed, and how many must be recomputed when a reaction occurs. |
8 | 8 |
|
9 |
| -Typically, propensities are recomputed only whenever a jump occurs. This means that jump simulations' run-times generally scale with the number of jumps that occur. Run-times typically become prohibitively expensive for: |
| 9 | +Typically, propensities are recomputed only at jump events. This means that jump simulations' run times generally scale with the number of jumps that occur. Run times typically become prohibitively expensive for: |
10 | 10 | - Simulations of large models (i.e. with many different species, where some species occur in large copy numbers, or where many reactions are present in a system).
|
11 | 11 | - Simulations over long time spans.
|
12 | 12 | - Simulations that are performed a large number of times.
|
@@ -53,11 +53,14 @@ we can now plot the new solution:
|
53 | 53 | ```
|
54 | 54 | 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.
|
55 | 55 |
|
| 56 | +!!! note |
| 57 | + With the default saving behavior, evaluating `sol(t)` at any time `t` within the `tspan` gives the *exact* value of the solution at that time (i.e. the vector of the exact number of each species). When using `saveat` and `save_positions = (false,false)` the solution is only saved at the selected time points, and as such `sol(t)` should not be evaluated except at times at which the solution was saved. While evaluating the solution at other times will return values they will not be the exact value of the jump process! |
| 58 | + |
56 | 59 | ## Types of jumps
|
57 | 60 | 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:
|
58 | 61 | - `MassActionJump`s: These correspond to reactions which rates are constant throughout the simulation. They are typically generated when the rate contains parameters only.
|
59 | 62 | - `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.
|
60 |
| -- `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$). |
| 63 | +- `VariableRateJump`s: These correspond to reactions which rates may change at any time during the simulation. They are typically generated when a reaction's rate depends on the time variable ($t$). |
61 | 64 |
|
62 | 65 | Here are some example reactions for the different types of jumps:
|
63 | 66 | ```@example jump_simulation_performance_1
|
|
0 commit comments