Multi-variable branching stochastic processes

The single-particle or single-cell stochastic process undergoing branching can be a multi-variable process. To illustrate this, consider the example of Michaelis-Menten enzyme kinetics, where an enzyme ($E$) transforms a substrate ($S$) into a product ($P$). The reaction network is defined in the Catalyst library of basic chemical reaction network models:

using Catalyst
mm_system = @reaction_network begin
    kB, S + E --> SE
    kD, SE --> S + E
    kP, SE --> P + E
end

\[ \begin{align*} \mathrm{S} + \mathrm{E} &\xrightleftharpoons[\mathtt{kD}]{\mathtt{kB}} \mathrm{\mathtt{SE}} \\ \mathrm{\mathtt{SE}} &\xrightarrow{\mathtt{kP}} \mathrm{P} + \mathrm{E} \end{align*} \]

We set up a single-cell process (with four variables, $S$. $E$. $SE$. and $P$) as in the model library:

using DifferentialEquations, JumpProcesses
u0 = [:S => 30, :E => 10, :SE => 0, :P => 0]
tspan = (0., 100.)
ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]

jinput = JumpInputs(mm_system, u0, tspan, ps)
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: 3

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 usual:

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

To sample a tranjectory of the branching process, we call the solve function as usual:

using AbstractTrees
tree = solve(bjprob,  SSAStepper());
treeheight(tree)
11

If we plot the solution, by default the trajectory for the first variable is shown:

plot(tree; branchpoints=true)
Example block output

Because the plot recipe piggybacks on the plot recipe for differential equation solutions, variables for plotting can be chose in the same way:

plot(tree; branchpoints=true, idxs=[2])
Example block output

It is possible to plot multiple variables in the same plot, but since branching can result in many particles and line color is already used to distinguish different particles in the tree, this is not particularly illuminating:

plot(tree; branchpoints=true, idxs=[1,2])
Example block output