Ensemble simulations

In a Luria-Delbrück experiment, here called a fluctuation experiment, multiple clones are grown, each starting from a different single cell.

To simulate fluctuation experiments, the BranchingProcesses package supports SciML's parallel ensemble simulation methods.^[1]

Let's model the actual Luria-Delbrück experiment as a birth-death process where a cell can switch from wild-type to mutant, without any back-mutations:

using StochasticDiffEq, JumpProcesses, Catalyst
rn = @reaction_network begin
    μ, W --> M
end
u0 = [:W => 1, :M => 0]  # Always start in the wild-type state
p = [:μ => 0.01]  # Mutation rate
tspan = (0.0, 10.0)
jprob = JumpProblem(rn, u0, tspan, p)

JumpProblem with problem DiscreteProblem with aggregator Direct
Number of jumps with discrete aggregation: 0
Number of jumps with continuous aggregation: 0
Number of mass action jumps: 1

Define a branching process problem:

using BranchingProcesses
λ = 1.0
nchild = 2
bjprob = ConstantRateBranchingProblem(jprob, λ, nchild);

To simulate a fluctuation experiment with 100 clones, first set up an EnsembleProblem:

ensemble_bjprob = EnsembleProblem(bjprob)

EnsembleProblem with problem ConstantRateBranchingProblem

Then solve the problem:

ensemble_sol = solve(ensemble_bjprob, SSAStepper(), EnsembleThreads(), trajectories=100)

EnsembleSolution Solution of length 100 with uType:
BranchingProcessSolution{ConstantRateBranchingProblem{JumpProblem{true, DiscreteProblem{Vector{Int64}, Tuple{Float64, Float64}, true, MTKParameters{Vector{Float64}, StaticArraysCore.SVector{0, Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#321#322", Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Nothing, ModelingToolkitBase.ObservedFunctionCache{System, Nothing}, System, Nothing}, Base.Pairs{Symbol, DataType, Nothing, @NamedTuple{u0_eltype::DataType}}}, Direct, CallbackSet{Tuple{}, Tuple{DiscreteCallback{JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkitBase.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkitBase.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkitBase.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, typeof(SciMLBase.FINALIZE_DEFAULT), Nothing, Tuple{}}}}, JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkitBase.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, Tuple{}, Tuple{}, Nothing, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkitBase.JumpSysMajParamMapper{Vector{Num}, Vector{SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}, Float64}}, TaskLocalRNG, Base.Pairs{Symbol, Nothing, Nothing, @NamedTuple{callback::Nothing}}}, Exponential{Float64}, Int64}, ODESolution{Int64, 2, Vector{Vector{Int64}}, Nothing, Nothing, Vector{Float64}, Nothing, Nothing, DiscreteProblem{Vector{Int64}, Tuple{Float64, Float64}, true, MTKParameters{Vector{Float64}, StaticArraysCore.SVector{0, Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.AutoSpecialize, SciMLBase.var"#321#322", Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Nothing, ModelingToolkitBase.ObservedFunctionCache{System, Nothing}, System, Nothing}, Base.Pairs{Symbol, DataType, Nothing, @NamedTuple{u0_eltype::DataType}}}, SSAStepper, SciMLBase.ConstantInterpolation{Vector{Float64}, Vector{Vector{Int64}}}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing, Nothing, SciMLBase.ReturnCode.T}}

We can obtain the number of cells and the number of mutants in each clone using the tip_values function:

cell_counts = [sum(tip_values(sol)) for sol in ensemble_sol.u];
total_cell_counts = [sum(x) for x in cell_counts];
#mutant_cell_counts = [x[2] for x in cell_counts];

using Plots
histogram(total_cell_counts,label="",xlabel="Total cell counts", ylabel="Number of clones")

#histogram(mutant_cell_counts,label="",xlabel="Mutant cell counts", ylabel="Number of clones")