Ensemble/Monte Carlo Simulations
In many contexts, a single model is re-simulated under similar conditions. Examples include:
- Performing Monte Carlo simulations of a stochastic model to gain insight in its behaviour.
- Scanning a model's behaviour for different parameter values and/or initial conditions.
While this can be handled using for
loops, it is typically better to first create an EnsembleProblem
, and then perform an ensemble simulation. Advantages include a more concise interface and the option for automatic simulation parallelisation. Here we provide a short tutorial on how to perform parallel ensemble simulations, with a more extensive documentation being available here.
Monte Carlo simulations using unmodified conditions
We will first consider Monte Carlo simulations where the simulation conditions are identical in-between simulations. First, we declare a simple self-activation loop model
using Catalyst
sa_model = @reaction_network begin
v0 + hill(X,v,K,n), ∅ --> X
deg, X --> ∅
end
u0 = [:X => 10.0]
tspan = (0.0, 1000.0)
ps = [:v0 => 0.1, :v => 2.5, :K => 40.0, :n => 4.0, :deg => 0.01]
We wish to simulate it as an SDE. Rather than performing a single simulation, however, we want to perform multiple ones. Here, we first create a normal SDEProblem
, and use it as the single input to a EnsembleProblem
(EnsembleProblem
are created similarly for ODE and jump simulations, but the ODEProblem
or JumpProblem
is used instead).
using StochasticDiffEq
sprob = SDEProblem(sa_model, u0, tspan, ps)
eprob = EnsembleProblem(sprob)
Next, the EnsembleProblem
can be used as input to the solve
command. Here, we use exactly the same inputs that we use for single simulations, however, we add a trajectories
argument to denote how many simulations we wish to carry out. Here we perform 10 simulations:
sols = solve(eprob, STrapezoid(); trajectories = 10)
Finally, we can use our ensemble simulation solution as input to plot
(just like normal simulations):
using Plots
plot(sols)
Here, each simulation is displayed as an individual trajectory.
While not used here, the la
plotting option (which modifies line transparency) can help improve the plot visual when a large number of (overlapping) lines are plotted.
Various convenience functions are available for analysing and plotting ensemble simulations (a full list can be found here). Here, we use these to first create an EnsembleSummary
(retrieving each simulation's value at time points 0.0, 1.0, 2.0, ... 1000.0
). Next, we use this as an input to the plot
command, which automatically plots the mean $X$ activity across the ensemble, while also displaying the 5% and 95% quantiles as the shaded area:
e_sumary = EnsembleAnalysis.EnsembleSummary(sols, 0.0:1.0:1000.0)
plot(e_sumary)
Ensemble simulations using varying simulation conditions
Previously, we assumed that each simulation used the same initial conditions and parameter values. If this is not the case (when e.g. performing a parameter scan), a prob_func
optional argument must be supplied to the EnsembleProblem
, this describes how the problem should be modified for each individual simulation of the ensemble.
Here, we first create an ODEProblem
of our previous self-activation loop:
using OrdinaryDiffEqTsit5
oprob = ODEProblem(sa_model, u0, tspan, ps)
Next, we wish to simulate the model for a range of initial conditions of $X$`. To do this we create a problem function, which takes the following arguments:
prob
: The problem given to ourEnsembleProblem
(which is the problem thatprob_func
modifies in each iteration).i
: The number of this specific Monte Carlo iteration in the interval1:trajectories
.repeat
: The iteration of the repeat of the simulation. Typically1
, but potentially higher if the simulation re-running option is used.
Here we will use the following problem function (utilising remake), which will provide a uniform range of initial concentrations of $X$:
function prob_func(prob, i, repeat)
remake(prob; u0 = [:X => i * 5.0])
end
Next, we create our EnsembleProblem
, and simulate it 10 times:
eprob = EnsembleProblem(oprob; prob_func)
sols = solve(eprob, Tsit5(); trajectories = 10)
plot(sols)
Here we see that the deterministic model (unlike the stochastic one), only activates for some initial conditions (while other tends to an inactive state).
The EnsembleProblem
constructor accepts a few additional optional options (output_func
, reduction
, u_init
, and safetycopy
), which are described in more detail here. These can be used to e.g. rerun an individual simulation which does not fulfil some condition, or extract and save only some selected information from a simulation (rather than the full simulation).