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 DifferentialEquations, 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)
dprob = DiscreteProblem(rn, u0, tspan, p)
jprob = JumpProblem(rn, dprob, Direct())

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.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, DiscreteFunction{true, true, DiscreteFunction{true, SciMLBase.FullSpecialize, SciMLBase.var"#292#293", Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Nothing, ModelingToolkit.ObservedFunctionCache{JumpSystem{ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}, Vector{Equation}}}}}, JumpSystem{ArrayPartition{Any, Tuple{Vector{MassActionJump}, Vector{ConstantRateJump}, Vector{VariableRateJump}, Vector{Equation}}}}, Nothing}, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}}, Direct, CallbackSet{Tuple{}, Tuple{DiscreteCallback{JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{Any}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{Any}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, JumpProcesses.DirectJumpAggregation{Float64, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{Any}, 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}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{Any}, Float64}}, Tuple{}, Tuple{}, TaskLocalRNG}, Tuple{}, Tuple{}, Nothing, MassActionJump{Vector{Float64}, Vector{Vector{Pair{Int64, Int64}}}, Vector{Vector{Pair{Int64, Int64}}}, ModelingToolkit.JumpSysMajParamMapper{Vector{Num}, Vector{Any}, Float64}}, TaskLocalRNG, Base.Pairs{Symbol, Nothing, Nothing, @NamedTuple{callback::Nothing}}}, Exponential{Float64}, Int64}}

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];
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")