Diffusion Pseudotime from a Root Cell

From Diffusion States to Root-Based Pseudotime

A diffusion operator describes how probability mass moves among cells on a graph. At a given number of steps, each starting cell has a probability profile over all destination cells. Two cells have similar diffusion states if these profiles are similar; they have different states if their probability mass is distributed differently. Diffusion pseudotime (DPT) turns these differences into a root-relative score: cells whose diffusion profiles are more different from those of a chosen root receive larger scores.

The root supplies an orientation, not a measured time origin in physical units. A DPT score is a graph-derived ordering that can help describe progression from a cell designated as the starting state. Its numerical scale depends on the graph, transition operator, and diffusion times included.

Definition: Let $P$ be a row-stochastic transition matrix on a finite set of cells, so $P_{ij}$ is the probability of moving from cell $i$ to cell $j$ in one step. Let $\pi$ be a stationary distribution with $\pi_j>0$, and let $r$ be a chosen root cell. At diffusion time $t$, define the profile distance between cells $r$ and $i$ by $$ D_t^2(r,i)=\sum_j \frac{\{P^t(r,j)-P^t(i,j)\}^2}{\pi_j}. $$ For a finite horizon $T$ and fixed nonnegative time weights $\omega_1,\ldots,\omega_T$, not all zero, define the cumulative root-based diffusion pseudotime score $$ \operatorname{DPT}_{T,\omega}(r,i)= \left\{\sum_{t=1}^{T}\omega_tD_t^2(r,i)\right\}^{1/2}. $$ Cells are ordered from the root by increasing score. This finite-horizon definition makes the included diffusion times explicit; implementations can use different horizons, weights, or related diffusion-map distances.

The factor $1/\pi_j$ weights discrepancies at destination cells by their stationary probability. In a reversible diffusion operator, it also makes the profile distance agree with the usual diffusion-distance geometry. The square root gives the score the same units as a profile distance; the score itself is not in hours, days, or any other biological-time unit.

Assumptions and a Spectral Justification

This construction needs a well-defined transition process and a meaningful set of cells over which to compare profiles. For the standard spectral interpretation, assume a finite, connected, undirected weighted cell graph, with a row-stochastic transition matrix $P$ that is reversible with respect to its positive stationary distribution $\pi$. Reversibility means $\pi_iP_{ij}=\pi_jP_{ji}$. A lazy walk, which retains some probability at the current cell, is one way to avoid periodic oscillation; it is not required to calculate a finite-horizon score, but it helps make long-run diffusion behavior stable. The graph and transition rule are treated as fixed for the calculation.

For a reversible $P$, take eigenfunctions $\psi_\ell$ orthonormal with respect to $\pi$, with eigenvalues $\lambda_\ell$. The stationary eigenfunction is constant and has eigenvalue 1. Reversibility gives the expansion $P^t(x,j)/\pi_j=1+\sum_{\ell\geq 2}\lambda_\ell^t\psi_\ell(x)\psi_\ell(j)$. Subtracting the profiles from $r$ and $i$ removes the constant term. Squaring and summing with weights $\pi_j$ eliminates cross-products between distinct orthonormal eigenfunctions. Thus

D_t^2(r,i) =\sum_{\ell\geq 2}\lambda_\ell^{2t} \{\psi_\ell(r)-\psi_\ell(i)\}^2, \qquad \operatorname{DPT}_{T,\omega}^2(r,i) =\sum_{\ell\geq 2} \left(\sum_{t=1}^{T}\omega_t\lambda_\ell^{2t}\right) \{\psi_\ell(r)-\psi_\ell(i)\}^2.

This derivation shows how time weights change the contribution of each spectral direction: directions that decay slowly can continue to distinguish cells over several steps. It also clarifies why DPT is a distance between diffusion profiles, not simply the number of graph edges from the root. The score is zero for the root, since every profile difference from the root to itself is zero. Other cells can also tie if their included profiles are identical.

Worked Example: Calculating Scores on a Three-Cell Graph

Consider three hypothetical cells $A,B,C$ arranged in a connected path. Use the transition matrix below, whose rows give the probabilities of moving from the cell named in that row. The chain is reversible with stationary distribution $\pi=(1/3,1/2,1/3)$: for example, $\pi_A P_{AB}=(1/3)(1/2)=1/6=\pi_B P_{BA}$, and the same balance holds on the $B$-to-$C$ edge.

P= \begin{pmatrix} 1/2&1/2&0\\ 1/3&1/3&1/3\\ 0&1/2&1/2 \end{pmatrix}, \qquad P^2= \begin{pmatrix} 5/12&5/12&1/6\\ 5/18&4/9&5/18\\ 1/6&5/12&5/12 \end{pmatrix}.

Take $A$ as the root, $T=2$, and equal weights $\omega_1=\omega_2=1$. At one step, the profile difference between $A$ and $B$ is $(1/6,1/6,-1/3)$. Therefore $$ D_1^2(A,B)=\frac{(1/6)^2}{1/3}+\frac{(1/6)^2}{1/2}+\frac{(-1/3)^2}{1/3} =\frac{17}{36}. $$ For $A$ and $C$, the one-step difference is $(1/2,0,-1/2)$, giving $$ D_1^2(A,C)=\frac{(1/2)^2}{1/3}+\frac{0^2}{1/2}+\frac{(-1/2)^2}{1/3} =\frac{3}{2}. $$

At two steps, the $A$-versus-$B$ difference is $(5/36,-1/36,-1/9)$. Hence $$ D_2^2(A,B)=\frac{(5/36)^2}{1/3}+\frac{(-1/36)^2}{1/2}+\frac{(-1/9)^2}{1/3} =\frac{125}{1296}. $$ The two-step difference between $A$ and $C$ is $(1/4,0,-1/4)$, so $$ D_2^2(A,C)=\frac{(1/4)^2}{1/3}+\frac{0^2}{1/2}+\frac{(-1/4)^2}{1/3} =\frac{3}{8}. $$

The root score for $B$ is $$ \operatorname{DPT}_{2,(1,1)}(A,B) =\sqrt{\frac{17}{36}+\frac{125}{1296}} =\sqrt{\frac{737}{1296}} \approx 0.7541. $$ For $C$, it is $$ \operatorname{DPT}_{2,(1,1)}(A,C) =\sqrt{\frac{3}{2}+\frac{3}{8}} =\sqrt{\frac{15}{8}} \approx 1.369. $$ Thus the ordering from root $A$ is $A$, then $B$, then $C$. This is a relative ordering under the chosen transition rule and horizon, not evidence that $C$ is a specified amount of biological time beyond $B$.

Root Choice and Time Horizon

Changing the root changes which profiles are compared and can change the ordering. This is why the root should represent the intended starting state, not merely a convenient cell. Even with a fixed root, including more diffusion times or changing their weights can alter the score magnitudes and potentially their ranking. The next example isolates the effect of the root on the same graph.

Worked Example: Moving the Root to the Middle Cell

Use the same hypothetical graph, but choose $B$ as the root and use only one diffusion step. The one-step profiles are the rows of $P$. The difference between $B$ and $A$ is $(-1/6,-1/6,1/3)$, so $$ D_1^2(B,A)=\frac{(-1/6)^2}{1/3}+\frac{(-1/6)^2}{1/2}+\frac{(1/3)^2}{1/3} =\frac{17}{36}, \qquad D_1(B,A)\approx0.6872. $$ The difference between $B$ and $C$ is $(1/3,-1/6,-1/6)$, which gives $$ D_1^2(B,C)=\frac{(1/3)^2}{1/3}+\frac{(-1/6)^2}{1/2}+\frac{(-1/6)^2}{1/3} =\frac{17}{36}, \qquad D_1(B,C)\approx0.6872. $$

The resulting one-step ordering is a tie: $B$ is the root and $A$ and $C$ have equal scores. The tie follows from the graph's symmetry around $B$, not from an algorithmic failure. If the biological question requires ordering $A$ and $C$, this root and these data do not provide that distinction.

A score often is reported without its units because it is a graph-derived quantity. If a display rescales scores—for example, dividing each by the largest score—that changes the numerical scale but not the ordering, provided the divisor is positive and common to every cell. Such rescaling does not turn pseudotime into elapsed time.

Choosing and Interpreting Time Weights

A finite horizon is explicit and avoids requiring an infinite sum to converge. Equal weights treat each included step's squared profile distance equally. Other nonnegative weights can emphasize early or later diffusion behavior. The choice should be stated because scores from different horizons or weights are not directly interchangeable. If a weight is zero, that time contributes nothing; if all weights are zero, the score is undefined as an ordering device.

Worked Example: Giving Later Diffusion Less Weight

Return to root $A$ in the three-cell example. Let $T=2$, $\omega_1=1$, and $\omega_2=1/4$, so that the second-step squared distance contributes one quarter as much as the first-step squared distance. For $B$, $$ \operatorname{DPT}_{2,(1,1/4)}^2(A,B) =\frac{17}{36}+\frac{1}{4}\frac{125}{1296} =\frac{2573}{5184}, \qquad \operatorname{DPT}_{2,(1,1/4)}(A,B)\approx0.7045. $$ For $C$, $$ \operatorname{DPT}_{2,(1,1/4)}^2(A,C) =\frac{3}{2}+\frac{1}{4}\frac{3}{8} =\frac{51}{32}, \qquad \operatorname{DPT}_{2,(1,1/4)}(A,C)\approx1.262. $$

The ordering remains $A,B,C$, but both nonroot scores are smaller than under equal weights. That outcome is specific to this graph and these profiles; a different weighting scheme need not preserve every ranking in a more complex graph.

Diagnostics, Alternatives, and Limitations

Before interpreting a root-based ordering, check the objects that define it. Each row of $P$ should sum to one, the stated stationary distribution should satisfy $\pi P=\pi$, and the weights in the profile distance must use that same $\pi$. For a reversible operator, check detailed balance as well. A mismatch between $P$ and $\pi$ changes the profile metric and can invalidate a claimed calculation. It is also useful to inspect whether the root has a plausible position and whether the ordering is stable to reasonable, explicitly stated choices of horizon and weights.

Common failure modes include a graph with weak links that make transition profiles misleading, a root chosen from the wrong cell state, and many cells receiving nearly indistinguishable scores. A disconnected graph prevents diffusion from reaching cells outside the root's component; profile comparisons across components should not be treated as an ordinary continuous trajectory ordering. Periodicity can make profiles oscillate at successive steps, so longer-horizon behavior deserves particular care when the operator is not made lazy. These checks flag concerns; a numerical score by itself cannot diagnose or repair them.

DPT differs from a shortest-path or minimum-spanning-tree pseudotime. A tree-path score accumulates edge lengths along a selected route, while DPT compares transition profiles and can reflect many possible walks through the graph. It also differs from using a single diffusion distance: the finite-horizon score combines distances across the specified time steps. Neither approach is universally preferable. Tree paths can be easier to interpret as route lengths but depend on the chosen tree; diffusion profiles capture broader graph connectivity but depend on the transition kernel and time weighting. Like other graph-based orderings, DPT is sensitive to how the graph is constructed and does not establish causation, cell ancestry, or actual elapsed time.

A frequent mistake is to read a DPT score of 2 as “twice as much developmental time” as a score of 1. The score is a weighted distance in diffusion-profile space, and its scale changes with the graph, stationary probabilities, weights, and horizon. The defensible interpretation is that, under the specified construction, one cell's diffusion profiles are more different from the root's than another's. A larger score supports a more root-distant position in this representation, not a direct claim about biological duration.

Key takeaway: DPT assigns cells root-relative scores by comparing their transition profiles across specified diffusion times. State the root, transition operator, stationary weights, horizon, and time weights; interpret the resulting order as a graph-dependent pseudotime, not a clock.

Check Your Understanding

Answer each question, then compare your reasoning with the brief solutions.

What does $P^t(i,\cdot)$ represent, and why compare it with the root's profile?
In the three-cell example with root $A$, what is $D_1^2(A,C)$, and which stationary probabilities weight the nonzero terms?
Why can two nonroot cells have the same DPT score?
What must be reported to make a root-based DPT calculation interpretable and reproducible?

Solutions. (1) It is the distribution over destination cells after $t$ transitions starting at $i$; comparison measures how differently the graph disperses probability from $i$ and the root. (2) It is $3/2$; the differences occur at $A$ and $C$, each weighted by $1/\pi=3$. (3) Their included profiles can be equally distant from the root, for example by graph symmetry, or can coincide at the selected times. (4) Report the root, transition operator and graph construction, stationary distribution used in the metric, and the diffusion horizon and time weights.