Utilities

tip_values(tree::BranchingProcessSolution)

Return the values at the tips of a [BranchingProcessSolution](@docs) tree.

See also: nodevalue

reduce_tree(tree::BranchingProcessSolution)

Reduce BranchingProcessSolution tree to an ordinary time series by combining the values of all particles alive at each time point. The timespan of the resulting time series is from the initial time of the root to the latest final time of any of the leaves of the input tree; the time points are spaced by dt (default: 0.01). The function transform (default: identity) is applied to the values of each particle before combining them. The function reduction (default: sum) is used to combine the (transformed) values of all particles alive at each time point. It can be any callable that accepts either a vector of particle values (reduction(vals)) or both the current time step and a vector of particle values (reduction(t, vals)), and returns a single aggregated value. The two-argument form allows time-dependent rescaling to be applied during reduction. For backward compatibility, the strings "sum" and "prod" are also accepted.

Note that if the resulting time series has different time points than the original trajectories in the input tree, the interpolating function of the solver used to sample the original trajectories is used to compute the reduced time series. Any keyword arguments kwargs... (for instance idxs=[1,3,5] to summarize only a subset variables) are passed to the interpolating function.

rescale(sol::ReducedBranchingProcessSolution, f)

Rescale the solution values of a ReducedBranchingProcessSolution elementwise by a function f of the corresponding time step values. Returns a new ReducedBranchingProcessSolution with values f(t) .* u at each time point t.

The scaling function f should accept a single time value and return either a scalar or a vector compatible with the solution values u. This is useful, for example, to normalize the particle sum by the expected number of particles in order to study fluctuations around the mean field.

The transform field of the returned solution is set to f ∘ sol.transform to record the full chain of transformations applied.

Examples

# Normalize by the square root of the expected number of particles exp(lambda*t) to study fluctuations
lambda = 1.0
rescaled = rescale(sol, t -> exp(-0.5 * lambda * t))

See also: rescale!, reduce_tree, ReducedBranchingProcessSolution

rescale!(sol::ReducedBranchingProcessSolution, f)

Rescale the solution values of a ReducedBranchingProcessSolution elementwise in-place by a function f of the corresponding time step values. Modifies sol.u directly so that each element becomes f(t) .* u at the corresponding time point t. Returns the modified solution.

The scaling function f should accept a single time value and return either a scalar or a vector compatible with the solution values u.

Note: because ReducedBranchingProcessSolution is an immutable struct, only the mutable u field is updated in-place; the transform field is not changed.

Examples

lambda = 1.0
rescale!(sol, t -> exp(-lambda * t))

See also: rescale, reduce_tree, ReducedBranchingProcessSolution

node_generations(root::BranchingProcessNode)

Compute the generation (distance from root) for all nodes in the tree. Returns a dictionary mapping each node to its generation, where the root has generation 0, its direct children have generation 1, etc.

fluctuation_experiment(bp::ConstantRateBranchingProblem, u0_dist, nclone::Integer; reduction=sum, ensemble_alg=EnsembleThreads(), alg=nothing, solver_kwargs=NamedTuple(), reduce_kwargs=NamedTuple(), rescale=nothing)

Simulate a Luria-Delbrück fluctuation experiment by running nclone independent branching process simulations, each starting from a different initial state sampled from u0_dist.

For each clone, the initial condition of the inner problem in bp is replaced with a sample drawn from u0_dist using SciMLBase.remake. The simulations are executed using the SciML EnsembleProblem interface, and each clone's branching tree is reduced to a time series via reduce_tree.

Arguments

• bp: A ConstantRateBranchingProblem defining the single-particle dynamics, lifetime distribution, and offspring distribution.
• u0_dist: A distribution (from Distributions.jl) from which the initial state for each clone is independently sampled. Must be compatible with the initial-condition type of bp.prob (e.g. a UnivariateDistribution for scalar problems, a MultivariateDistribution for vector-valued problems). Alternatively, a deterministic initial state of the same type as the initial state of bp.prob, in which case all clones start from this same state.
• nclone: The number of independent clones (trajectories) to simulate.

Keyword Arguments

• reduction=sum: The reduction function used in reduce_tree to aggregate particle values at each time point. Any callable that accepts a vector of particle values and returns a single aggregated value is accepted; defaults to sum.
• ensemble_alg=EnsembleThreads(): The SciML ensemble algorithm controlling how trajectories are parallelised. Common choices are EnsembleSerial(), EnsembleThreads() (default), and EnsembleDistributed().
• alg=nothing: The algorithm passed to solve for each individual trajectory. Defaults to nothing (automatic algorithm selection).
• solver_kwargs=NamedTuple(): Additional keyword arguments passed to the solve function for each individual trajectory (e.g. (; dt=0.1, saveat=0:0.1:5, reltol=1e-6)).
• reduce_kwargs=NamedTuple(): Additional keyword arguments passed to solve_and_reduce (e.g. (; output_dt=0.01, transform=log)).
• rescale=nothing: An optional scaling function applied to each clone's reduced solution via rescale! before returning. The function should accept a single time value and return a scalar or vector compatible with the solution values.

Returns

An EnsembleSolution whose elements are ReducedBranchingProcessSolution objects, one per clone.

Examples

using BranchingProcesses, Distributions, StochasticDiffEq

# Define a simple branching diffusion (Brownian motion with branching rate 1)
f(u, p, t) = 0.0
g(u, p, t) = 0.5
prob = SDEProblem(f, g, 1.0, (0.0, 5.0))
bp = ConstantRateBranchingProblem(prob, 1.0, 2)

# Simulate 100 clones with initial states drawn from a log-normal distribution
results = fluctuation_experiment(bp, LogNormal(0.0, 0.5), 100)

# Use different dt for solver and solve_and_reduce
results = fluctuation_experiment(bp, LogNormal(0.0, 0.5), 100;
                                solver_kwargs=(; dt=0.1, saveat=0:0.1:5),
                                reduce_kwargs=(; dt=0.02, transform=log))

# Use a custom reduction function with specific solver tolerances
results_max = fluctuation_experiment(bp, LogNormal(0.0, 0.5), 100; 
                                   reduction=maximum,
                                   solver_kwargs=(; reltol=1e-8, abstol=1e-10))

# Rescale each clone's solution by exp(-lambda*t) before returning
lambda = 1.0
results_rescaled = fluctuation_experiment(bp, LogNormal(0.0, 0.5), 100;
                                         rescale=t -> exp(-lambda * t))

See also: ConstantRateBranchingProblem, solve_and_reduce, rescale!, ReducedBranchingProcessSolution

timestep_crosscov(sim, i)

Compute the cross-covariance matrix between variables at time step i of an ensemble solution. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects (i.e. the .u field of an EnsembleSolution). Given N ensemble members each with a d-dimensional state vector, this function returns the d×d sample covariance matrix whose (j,k) entry is the sample covariance between variable j and variable k across all trajectories at time step i.

Returns the covariance matrix as a vector of length d² (using C[:]).

See also: timeseries_steps_crosscov, timestep_crosscor

timeseries_steps_crosscov(sim)

Compute the cross-covariance matrix between variables at each time step of an ensemble solution. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects (i.e. the .u field of an EnsembleSolution). For each time step i, returns the d×d sample covariance matrix whose (j,k) entry is the sample covariance between variable j and variable k across all trajectories.

Returns a DiffEqArray containing the covariance vectors (each of length d²) at each time step, together with the corresponding time values.

See also: timestep_crosscov, timeseries_steps_crosscor

timestep_crosscor(sim, i)

Compute the cross-correlation matrix between variables at time step i of an ensemble solution. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects (i.e. the .u field of an EnsembleSolution). Given N ensemble members each with a d-dimensional state vector, this function returns the d×d sample correlation matrix whose (j,k) entry is the sample correlation between variable j and variable k across all trajectories at time step i.

Returns the correlation matrix as a vector of length d² (using R[:]). Entries where a variable has zero variance are set to NaN.

See also: timeseries_steps_crosscor, timestep_crosscov

timeseries_steps_crosscor(sim)

Compute the cross-correlation matrix between variables at each time step of an ensemble solution. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects (i.e. the .u field of an EnsembleSolution). For each time step i, returns the d×d sample correlation matrix whose (j,k) entry is the sample correlation between variable j and variable k across all trajectories.

Returns a DiffEqArray containing the correlation vectors (each of length d²) at each time step, together with the corresponding time values. Entries where a variable has zero variance are set to NaN.

See also: timestep_crosscor, timeseries_steps_crosscov

timeseries_steps_crosscov_variance_explained_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of percent variance explained by the eigenvectors of the cross-covariance matrix at each time step. The percent variance explained is computed as λᵢ / tr(C), where λᵢ is an eigenvalue of the cross-covariance matrix C. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects.

The returned BootstrappedTimeSeriesSolution stores the bootstrapped expected percent variance explained in u and the corresponding lower and upper confidence bounds in lower and upper.

timestep_particle_number(sim, i)

Return the particle numbers across clones at time step i.

timeseries_steps_particle_number(sim)

Return the particle numbers across clones at each time step.

timeseries_steps_particle_number_mean(sim)

Return the time series of mean particle number across clones.

timeseries_steps_particle_number_var(sim)

Return the time series of particle-number variance across clones.

timestep_mean(sim, i)

Compute clonal means at time step i, defined as the cross-clone mean of the reduced state divided by the particle number.

timeseries_steps_mean(sim)

Compute clonal means at each time step.

timestep_intrinsic_crosscov(sim, i)

Compute intrinsic clonal cross-covariance at time step i by subtracting μμ' * var(N) from the clonal cross-covariance, where μ is the clonal mean vector and N the particle number across clones.

timeseries_steps_intrinsic_crosscov(sim)

Compute intrinsic clonal cross-covariance at each time step.

timestep_intrinsic_var(sim, i)

Compute intrinsic clonal variances at time step i.

timeseries_steps_intrinsic_var(sim)

Compute intrinsic clonal variances at each time step.

timeseries_steps_particle_number_mean_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of mean particle number across clones.

timeseries_steps_particle_number_var_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of particle-number variance across clones.

timeseries_steps_mean_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of clonal means.

timeseries_steps_intrinsic_crosscov_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of intrinsic clonal cross-covariance vectors.

timeseries_steps_intrinsic_var_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of intrinsic clonal variances.

timeseries_steps_crosscov_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of cross-covariance vectors returned by timeseries_steps_crosscov. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects.

The returned BootstrappedTimeSeriesSolution stores the bootstrapped expected cross-covariances in u and the corresponding lower and upper confidence bounds in lower and upper.

timeseries_steps_crosscor_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of cross-correlation vectors returned by timeseries_steps_crosscor. sim can be either an EnsembleSolution or a vector of ReducedBranchingProcessSolution objects.

The returned BootstrappedTimeSeriesSolution stores the bootstrapped expected cross-correlations in u and the corresponding lower and upper confidence bounds in lower and upper.

timeseries_steps_intrinsic_crosscor(sim)

Compute intrinsic clonal cross-correlations at each time step.

timeseries_steps_intrinsic_crosscor_bootstrap(sim; sampling=BasicSampling(1000), confint_method=BasicConfInt, level=0.95)

Bootstrap the time series of intrinsic clonal cross-correlation vectors.

timestep_intrinsic_crosscor(sim, i)

Compute intrinsic clonal cross-correlations at time step i. The intrinsic cross-correlation matrix is obtained by normalizing the intrinsic cross-covariance matrix by the product of intrinsic standard deviations. Entries where the intrinsic variance of a variable is non-positive are set to NaN.

get_transform(sol::ReducedBranchingProcessSolution)

Get the transformation function used in the reduction.

get_reduction_method(sol::ReducedBranchingProcessSolution)

Get the reduction method used to combine particle values.

get_original_solution(sol::ReducedBranchingProcessSolution)

Get the original BranchingProcessSolution that was reduced (if available).

has_original_solution(sol::ReducedBranchingProcessSolution)

Check if the original BranchingProcessSolution is stored.