Advanced Simulation Options

Throughout the preceding tutorials, we have shown the basics of how to solve ODE, SDE, and jump process models generated from Catalyst ReactionSystems. In this tutorial we'll illustrate some more advanced functionality that can be useful in many modeling contexts, and that provide conveniences for common workflows. For a comprehensive overview of solver properties, parameters, and manipulating solution objects, please read the documentation of the DifferentialEquations package, which Catalyst uses for all simulations.

Monte Carlo simulations using EnsembleProblems

In many contexts one needs to run multiple simulations of a model, for example to collect statistics of SDE or jump process solutions, or to systematically vary parameter values within a model. While it is always possible to manually run such ensembles of simulations via a for loop, DifferentialEquations.jl provides the EnsembleProblem as a convenience to manage structured collections of simulations. EnsembleProblems provide a simple interface for modifying a problem between individual simulations, and offers several options for batching and/or parallelizing simulation runs. For a more thorough description, please read the Parallel Ensemble Simulations section of the DifferentialEquations documentation. Here, we will give a brief introduction to the use of EnsembleProblems from Catalyst-generated models.

Let's look at a single-component bistable self-activation model:

using Catalyst, DifferentialEquations, Plots

rn = @reaction_network begin
    v0 + hill(X,v,K,n), ∅ --> X
    deg, X --> ∅
end
u0 = [:X => 0.0]
tspan = (0.0,1000.0)
p = [:v0 => 0.1, :v => 2.5, :K => 75.0, :n => 2.0, :deg => 0.01];
sprob = SDEProblem(rn, u0, tspan, p)

we can then use our SDEProblem as input to an EnsembleProblem:

eprob = EnsembleProblem(sprob)
EnsembleProblem with problem SDEProblem

The EnsembleProblem can now be used as input to the solve command. It has the same options as when simulating the SDEProblem directly, however, it has an additional argument trajectories to determine how many simulations should be performed.

esol = solve(eprob; trajectories=5)
EnsembleSolution Solution of length 5 with uType:
SciMLBase.RODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, DiffEqNoiseProcess.NoiseProcess{Float64, 2, Float64, Vector{Float64}, Nothing, Nothing, typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_DIST), typeof(DiffEqNoiseProcess.INPLACE_WHITE_NOISE_BRIDGE), true, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, ResettableStacks.ResettableStack{Tuple{Float64, Vector{Float64}, Nothing}, true}, DiffEqNoiseProcess.RSWM{Float64}, Nothing, RandomNumbers.Xorshifts.Xoroshiro128Plus}, SDEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, Nothing, SDEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.var"#f#451"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9f3bcf4d, 0x6f529616, 0x4e1c670d, 0x1896f613, 0x7d2685a4), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x511bdce0, 0xdef82bac, 0xe6d90fc9, 0xfe8f9322, 0x8a889647), Nothing}}, ModelingToolkit.var"#g#452"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb3d4032, 0xb05b474d, 0xe65682e9, 0x10c1823b, 0xd42a14e7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9941be3a, 0x1fda6b39, 0xd8a6d183, 0x31c8b745, 0x3d888d98), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Symbol}, ModelingToolkit.var"#generated_observed#461"{Bool, SDESystem, Dict{Any, Any}}, Nothing, SDESystem}, ModelingToolkit.var"#g#452"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xeb3d4032, 0xb05b474d, 0xe65682e9, 0x10c1823b, 0xd42a14e7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9941be3a, 0x1fda6b39, 0xd8a6d183, 0x31c8b745, 0x3d888d98), Nothing}}, Base.Pairs{Symbol, Nothing, Tuple{Symbol}, NamedTuple{(:callback,), Tuple{Nothing}}}, Matrix{Float64}}, StochasticDiffEq.LambaEM{true}, StochasticDiffEq.LinearInterpolationData{Vector{Vector{Float64}}, Vector{Float64}}, SciMLBase.DEStats, Nothing}

This simulation is automatically multithreaded over all available threads. Please read this documentation for more information on parallelisation alternatives. The ensemble simulations can be plotted using the plot function, which by default displays all trajectories:

plot(esol)

Sometimes when performing a large number of ensemble simulations, the plots get very dense. In these cases, the plot argument linealpha (which sets trajectory transparency) may be useful:

esol = solve(eprob; trajectories = 100)
plot(esol; linealpha = 0.5)

Sometimes, one wishes to perform the same simulation a large number of times, while making minor modifications to the problem each time. This can be done by giving a problem function, prob_func, argument to the EnsembleProblem. Let us consider ODE simulations of a simple birth/death process:

rn = @reaction_network begin
    (b,1.0), ∅ <--> X
end
u0 = [:X => 1.0]
tspan = (0.0, 1.0)
p = [:b => 1.];
oprob = ODEProblem(rn, u0, tspan, p)

We wish to simulate this model for a large number of values of b. We do this by creating a prob_func that will make a modification to the problem at the start of each Monte Carlo simulation:

b_values = 1.0:0.1:2.0
function prob_func(prob, i, repeat)
    @unpack b = prob.f.sys    # Fetches the b parameter to be used in the local scope.
    remake(prob; p = [b => b_values[i]])
end

Here, prob_func takes three arguments:

  • prob: The problem given to our EnsembleProblem, this is the problem that prob_func modifies in each iteration.
  • i: The number of this specific Monte Carlo iteration in the interval 1:trajectories.
  • repeat: The repeat of this specific Monte Carlo simulation (We will ignore

this argument in this brief overview). In our case, for each Monte Carlo simulation, our prob_func takes our original ODEProblem and uses the remake function to change the parameter vector. Here, for the ith Monte Carlo simulation, the value of b is also the ith value of our b_values vector. Finally, we can simulate and plot our problem:

eprob = EnsembleProblem(oprob; prob_func = prob_func)
esol = solve(eprob; trajectories = length(b_values))
plot(esol)

Note that plot legends are disabled when plotting ensemble solutions. These can be re-enabled using the legend plotting keyword. However, when plotting a large number of trajectories, each will generate a label. Sometimes the best approach is to remove these and add a label manually:

p = plot(esol; label = nothing)
plot!(p, Float64[], Float64[]; label = "X", legend = :topleft)