Model Simulation Introduction

Catalyst's core functionality is the creation of chemical reaction network (CRN) models that can be simulated using ODE, SDE, and jump simulations. How such simulations are carried out has already been described in Catalyst's introduction. This page provides a deeper introduction, giving some additional background and introducing various simulation-related options.

Here we will focus on the basics, with other sections of the simulation documentation describing various specialised features, or giving advice on performance. Anyone who plans on using Catalyst's simulation functionality extensively is recommended to also read the documentation on solution plotting, and on how to interact with simulation problems, integrators, and solutions. Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for ODEs or SDEs.

Background to CRN simulations

This section provides some brief theory on CRN simulations. For details on how to carry out these simulations in actual code, please skip to the following sections.

CRNs are defined by a set of species (with the amounts of these determining the system's state during simulations) and a set of reaction events (rules for how the state of the system changes). In real systems, the species amounts are discrete copy-numbers, describing the exact numbers of each species type present in the system (in systems biology this can e.g. be the number of a specific molecule present in a cell). Given rates for these reaction events, stochastic chemical kinetics provides a formula for simulating the system that recreates its real reaction process. During stochastic chemical kinetics simulations, the system's state is defined by discrete copy-numbers (denoting the number of each species present in the system). Next, at the occurrence of individual reaction events, the system's state is updated according to the occurred reaction. The result is a stochastic process. The most well-known approach for simulating stochastic chemical kinetics is Gillespie's algorithm.

In practice, these jump simulations are computationally expensive. In many cases, copy-numbers are so large that they can be approximated as continuous concentrations, and the time-development of the system as a deterministic process. This creates an ordinary differential equation (ODE), and is the chemical reaction network form most people are most familiar with. The rule for how ODEs are generated from CRNs is called the reaction rate equation (RRE).

Here, the RRE enables fast, approximate, and deterministic simulations of CRNs, while stochastic chemical kinetics enables exact, stochastic, simulations of the true process. An intermediary approach is to use the chemical Langevin equation (CLE) to formulate a stochastic differential equation (SDE). This approximates the system's state as continuous concentrations, but does not assume that its time development is deterministic. Generally, the CLE is used when copy-numbers are large enough that the continuous approximation holds, but not so large that the system's behaviour is deterministic. Generally, the advantage of SDE simulations (compared to jump ones) is that they are faster. Also, since the system state is continuous, interpretation of e.g. stability and steady state results from the deterministic (also continuous) domain is easier for SDEs (however one should be careful when making such interpretations).

These three different approaches are summed up in the following table:

Performing (ODE) simulations

The following section gives a (more complete introduction of how to simulate Catalyst models than our previous introduction). This is illustrated using ODE simulations (some ODE-specific options will also be discussed). Later on, we will describe things specific to SDE and jump simulations. All ODE simulations are performed using the OrdinaryDiffEq.jl package, which full documentation can be found here. A dedicated section giving advice on how to optimise ODE simulation performance can be found here

To perform any simulation, we must first define our model, as well as the simulation's initial conditions, time span, and parameter values. Here we will use a simple two-state model:

using Catalyst
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
u0 = [:X1 => 100.0, :X2 => 200.0]
tspan = (0.0, 5.0)
ps = [:k1 => 2.0, :k2 => 5.0]

To simulate the model we first bundle these up into an ODEProblem:

oprob = ODEProblem(two_state_model, u0, tspan, ps)

Next, we can simulate the model (requires loading the OrdinaryDiffEq.jl package). Simulations are performed using the solve function.

using OrdinaryDiffEqDefault
sol = solve(oprob)

Finally, the result can be plotted using the Plots.jl package's plot function:

using Plots
plot(sol)

More information on how to interact with solution structures is provided here and on how to plot them here.

Some additional considerations:

• If a model without parameters has been declared, only the first three arguments must be provided to ODEProblem.
• While the first value of tspan will almost always be 0.0, other starting times (both negative and positive) are possible.

Designating solvers and solver options

While good defaults are generally selected, OrdinaryDiffEq enables the user to customise simulations through a long range of options that can be provided to the solve function. This includes specifying a solver algorithm, which can be provided as a second argument to solve (if none is provided, a suitable choice is automatically made). E.g. here we specify that the Rodas5P method should be used:

using OrdinaryDiffEqRosenbrock
sol = solve(oprob, Rodas5P())

A full list of available solvers is provided here, and a discussion on optimal solver choices here.

Additional options can be provided as keyword arguments. E.g. the maxiters arguments determines the maximum number of simulation time steps (before the simulation is terminated). This defaults to 1e5, but can be modified through:

sol = solve(oprob; maxiters = 1e4)

Here follows a list of solver options which might be of interest to the user.

• adaptive: Toggles adaptive time stepping for valid methods. Default to true.
• dt: For non-adaptive simulations, sets the step size (also sets the initial step size for adaptive methods).
• saveat: Determines the time points at which the simulation is saved. E.g. for saveat = 2.0 the simulation is saved every second time unit. If not given, the solution is saved after each time step.
• save_idxs: Provides a vector of species whose values should be saved during the simulation. E.g. for save_idxs = [:X1], only the value of species $X1$ is saved.
• maxiters: The maximum number of time steps of the simulation. If this number is reached, the simulation is terminated.
• seed: Sets a seed for stochastic simulations. Stochastic simulations with the same seed generate identical results.

A full list of solver options can be found here.

Alternative problem input forms

Throughout Catalyst's documentation, we typically provide initial condition and parameter values as vectors. However, these can also be provided as tuples:

u0 = (:X1 => 100.0, :X2 => 200.0)
tspan = (0.0, 5.0)
ps = (:k1 => 2.0, :k2 => 5.0)
oprob = ODEProblem(two_state_model, u0, tspan, ps)

or dictionaries:

u0 = Dict([:X1 => 100.0, :X2 => 200.0])
tspan = (0.0, 5.0)
ps = Dict([:k1 => 2.0, :k2 => 5.0])
oprob = ODEProblem(two_state_model, u0, tspan, ps)

The forms used for u0 and ps does not need to be the same (but can e.g. be a vector and a tuple).

Throughout Catalyst's documentation, we typically provide the time span as a tuple. However, if the first time point is 0.0 (which is typically the case), this can be omitted. Here, we supply only the simulation endpoint to our ODEProblem:

tend = 5.0
oprob = ODEProblem(two_state_model, u0, tend, ps)

Performing SDE simulations

Catalyst uses the StochasticDiffEq.jl package to perform SDE simulations. This section provides a brief introduction, with StochasticDiffEq's documentation providing a more extensive description. By default, Catalyst generates SDEs from CRN models using the chemical Langevin equation. A dedicated section giving advice on how to optimise SDE simulation performance can be found here.

SDE simulations are performed in a similar manner to ODE simulations. The only exception is that an SDEProblem is created (rather than an ODEProblem). Furthermore, the StochasticDiffEq.jl package (rather than the OrdinaryDiffEq package) is required for performing simulations. Here we simulate the two-state model for the same parameter set as previously used:

using Catalyst, StochasticDiffEq, Plots
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
u0 = [:X1 => 100.0, :X2 => 200.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 5.0]

sprob = SDEProblem(two_state_model, u0, tspan, ps)
sol = solve(sprob, STrapezoid())
plot(sol)

we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations.

Common SDE simulation pitfalls

Next, let us reduce species amounts (using remake), thereby also increasing the relative amount of noise, we encounter a problem when the model is simulated:

sprob = remake(sprob; u0 = [:X1 => 0.33, :X2 => 0.66])
sol = solve(sprob, STrapezoid())

┌ Warning: dt(2.220446049250313e-16) <= dtmin(2.220446049250313e-16) at t=0.7083659063508818, and step error estimate = 2.5896054312696957. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ SciMLBase ~/.julia/packages/SciMLBase/Ha7rZ/src/integrator_interface.jl:612

Here, we receive a warning that the simulation was terminated. next, if we plot the solution:

plot(sol)

we note that the simulation didn't reach the designated final time point ($t = 1.0$). In this case we also note that species concentrations are very low (and sometimes, due to the relatively high amount of noise, even negative). This, combined with the early termination, suggests that we are simulating our model for too low species concentration for the assumptions of the CLE to hold. Instead, jump simulations should be used.

Next, let us consider a simulation for another parameter set:

sprob = remake(sprob; u0 = [:X1 => 100.0, :X2 => 200.0], p = [:k1 => 200.0, :k2 => 500.0])
sol = solve(sprob, STrapezoid())

plot(sol)

Again, the simulation is aborted. This time, however, species concentrations are relatively large, so the CLE might still hold. What has happened this time is that the accuracy of the simulations has not reached its desired threshold. This can be deal with by reducing simulation tolerances:

sol = solve(sprob, STrapezoid(), abstol = 1e-1, reltol = 1e-1)
plot(sol)

SDE simulations with fixed time stepping

StochasticDiffEq implements SDE solvers with adaptive time stepping. However, when using a non-adaptive solver (or using the adaptive = false argument to turn adaptive time stepping off for an adaptive solver) a fixed time step dt must be designated. Here we simulate the same SDEProblem which we struggled with previously, but using the non-adaptive EM solver and a fixed dt:

sol = solve(sprob, EM(); dt = 0.001)
plot(sol)

We note that this approach also enables us to successfully simulate the SDE we previously struggled with.

Generally, using a smaller fixed dt provides a more exact simulation, but also increases simulation runtime.

Scaling the noise in the chemical Langevin equation

When using the CLE to generate SDEs from a CRN, it can sometimes be desirable to scale the magnitude of the noise. This can be done by introducing a noise scaling term, with each noise term generated by the CLE being multiplied with this term. A noise scaling term can be set using the @default_noise_scaling option:

two_state_model = @reaction_network begin
    @default_noise_scaling 0.1
    (k1,k2), X1 <--> X2
end

Here, we set the noise scaling term to 0.1, reducing the noise with a factor $10$ (noise scaling terms $>1.0$ increase the noise, while terms $<1.0$ reduce the noise). If we re-simulate the model using the low-concentration settings used previously, we see that the noise has been reduced (in fact by so much that the model can now be simulated without issues):

u0 = [:X1 => 100.0, :X2 => 200.0]
tspan = (0.0, 1.0)
ps = [:k1 => 200.0, :k2 => 500.0]
sprob = SDEProblem(two_state_model, u0, tspan, ps)
sol = solve(sprob, STrapezoid())
plot(sol)

The @default_noise_scaling option can take any expression. This can be used to e.g. designate a noise scaling parameter:

two_state_model = @reaction_network begin
    @parameters η
    @default_noise_scaling η
    (k1,k2), X1 <--> X2
end

Now we can tune the noise through $η$'s value. E.g. here we remove the noise entirely by setting $η = 0.0$ (thereby recreating an ODE simulation's behaviour):

u0 = [:X1 => 0.33, :X2 => 0.66, :η => 0.0]
tspan = (0.0, 1.0)
ps = [:k1 => 2.0, :k2 => 5.0]
sprob = SDEProblem(two_state_model, u0, tspan, ps)
sol = solve(sprob, STrapezoid())
plot(sol)

It is possible to designate specific noise scaling terms for individual reactions through the noise_scalingreaction metadata. Here, CLE noise terms associated with a specific reaction are multiplied by that reaction's noise scaling term. Here we use this to turn off the noise in the $X1 \to X2$ reaction:

two_state_model = @reaction_network begin
    k1, X1 --> X2, [noise_scaling = 0.0]
    k2, X2 --> X1
end

If the @default_noise_scaling option is used, that term is only applied to reactions withoutnoise_scaling metadata.

While the @default_noise_scaling option is unavailable for programmatically created models, the set_default_noise_scaling function can be used to achieve a similar effect.

Performing jump simulations using stochastic chemical kinetics

Catalyst uses the JumpProcesses.jl package to perform jump simulations. This section provides a brief introduction, with JumpProcesses's documentation providing a more extensive description.

Jump simulations are performed using so-called JumpProblems. Unlike ODEs and SDEs (for which the corresponding problem types can be created directly), jump simulations require first processing inputs into a correct format creating an intermediary JumpInputs. In this example, we first declare our two-state model and its initial conditions, time span, and parameter values.

using Catalyst, JumpProcesses, Plots
two_state_model = @reaction_network begin
    (k1,k2), X1 <--> X2
end
u0 = [:X1 => 5, :X2 => 10]
tspan = (0.0, 5.0)
ps = [:k1 => 2.0, :k2 => 5.0]

Next, we process these into a JumpInputs:

jinput = JumpInputs(two_state_model, u0, tspan, ps)

This is then used as input to a JumpProblem:

jprob = JumpProblem(jinput)

The JumpProblem can now be simulated using solve (just like any other problem type).

sol = solve(jprob)

If we plot the solution we can see how the system's state does not change continuously, but instead in discrete jumps (due to the occurrence of the individual reactions of the system).

using Plots
plot(sol)

Designating aggregators and simulation methods for jump simulations

Jump simulations (just like ODEs and SDEs) are performed using stochastic simulation algorithms (SSAs) to generate exact samples of the underlying jump process. In JumpProcesses.jl and Catalyst, we call SSAs aggregators. These methods determine the time until, and type of, the next reaction in a system. A separate time-stepping method is then used to actually step from one reaction instance to the next.

Several different aggregators are available (a full list is provided here). To designate a specific one, provide it as the second argument to the JumpProblem. E.g. to designate that the sorting direct method (SortingDirect) should be used, use:

jprob = JumpProblem(jinput, SortingDirect())

Especially for large systems, the choice of aggregator can dramatically impact simulation performance.

Jump simulations where some rate depends on time

For some models, the rate terms of reactions may explicitly depend on time. E.g. consider the following circadian clock (inspired) model, where the production rate of some protein ($P$) depends on a sinusoid function:

circadian_model = @reaction_network begin
    A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P
    d, P --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[d]{\frac{1}{2} A \left( 1 + \sin\left( - \phi + 6.283185307179586 f t \right) \right)} \mathrm{P} \end{align*} \]

This type of model will generate so called variable rate jumps (VariableRateJumps in JumpProcesses.jl). Such models can be simulated in Catalyst too, but note that now a method for time-stepping the solver must be provided to solve. Here ODE solvers should be given as they are used to handle integrating the explicitly time-dependent propensities for problems with variable rates, i.e. the proceeding example can be solved like

using OrdinaryDiffEqTsit5
u0map = [:P => 0]
pmap = [:f => 1.0, :A => 2.0, :ϕ => 0.0, :d => 1.0]
tspan = (0.0, 24.0)
jinputs = JumpInputs(circadian_model, u0map, tspan, pmap)
jprob = JumpProblem(jinputs)
sol = solve(jprob, Tsit5())  # use the Tsit5 ODE solver to time-step
plot(sol; idxs = :P, lw = 2)

Citation

When you simulate Catalyst models in your research, please cite the corresponding paper(s) to support the simulation package authors. For ODE simulations:

@article{DifferentialEquations.jl-2017,
 author = {Rackauckas, Christopher and Nie, Qing},
 doi = {10.5334/jors.151},
 journal = {The Journal of Open Research Software},
 keywords = {Applied Mathematics},
 note = {Exported from https://app.dimensions.ai on 2019/05/05},
 number = {1},
 pages = {},
 title = {DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia},
 url = {https://app.dimensions.ai/details/publication/pub.1085583166 and http://openresearchsoftware.metajnl.com/articles/10.5334/jors.151/galley/245/download/},
 volume = {5},
 year = {2017}
}

For SDE simulations:

@article{rackauckas2017adaptive,
  title={Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory},
  author={Rackauckas, Christopher and Nie, Qing},
  journal={Discrete and continuous dynamical systems. Series B},
  volume={22},
  number={7},
  pages={2731},
  year={2017},
  publisher={NIH Public Access}
}

For jump simulations:

@misc{2022JumpProcesses,
  author       = {Isaacson, S. A. and Ilin, V. and Rackauckas, C. V.},
  title        = {{JumpProcesses.jl}},
  howpublished = {\url{https://github.com/SciML/JumpProcesses.jl/}},
  year         = {2022}
}
@misc{zagatti_extending_2023,
	title = {Extending {JumpProcess}.jl for fast point process simulation with time-varying intensities},
	url = {http://arxiv.org/abs/2306.06992},
	doi = {10.48550/arXiv.2306.06992},
	publisher = {arXiv},
	author = {Zagatti, Guilherme Augusto and Isaacson, Samuel A. and Rackauckas, Christopher and Ilin, Vasily and Ng, See-Kiong and Bressan, Stéphane},
	year = {2023},
}