Branching Ornstein-Uhlenbeck process

The Ornstein-Uhlenbeck process (OUP) is defined by the stochastic differential equation

\[dX_t = \alpha (\mu - X_t)dt + \sigma dW_t\]

and describes a variable $X_t$ which fluctuates stochastically around its average or steady-state value $\mu$; $\alpha>0$ is the rate at which fluctuations return to their steady-state value, $\sigma>0$ defines the relative magnitude of the stochastic fluctuations, and $W_t$ is the Wiener process.

If $X_t$ represents the expression level of a gene, then the OUP is the simplest model of stochastic gene expression. The branching OUP can then be used to model the evolution of gene expression during cell proliferation.^[1] At the opposite end of biological time scales, the branching OUP is also used as a model for the evolution of traits under selection in a phylogenetic tree.^[2]

Setting up the problem

Similar to the branching Brownian motion tutorial, we cannot (yet) use the distributionally exact OUP implementation from the DiffEqNoiseProcess package, but must define the OUP as a SDEProblem:

using DifferentialEquations
f(u,p,t) = p[2]*(p[1]-u)
g(u,p,t) = p[3]
u0 = 0.0
tspan = (0.0, 5.0)
μ = 2.0
α = 5.0
σ = 0.5
oup = SDEProblem(f,g, u0, tspan, (μ, α, σ))

SDEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 5.0)
u0: 0.0

We can first verify that after an initial "burn-in" period, the OUP indeed fluctuates around its steady-state value $\mu$:

using Plots
sol = solve(oup, EM(), dt=0.01)
plot(sol; linewidth=2, legend=false)

We can now set up the branching OUP problem in the same way as in the branching Brownian motion tutorial:

using BranchingProcesses
λ = 1.0         # branching rate
nchild = 2      # deterministic number of offspring
boup = ConstantRateBranchingProblem(oup, λ, nchild)

ConstantRateBranchingProblem{SDEProblem{Float64, Tuple{Float64, Float64}, false, Tuple{Float64, Float64, Float64}, Nothing, SDEFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), typeof(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.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.FullSpecialize, typeof(Main.f), typeof(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.g), Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}, Nothing}(SDEFunction{false, SciMLBase.FullSpecialize, typeof(Main.f), typeof(Main.g), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing, Nothing}(Main.f, Main.g, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, SciMLBase.DEFAULT_OBSERVED, nothing, nothing, nothing), Main.g, 0.0, (0.0, 5.0), (2.0, 5.0, 0.5), nothing, Base.Pairs{Symbol, Union{}, Nothing, @NamedTuple{}}(), nothing, 0x0000000000000000), Exponential{Float64}(θ=1.0), 2)

Sampling a trajectory

To sample a tranjectory of the branching process, we call the solve function, resulting in a BranchingProcessSolution tree, that can be visualized using the standard plot function:

using Random
Random.seed!(123)
sol = solve(boup, EM(); dt=0.01);
plot(sol; linewidth=2, branchpoints=true)

We observe that after a few generations, all cells fluctuate around their steady state value $\mu=2.0$.

Memory in the branching Ornstein-Uhlenbeck process

The dynamics of the branching OUP is driven by two parallel processes, the return to steady-state with rate $\alpha$ and the branching with rate $\lambda$. If the half-life $\ln 2/\alpha$ of the stochastic fluctuations is short compared to the average lifetime $1/\lambda$ of a cell, each cell has sufficient time to equilibrate and the population as a whole will be distributed around the steady-state value, as in the figure above.

However, if the half-life of the stochastic fluctuations is long compared to the average lifetime of a cell, large fluctuations early in the expansion can have a long-lasting effect on the population as a whole. This is the memory phenomenon referred to in the BranchingProcesses package introduction. For the branching OUP specifically, this phenomenon has been analyzed in great mathematical detail.

To illustrate this phenomenon, we set the same random seed as above (to ensure exactly the same branching times) and only change the value of $\alpha$:

α = 0.5
oup = SDEProblem(f,g, u0, tspan, (μ, α, σ))
boup = ConstantRateBranchingProblem(oup, λ, nchild)
Random.seed!(123)
sol = solve(boup, EM(); dt=0.01)
plot(sol; linewidth=2, branchpoints=true)

We observe that the system now did not have sufficient time to equilibrate within the typical life-time of a cell, and the population as a whole is shifted towards, that is, remembers, its initial state.

That the memory is due to the slow return to steady-state of the fluctuations and not the value of the initial state can be seen by starting the same problem in the steady state itself:

α = 0.5
u0 = μ
oup = SDEProblem(f,g, u0, tspan, (μ, α, σ))
boup = ConstantRateBranchingProblem(oup, λ, nchild)
Random.seed!(123)
sol = solve(boup, EM(); dt=0.01)
plot(sol; linewidth=2, branchpoints=true)

Again we observe a shifted distribution due to large fluctuations early in the expansion.