Tree reduction

In experiments it is often the case that a single measurement (per variable) is obtained for a clone that is grown for set time or number of generations from a single cell. For instance in the MemorySeq paper, total RNA counts were obtained for each clone. In the BranchingProcesses package this can be easily simulated by applying a function to (for instance summing) the Values at the tips of a branching process. However in simulations we may want to do a bit more and obtain the entire time series of a summary statistic from a BranchingProcessSolution. This can be done with the reduce_tree function.

Setting up the problem and sampling a tree

Start by creating a branching birth-death process. First define the single-cell process,

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

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

and then the branching process,

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

Sample and plot a tree:

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

Obtaining a reduced time series

By default, reduce_tree sums the values of all cells alive at a given time, starting from the initial time of the root cell and stopping at the final time of the last living cell, with a time step of dt:

sol_red = reduce_tree(sol; dt=0.01);

The reduced time series has the type ReducedBranchingProcessSolution and a plot recipe is available for this type:

plot(sol_red)
Example block output

It is possible to replace the default summing of values to taking a product, although use cases for this in practice may be rather limited:

reduce_tree(sol; dt=0.01, reduction="prod");

Applying transformations

We can apply a transformation to the values of the process before summing. For instance, to sum the squared values:

sol_red = reduce_tree(sol; dt=0.01, transform=(x -> x.^2));
plot(sol_red)
Example block output

or to count the number of cells alive at any given time:

sol_red = reduce_tree(sol; dt=0.01, transform=(x -> 1));
plot(sol_red)
Example block output

Reducing and transforming multivariable processes

The default for reducing multivariable processes is to sum all variables separately. Let's recreate the multivariable branching process example:

mm_system = @reaction_network begin
    kB, S + E --> SE
    kD, SE --> S + E
    kP, SE --> P + E
end
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)
λ = 0.05
nchild = 2
bjprob = ConstantRateBranchingProblem(jprob, λ, nchild);
sol = solve(bjprob,  SSAStepper());

Obtain a reduced time series for the sum of all variables:

var_names = string.(unknowns(mm_system))
sol_red = reduce_tree(sol; dt=0.1);
plot(sol_red)
Example block output

To summarize only one or a subset of variables, we can either use an idxs keyword argument:

subs = [1,4]
sol_red = reduce_tree(sol; dt=0.1, idxs=subs);
plot(sol_red)
Example block output

or use a transformation:

sol_red = reduce_tree(sol; dt=0.1, transform=(x -> x[subs]));
plot(sol_red);
Example block output

Transformations can also be used to create summaries of summaries, for instance, for the total number of molecules:

sol_red = reduce_tree(sol; dt=0.1, transform=(x -> sum(x)));
plot(sol_red)
Example block output

Tree reduction and ensemble simulations

Tree reduction is particularly useful in the context of ensemble simulations.

ensemble_bjprob = EnsembleProblem(bjprob)
ensemble_tree = solve(ensemble_bjprob, SSAStepper(), EnsembleThreads(), trajectories=100)