Simulation of branching one-dimensional Ornstein-Uhlenbeck processes

Set up the environment

using BranchingProcesses
using StochasticDiffEq
using Distributions
using LaTeXStrings
using Plots
using Random

Define the model parameters

f(u,p,t) = p[2]*(p[1]-u)	# drift function with rate p[2] around mean p[1]
g(u,p,t) = p[3]             # noise function with standard deviation p[3]
tspan_short = (0.0, 5.0)	# short time span for full trajectory sampling
tspan = (0.0, 5.0)			# time span for fluctuation experiments
dt = 0.01
μ = 0.0 					# OUP mean
σst = 1.0 					# OUP stationary state variance
dst = Normal(μ, σst); 		# OUP stationary distribution
λ = 1.0 					# branching rate
np = 2 						# deterministic number of offspring
vp = 2 						# deterministic number of offspring
r = λ*(np-1) 				# population growth rate
α = [r, 0.5r, 0.025r] 		# drift parameters for fast, critical, and slow fluctuations
σ = σst * sqrt.(2α)			# diffusion parameters for fast, critical, and slow fluctuations
nclone = 100 				# number of clones

100

Compute the theoretical (population average) clonal variance

Ratio of clonal variance to variance in the same number of independent cells for the observable $f(x)=x-\mu$:

\[\begin{aligned} \frac{\mathrm{Var}_{\mathrm{int}}(S_t)}{e^{\lambda t}\sigma^2_{\mathrm{st}}} &= 1 + v_p \int_0^tds\, e^{\lambda s} e^{-2\alpha s}\\ &= \begin{cases} 1 + \frac{v_p}{\lambda-2\alpha}(e^{(\lambda-2\alpha)t}-1) & 2\alpha\neq\lambda\\ 1 + v_p t & 2\alpha=\lambda \end{cases} \end{aligned}\]

function varratio(trange, α, λ, vp)
	if λ == 2α
		return (1. .+ vp .* trange);
	elseif λ > 2α
		return (1. .+ vp .* (exp.((λ-2α).*trange) .- 1.) ./ (λ-2α))
	else
		return (1. .+ vp .* (exp.((λ-2α).*trange) .- 1.) ./ (λ-2α))
	end
end;

varratio (generic function with 1 method)

trange = range(tspan[1], tspan[2], 200);
varratio_theo = σst^2 * [varratio(trange, αi, λ, vp) for αi in α]

3-element Vector{AbstractVector{Float64}}:
 [1.0, 1.0496252132929347, 1.0980190956886848, 1.1452121997164537, 1.191234319817556, 1.236114511155553, 1.279881107959662, 1.3225617414130162, 1.3641833570970703, 1.4047722320031653  …  2.983104739939309, 2.9835239553813846, 2.9839327689955963, 2.9843314388784066, 2.9847202167222324, 2.985099347974347, 2.9854690719918358, 2.9858296221927105, 2.986181226203274, 2.986524106001829]
 1.0:0.05025125628140704:11.0
 [1.0, 1.050855788962096, 1.1029400757778969, 1.1562825366555074, 1.2109135646763451, 1.2668642871122924, 1.3241665831611678, 1.382853102110623, 1.4429572819408194, 1.5045133683764789  …  195.189290193753, 199.93107952409161, 204.7874138872364, 209.76106029024425, 214.85485258137768, 220.07169306475427, 225.41455415399963, 230.88648006584924, 236.49058855465856, 242.23007268881614]

Simulations

Single trajectory sampling

Initialize dummy OUP and branching OUP problems:

u0 = rand(dst)
oup = SDEProblem(f, g, u0, tspan_short, (μ, α[1], σ[1]))
boup = ConstantRateBranchingProblem(oup, λ, np)

ConstantRateBranchingProblem{SDEProblem{Float64, Tuple{Float64, Float64}, false, Tuple{Float64, Float64, Float64}, Nothing, SDEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.var"Main".f), typeof(Main.var"Main".g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, typeof(Main.var"Main".g), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, Nothing}, Exponential{Float64}, Int64}(SDEProblem{Float64, Tuple{Float64, Float64}, false, Tuple{Float64, Float64, Float64}, Nothing, SDEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.var"Main".f), typeof(Main.var"Main".g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}, typeof(Main.var"Main".g), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, Nothing}(SDEFunction{false, SciMLBase.AutoSpecialize, typeof(Main.var"Main".f), typeof(Main.var"Main".g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}(Main.var"Main".f, Main.var"Main".g, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing, nothing), Main.var"Main".g, -0.07979991959098581, (0.0, 5.0), (0.0, 1.0, 1.4142135623730951), nothing, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}(), nothing, 0x0000000000000000), Exponential{Float64}(θ=1.0), 2, 0)