Estimating Branchpoint Position and Timing

What a Branchpoint Estimate Represents

A branching trajectory has a shared path before a split and branch-specific paths afterward. Estimating the split means locating both its time and its state: the time is the point on the trajectory at which the paths separate, and the state is the predicted cell state at that time. Because observed cell states are noisy, the estimate is not usually a visible intersection of two perfectly observed curves. It is a parameter value that makes the observed cells most compatible with the assumed branching model.

Here, time \(t_i\) denotes a supplied trajectory coordinate for cell \(i\), and \(\mathbf{y}_i\) is that cell’s measured state in a \(d\)-dimensional state space. The coordinate could be a time measurement or an estimated progression coordinate. Treating it as known is an assumption: uncertainty in that coordinate is not included in the calculations below.

A Likelihood Model for Timing and Position

Consider two branches, labelled \(z_i\in\{A,B\}\), with known labels for the purpose of fitting. Let \(\mathbf{a}\) be the shared trajectory’s starting state, \(\mathbf{v}_0\) its velocity before the split, and \(\mathbf{v}_{z_i}\) the velocity on cell \(i\)’s branch afterward. The continuous piecewise-linear mean state is

For \(t_i\leq\tau\), this gives the shared trajectory \(\mathbf{a}+\mathbf{v}_0t_i\). For \(t_i>\tau\), it gives \(\mathbf{a}+\mathbf{v}_0\tau+\mathbf{v}_{z_i}(t_i-\tau)\). Thus the branch trajectories meet continuously at

A common observation model is \(\mathbf{y}_i\sim N_d(\mathbf{m}_i(\tau),\boldsymbol{\Sigma})\), independently across cells conditional on the model parameters, with a positive-definite measurement covariance matrix \(\boldsymbol{\Sigma}\). This assumes that the mean trajectory has the stated form, the supplied times and branch labels are correct, and the residual distribution and covariance are adequately specified. It also treats cells as independent observations; shared biological sources of variation can make that assumption inappropriate.

If \(\mathbf{a}\), the velocities, and \(\boldsymbol{\Sigma}\) are fixed, the log likelihood for \(\tau\), up to terms that do not depend on \(\tau\), is minus one half the weighted residual sum of squares:

The maximum-likelihood estimate \(\hat{\tau}\) minimizes \(Q(\tau)\) over the allowed time range. The estimated branchpoint is \(\hat{\mathbf{b}}=\mathbf{a}+\mathbf{v}_0\hat{\tau}\). If the trajectory parameters are also unknown, they must be estimated jointly or treated as nuisance parameters when assessing the likelihood for \(\tau\). The calculations in the first example hold them fixed to show how timing estimation works.

Worked Examples

Uncertainty in the Estimated Time

For a fixed set of cells that remain post-split over the candidate range of \(\tau\), the model is linear in \(\tau\). In one dimension, a post-split mean can be written \(a+v_{z_i}t_i+c_i\tau\), where \(c_i=v_0-v_{z_i}\). Define \(r_i=y_i-a-v_{z_i}t_i\). Then \(r_i=c_i\tau+\epsilon_i\). If the error variances \(\sigma_i^2\) are known and errors are independent and Gaussian, maximizing the likelihood is weighted least squares. Differentiating the weighted squared residuals and setting the derivative to zero gives

The curvature of the log likelihood in this case is \(-\sum_i c_i^2/\sigma_i^2\). Consequently, the model-based standard error is

This formula explains why precision depends both on measurement noise and on how differently the branches move from the shared path. It applies directly only while the set of pre- and post-split cells stays fixed and the trajectory parameters are treated as known. In general, evaluate the likelihood across the allowed range, since cells can change from pre-split to post-split as \(\tau\) changes.

A general likelihood-based interval can be formed by retaining values of \(\tau\) whose likelihood is sufficiently close to its maximum. With known covariance and standard regularity conditions, a 95% profile-likelihood confidence set uses values satisfying \(2\{\ell(\hat{\tau})-\ell(\tau)\}\leq3.84\), the approximate 95th percentile of a chi-squared distribution with one degree of freedom. If nuisance trajectory parameters are estimated, the likelihood at each candidate \(\tau\) must be maximized over those parameters before applying the comparison. The approximation can be unreliable near parameter boundaries or when the likelihood is very flat; in such cases, report the profile itself or assess uncertainty by simulating data from the fitted model and refitting it.

Diagnostics, Alternatives, and Limitations

Inspect the profile likelihood or residuals rather than relying only on a point estimate. A clear interior minimum and residuals without systematic patterns support the fitted timing more than a flat profile or a minimum at the edge of the allowed time range. Check whether residual spread changes by branch or time, whether a few cells dominate the fit, and whether the assumed linear paths describe the observed trends. A split estimate that moves substantially when plausible cells are removed or the time range is changed is sensitive, even if a conditional standard error looks small.

The likelihood approach uses the full state-space residuals and can incorporate unequal measurement variances through \(\boldsymbol{\Sigma}\). A simpler alternative is to fit separate branch curves and take their intersection. That can be useful as a visual or exploratory summary, but it may not enforce a shared pre-split path or account for observation noise in a coherent likelihood. Searching a grid of candidate split times and computing \(Q(\tau)\) at each is a practical way to find a minimum; it approximates the same likelihood fit, provided the grid is fine enough and the candidate range is appropriate.

This model’s uncertainty is conditional on its assumptions. Incorrect branch labels, uncertainty in the supplied times, non-linear trajectories, correlated measurements, or a misspecified noise covariance can distort both the estimated split and its interval. If the branch-specific velocities are close to the shared velocity, observations may change little as the proposed split time changes, leaving little information about \(\tau\). Such a result should be reported as weakly determined rather than as a precise biological event. The estimated branchpoint is a feature of the fitted trajectory, not proof that a discrete biological decision occurred at that exact state or time.