Branching birth-death processes

In the branching Ornstein-Uhlenbeck process example, we considered the Ornstein-Uhlenbeck process (OUP) as a simple model for stochastically fluctuating gene expression. However, gene expression is modelled more realistically as a process taking discrete values representing molecule counts.

If we assume molecules are produced and degraded one at a time, we obtain a birth-death process. Such processes, and many other more complicated processes, are most easily defined using the Catalyst reaction network modelling language.

Setting up the problem

We first define the single-particle, in our case single-cell, dynamics as a reaction network:

using Catalyst
rn = @reaction_network begin
    kp, 0 --> X
    kd, X --> 0
end

\[ \begin{align*} \varnothing &\xrightleftharpoons[\mathtt{kd}]{\mathtt{kp}} \mathrm{X} \end{align*} \]

where kp and kd are the parameters of the model, respectively the production and degradation rates.

As explained in the Catalyst tutorials, we can define a JumpProcess to simulate this reaction network as follows:

using DifferentialEquations, JumpProcesses
u0 = [200]
tspan = (0.0, 3.0)
p = [:kp => 50.0, :kd => 0.25]
jinput = JumpInputs(rn, u0, tspan, p)
jprob = JumpProblem(jinput)
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: 2

A trajectory for a single cell can be sampled and plotted:

using Plots
jsol = solve(jprob, SSAStepper())
plot(jsol)
Example block output

A branching jump problem is set up as in the branching Brownian motion and branching Ornstein-Uhlenbeck process examples:

using BranchingProcesses
λ = 1.0         # branching rate
nchild = 2      # deterministic number of offspring
bjprob = ConstantRateBranchingProblem(jprob, λ, nchild);

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:

sol = solve(bjprob,  SSAStepper());
plot(sol; linewidth=2, branchpoints=true)
Example block output