Worked Example of a Two-Branch Lineage Model

From a lineage diagram to a count model

A branching diagram describes a biological hypothesis, but estimation requires a probability model for the observed counts. This tutorial works through a deliberately small, hypothetical example: cells follow a shared trajectory until a split, after which they belong to one of two branches with different expression rates for a gene. The aim is to show how the model links counts to branch assignments and how candidate split points and parameter uncertainty can be examined. The example is simplified; real expression trajectories can change continuously along each branch.

Let \(Y_i\) be the count for a gene in cell \(i\), \(t_i\) its position along an ordered trajectory, and \(\tau\) the split position. Let \(e_i\) be a known exposure, such as a library-size offset on a suitable scale. A latent variable \(Z_i\) records the branch for a post-split cell: \(Z_i=A\) or \(Z_i=B\). Before the split, the model uses a shared rate \(\lambda_0\); after it, the branch-specific rates are \(\lambda_A\) and \(\lambda_B\).

This model assumes that the trajectory positions and exposures are available, the split is represented by a single \(\tau\), and a Poisson distribution is an adequate conditional count model. It also assumes a constant rate within each modeled segment: one shared pre-split rate and one rate per post-split branch. Those assumptions are restrictive, but they make the relationship between the data and the estimates transparent. The earlier material on generative models for tree-structured cell fates and likelihoods for branching trajectories establishes the general latent-branch framework; here the emphasis is on carrying out a complete numerical fit.

For a post-split cell, its marginal count probability is a mixture: \(P(Y_i=y)=\pi f_{\mathrm{Pois}}(y;e_i\lambda_A)+(1-\pi)f_{\mathrm{Pois}}(y;e_i\lambda_B)\). The posterior probability of branch \(B\), sometimes called its responsibility, follows from Bayes' rule:

This is a probability for the cell's branch conditional on the observed count and the current model parameters. It is not the same quantity as the population branch probability \(1-\pi\). A count that is much more likely under branch \(B\) than under branch \(A\) can have a high branch-\(B\) responsibility even when branch \(B\) is less common overall.

Updating the branch probabilities and rates

The expectation-maximization (EM) algorithm alternates between estimating uncertain branch memberships and updating model parameters. At the E-step, compute each post-split cell's responsibility using the current parameters. At the M-step, maximize the expected complete-data log likelihood, treating each cell as fractionally assigned according to its responsibility. For cells with equal exposure \(e_i=1\), the rate update for branch \(k\) is the responsibility-weighted count total divided by the responsibility-weighted cell total.

where \(m\) is the number of post-split cells and \(r_{iA}=1-r_{iB}\). With unequal known exposures, the rate update is \(\lambda_k^{\mathrm{new}}=\sum_i r_{ik}Y_i/\sum_i r_{ik}e_i\). That denominator accounts for the amount of exposure contributing to each branch's counts.

Estimating the split point by profiling

The split point is a model parameter, not a label that can be read directly from a branching diagram. A simple approach is to evaluate a prespecified set of candidate split positions. For each candidate, classify cells before the split as shared-trajectory observations, fit the shared and branch-specific rates under that candidate, and calculate the maximized likelihood. The candidate with the largest profile likelihood is preferred within the candidate set. If branch labels are latent, the mixture likelihood must be used; if labels are observed, the corresponding branch-specific likelihood can be used directly.

The following example uses known hypothetical branch labels to keep the profile calculation explicit. For independent Poisson counts with total count \(Y\) and total exposure \(E\) in a group, the maximum-likelihood rate is \(\widehat{\lambda}=Y/E\). Omitting the factorial terms, which are constant across the split candidates for the same observed cells, that group's maximized log likelihood is \(Y\log(Y/E)-Y\). For \(Y=0\), the term \(Y\log(Y/E)\) is defined as zero.

What the uncertainty statements mean

Parameter uncertainty has more than one source. Even with the split point and branch labels fixed, the rates are estimated from finite counts. In the last example, conditional on \(\tau_1\) and the known labels, the branch-\(A\) estimate is \(2/2=1\) count per exposure unit and the branch-\(B\) estimate is \(11/2=5.5\). A large-sample log-scale approximation gives a standard error for \(\log(\widehat{\lambda})\) of \(1/\sqrt{Y}\), when the Poisson count total is positive. An approximate 95% interval is \(\exp[\log(\widehat{\lambda})\pm1.96/\sqrt{Y}]\).

For branch \(A\), this gives \(\exp[0\pm1.96/\sqrt{2}]\approx[0.250,4.00]\) counts per exposure unit. For branch \(B\), it gives \(\exp[\log(5.5)\pm1.96/\sqrt{11}]\approx[3.05,9.93]\) counts per exposure unit. These are approximate intervals conditional on the selected split and known assignments; with counts this small, the approximation may be rough. They do not account for uncertainty in branch labels or in \(\tau\), which can make overall uncertainty wider.

For latent labels, responsibilities summarize assignment uncertainty under the fitted parameters. For the split point, a profile likelihood across candidate values shows how much support changes with \(\tau\), but choosing a confidence set requires an appropriate calibration. In particular, a split point can be weakly identified when branch-specific rates are similar, few cells lie near the split, or the candidate positions give nearly equal likelihoods. The earlier material on identifiability of branching trajectory models explains why such ambiguity cannot always be resolved by collecting more counts from the same uninformative region.

Diagnostics, alternatives, and limitations

A numerical fit should be examined before its branch assignments are treated as meaningful. Check that estimated rates are finite and nonnegative, that branch proportions lie between zero and one, and that the fitted likelihood is not decreasing across valid EM iterations. Because mixture likelihoods can have local maxima, repeat the fit from multiple initial values. Examine whether responsibilities are nearly all zero or one, whether many are intermediate, and whether the split profile is sharply peaked or nearly flat. Extreme responsibilities are not automatically wrong, but a fit that changes substantially across starting values signals weak identification or optimization difficulty.

The Poisson assumption entails equality of conditional mean and variance. Substantial extra variation, unmodeled library-size differences, or a changing rate within a branch can make the model's assignment probabilities overconfident or its rates misleading. A more flexible count distribution or a time-varying branch trajectory may be needed, but each adds assumptions and parameters. This simple model also treats the supplied trajectory positions and exposure offsets as known; uncertainty in their estimation is not included.

Hard assignment is an alternative: assign every cell to its most probable branch and then estimate rates. It is simple to communicate, but discards uncertainty and can distort estimates when assignments are ambiguous. Soft assignment retains that uncertainty through responsibilities, though it remains conditional on the model and its current parameters. If lineage labels are observed from an independent source, a labeled likelihood can be more direct; if labels are latent, a mixture likelihood is appropriate. Neither approach makes a one-gene, constant-rate model a complete description of cell development. Branch-specific expression association is not evidence that branch identity causes the observed count.