🧬 Dynamic RNA velocity model-- (4) latent time

Introduction

In this fourth installment of our blog series on effectively applying the dynamic model to infer RNA velocity from single-cell RNA-seq, we reveal the mathematical foundations of latent time in the context of RNA velocity analysis, which enables in-depth grasp and interpretation of the latent time of RNA velocity.

Latent time is a crucial interpretation of RNA velocity independent of pre-defined (never perfect) dimentional reduction and it embedding. And latent time provide expressive visulization of the temporal dynamics of cell differentiation according to RNA velocity.

Here specifically focuses on:

Identification of root cells in section A)
Discussion and numeric examples of transport maps in section B)
Gene and cell dependence of latent time in section C)
Neighborhood-convolved latent time in section D)

A) Identify Starting Points Using Markov Chain Analysis

1. What Are Root Cells in Differentiation

In single-cell RNA velocity analysis, root cells are the initial states in a developmental trajectory. These are cells that:
✅ Have low latent time (early progenitors).
✅ Are unlikely to have transitioned from other states.
✅ Serve as the starting points for differentiation pathways.

Identifying these root cells is crucial for understanding how differentiation originates.

2. Stationary State Equation in RNA Velocity

The root cells are found by solving the stationary distribution equation for the transition probability matrix $ \tilde{\pi} $ :

$$ \mu^* = \mu^* \tilde{\pi}^T $$

where:

$ \tilde{\pi}^T $ is the transpose of the transition probability matrix that encodes RNA velocity-driven transitions.
$ \mu^* $ represents the stationary distribution, which describes the long-term probabilities of cells staying in each state.
The left eigenvectors of $ \tilde{\pi}^T $ corresponding to the eigenvalue of 1 define these stationary states.

3. Why Does This Identify Root Cells

a. Probability Evolution in a Markov Process

In a discrete-time Markov chain, the state of the system at time $ t $ is represented by a probability distribution $ \mu_t $ , which evolves according to the transition probability matrix $ P $ :

$$ \mu_{t+1} = \mu_t P $$

where:

$ \mu_t $ is a row vector of probabilities over all possible states at time $ t $ .
$ P $ is an $ n \times n $ transition matrix, where $ P_{ij} $ represents the probability of transitioning from state $ i $ to state $ j $ .

Each step updates the probability distribution, progressively altering the likelihood of being in each state.

b. Convergence to a Stable Distribution

As the Markov chain progresses, repeated multiplications by $ P $ lead the probability distribution toward a stable equilibrium:

$$ \mu_{t+k} = \mu_t P^k, \quad \text{as } k \to \infty $$

After many steps, the probabilities stabilize, meaning the system reaches a final long-term behavior where further transitions do not significantly change the probabilities.

c. Definition of the Stationary Distribution

The stationary distribution $ \mu^* $ is a probability vector that remains unchanged under transitions:

$$ \mu^* = \mu^* P $$

This implies that if a system starts in $ \mu^* $ , applying the transition matrix does not alter its probabilities–it stays in equilibrium.

d. Interpretation of Stationary Distribution

✅ Predicts long-term state occupancy–the fraction of time spent in each state.
✅ Defines steady-state probabilities in applications like RNA velocity analysis.
✅ Identifies stable states (e.g., progenitor/root cells in differentiation) when applied to biological modeling.

4. Compute Stationary States Using Eigenvectors

Step 1: Define the Transition Probability Matrix

Consider a 3-state system with transition probabilities:

$$ \tilde{\pi} = \begin{bmatrix} 0.8 & 0.1 & 0.1 \\ 0.2 & 0.7 & 0.1 \\ 0.1 & 0.2 & 0.7 \end{bmatrix} $$

This matrix describes how probabilities shift between states.

Step 2: Find the Left Eigenvector for Eigenvalue 1

We solve:

$$ v^T \tilde{\pi}^T = v^T $$

This means we need the eigenvector corresponding to eigenvalue $ \lambda = 1 $.

Computing the eigenvectors of $ \tilde{\pi}^T $ , we obtain:

$$ v = \begin{bmatrix} 0.45 \\ 0.35 \\ 0.20 \end{bmatrix} $$

This left eigenvector represents the stationary distribution, meaning: ✅ State 1 holds 45% of the long-term probability.
✅ State 2 holds 35%.
✅ State 3 holds 20%.

These steady-state values describe how the system behaves in the long run, identifying the states least likely to transition further.

Summery & Interpretation

In single-cell analysis:

Cells corresponding to the stationary distribution act as root cells, initiating differentiation.
This eigenvector-driven method ensures an objective, data-driven way to define progenitor states.
Eigenvalue 1 captures states that remain stable over differentiation, meaning they are starting points in biological development.

B) Details and Numeric Examples of Transport Maps for Cell Development Analysis

1. Mathematical Notation

$\mu = (\mu_1, \ldots, \mu_n)$ : a distribution over cells

$\pi_e$ : a transport map, e.g., inferred from RNA velocity

$\tilde{\pi}^T$ : the transpose (or reverse) of the transport map

$\mu_{\text{des}} = \mu \cdot \pi_e$ : the descendant distribution

$\mu_{\text{anc}} = \mu \cdot \tilde{\pi}^T$ : the ancestor distribution

2. Core Concepts

Forward Transport: “Where are cells going?”

A distribution of cells $\mu$ can be pushed through the transport map to determine future cell states.

Starting with a distribution over some cells $\mu$ , for example:

All cells in a particular cluster or cell type
Or a uniform distribution over some selected cells

Pushing forward this distribution through the transport map $\pi_e$ gives the distribution of where these cells go in the future:

$$\mu_{\text{des}} = \mu \cdot \pi_e \tag{1}$$

This answers the question: “Given that my cells are here now, where are they going?”

Reverse Transport: “Where did cells come from?”

Conversely, a distribution $\mu$ can be pulled back to determine cell origins.

You take a set of cells at a later time or later state, and ask: “Where did they come from?”

You pull back the distribution via the transpose of the transport matrix:

$$\mu_{\text{anc}} = \mu \cdot \tilde{\pi}^T \tag{2}$$

The transpose reverses the direction of flow – it computes how much each earlier state contributed to the current cell distribution.

Indicator-Based Cell Selection

For a set of cells $S = \{s_1, \ldots, s_n\}$ , the distribution $\mu$ is defined using an indicator vector:

$$\mu_i = \mathbf{1}_{s_i \in S} \tag{3}$$

This means:

$\mu_i = 1$ if cell $i$ is in the set $S$
$\mu_i = 0$ otherwise

The input distribution puts mass 1 on each selected cell, and 0 elsewhere.

3. Biological Cases

Progenitor Analysis:

Select a set of progenitor cells (set $S$ )
Compute $\mu_{\text{des}} = \mu \cdot \pi_e$ → get the expected distribution of their descendants

Lineage Tracing:

Select differentiated cells
Compute $\mu_{\text{anc}} = \mu \cdot \tilde{\pi}^T$ → get their likely origins

This framework helps answer:

“Where are these cells going?”
“Where did these cells come from?”

4. Numerical Example: Forward Transport

Consider 4 cells ( $n=4$ ), selecting cells 1 and 3:

$$\mu = [1, 0, 1, 0] \tag{4}$$

Assume a learned transport map:

$$\pi_e = \begin{bmatrix} 0.6 & 0.2 & 0.1 & 0.1 \\ 0.1 & 0.7 & 0.1 & 0.1 \\ 0.2 & 0.1 & 0.6 & 0.1 \\ 0.3 & 0.2 & 0.1 & 0.4 \end{bmatrix} \tag{5}$$

Then the descendant distribution is:

$$\mu_{\text{des}} = \mu \cdot \pi_e = [1, 0, 1, 0] \cdot \pi_e = \text{row 1} + \text{row 3} \tag{6}$$

Computing the result:

$$\mu_{\text{des}} = [0.6 + 0.2, 0.2 + 0.1, 0.1 + 0.6, 0.1 + 0.1] = [0.8, 0.3, 0.7, 0.2] $$

This gives a distribution over where the selected cells (1 & 3) are headed – their descendant cell states.

5. Numerical Example: Reverse Transport

Step 1: Compute the Ancestor Distribution

To find the ancestor distribution, we use the transpose of $\pi_e$ :

$$(\pi_e)^T = \begin{bmatrix} 0.6 & 0.1 & 0.2 & 0.3 \\ 0.2 & 0.7 & 0.1 & 0.2 \\ 0.1 & 0.1 & 0.6 & 0.1 \\ 0.1 & 0.1 & 0.1 & 0.4 \end{bmatrix} \tag{3}$$

Step 2: Calculate the Ancestor Distribution

The ancestor distribution is computed as:

$$\mu_{\text{anc}} = \mu_{\text{des}} \cdot (\pi_e)^T \tag{4}$$

Multiplying the row vector $\mu_{\text{des}} = [0.8, 0.3, 0.7, 0.2]$ by the matrix $(\pi_e)^T$ :

$$\mu_{\text{anc},1} = 0.8 \times 0.6 + 0.3 \times 0.1 + 0.7 \times 0.2 + 0.2 \times 0.3 = 0.48 + 0.03 + 0.14 + 0.06 = 0.71 $$ $$\mu_{\text{anc},2} = 0.8 \times 0.2 + 0.3 \times 0.7 + 0.7 \times 0.1 + 0.2 \times 0.2 = 0.16 + 0.21 + 0.07 + 0.04 = 0.48 $$ $$\mu_{\text{anc},3} = 0.8 \times 0.1 + 0.3 \times 0.1 + 0.7 \times 0.6 + 0.2 \times 0.1 = 0.08 + 0.03 + 0.42 + 0.02 = 0.55 $$ $$\mu_{\text{anc},4} = 0.8 \times 0.1 + 0.3 \times 0.1 + 0.7 \times 0.1 + 0.2 \times 0.4 = 0.08 + 0.03 + 0.07 + 0.08 = 0.26 $$

Step 3: Resulting Ancestor Distribution

$$\mu_{\text{anc}} = [0.71, 0.48, 0.55, 0.26] \tag{9}$$

Interpretation: Given that the descendant distribution was $[0.8, 0.3, 0.7, 0.2]$ , pulling it back through the transpose of the transport map estimates the likely ancestor distribution of these descendants. The weights show which earlier states contributed more to this descendant cell population.

Normalization

We can normalize $\mu_{\text{anc}}$ to turn it into a proper probability distribution:

$$\text{Sum} = 0.71 + 0.48 + 0.55 + 0.26 = 2.00 $$ $$\mu_{\text{anc,normalized}} = [0.355, 0.24, 0.275, 0.13] $$

6. Why Perfect Inversion is Impossible

The descendant distribution $[0.8, 0.3, 0.7, 0.2]$ cannot lead exactly back to the original indicator vector $[1, 0, 1, 0]$ when using the transpose of the transport map for several fundamental reasons:

Non-invertibility of the Transport Map

The transport matrix $\pi_e$ typically represents a probabilistic, many-to-many transition between cell states. It’s a stochastic matrix where rows sum to 1, but it’s often not invertible or the inverse is not unique. The operation of pushing forward loses information through mixing and averaging, and pulling back is not a perfect inverse.

Smoothing and Mixing

When we multiplied $[1, 0, 1, 0]$ by $\pi_e$ , we obtained a weighted average of rows 1 and 3, resulting in the smoothed vector $[0.8, 0.3, 0.7, 0.2]$ . This vector is continuous-valued and represents a distribution over many possible descendant states. The reverse step with the transpose spreads these probabilities back to ancestors but cannot undo the smoothing perfectly.

Loss of Sharpness

The original vector $[1, 0, 1, 0]$ is a sharp indicator (only cells 1 and 3 are selected). After forward mapping, the distribution reflects uncertainty or spreading. Pulling back cannot perfectly re-sharpen because the transport matrix mixes information.

Biological Analogy

In biological systems, a population of cells can differentiate into multiple descendant states probabilistically. The future cell states reflect a blend of many ancestors. Hence, you cannot simply invert descendants back to a unique set of ancestors – instead, you get a probabilistic ancestor distribution that represents the uncertainty inherent in the biological process.

C) Gene and cell dependence in the calculation of latent time

scVelo’s latent time is calculated as:

$$t_{i;o} = Q^p_g \left( t_{i;g} - t_{o;g} \right)$$

Where:

$t_{i;o}$ : latent time for cell $i$ relative to root cell $o$
$Q^p_g$ : $p$ -quantile computed for gene $g$
$t_{i;g} - t_{o;g}$ : time shift between cell $i$ and root cell $o$ for gene $g$

Please note $t_{i;g}$ and $t_{o;g}$ have been normalized to a common overall time scale, which worths further discussion in the next blog.

1. Is p Gene-Specific?

❌ No – $p$ is not gene-specific.

Explanation

$Q^p_g$ means: for each gene $g$ , you compute the $p$ -quantile over the time shifts $t_{i;g} - t_{o;g}$ , across all cells $i$ .

But the value of $p$ itself is the same for all genes – it is not gene-specific.

So even though the quantile is applied per gene (i.e., each gene gets its own $Q^p_g$ ), the $p$ used in all those calculations is shared.

Mathematical Illustration

For a fixed value $p = 0.5$ (median):

$$Q^{0.5}_{gene1} = \text{median}(\{t_{i;gene1} - t_{o;gene1}\}_{i=1}^N)$$ $$Q^{0.5}_{gene2} = \text{median}(\{t_{i;gene2} - t_{o;gene2}\}_{i=1}^N)$$ $$\vdots$$ $$Q^{0.5}_{geneM} = \text{median}(\{t_{i;geneM} - t_{o;geneM}\}_{i=1}^N)$$

The same $p = 0.5$ is used for all genes, but each gene gets its own quantile value.

2. Is p Cell-Specific?

❌ Also no – $p$ is not adjusted for each cell. It is only used to transform the gene-wise temporal shifts into a consistent time estimate for each cell, relative to the root cell $o$ .

Explanation

The only place where indexing over cells happens is in computing the latent time $t_{i;o}$ , but the value of $p$ stays fixed for all cells.

Mathematical Process

Fixed $p$ across all computations: $p = p_{\text{fixed}}$
Gene-wise quantile computation: $$Q^{p_{\text{fixed}}}_g = \text{quantile}_{p_{\text{fixed}}}(\{t_{i;g} - t_{o;g}\}_{i=1}^N)$$
Cell-specific time assignment: $$t_{1;o} = Q^{p_{\text{fixed}}}_g \left( t_{1;g} - t_{o;g} \right)$$ $$t_{2;o} = Q^{p_{\text{fixed}}}_g \left( t_{2;g} - t_{o;g} \right)$$ $$\vdots$$ $$t_{N;o} = Q^{p_{\text{fixed}}}_g \left( t_{N;g} - t_{o;g} \right)$$

3. Is p Root-Specific?

✅ Yes – scVelo searches for the best $p$ per root candidate to optimize time smoothness.

Optimization Process

For each potential root cell $o$ , scVelo optimizes:

$$p^*_o = \arg\min_p \mathcal{L}_{\text{smoothness}}(p, o)$$

Where the loss function might be:

$$\mathcal{L}_{\text{smoothness}}(p, o) = \sum_{i,j \in \text{neighbors}} \left| t_{i;o}(p) - t_{j;o}(p) \right|^2$$

Root Selection

The final root and $p$ are chosen as:

$$o^*, p^* = \arg\min_{o,p} \mathcal{L}_{\text{smoothness}}(p, o)$$

4. Summary Table

Aspect	Is $p$ dependent?	Explanation
Per gene $g$	❌ No	The quantile $Q^p$ is computed per gene, but $p$ itself is fixed
Per cell $i$	❌ No	All cells use the same $p$ ; only their computed values differ
Per root $o$	✅ Yes	scVelo searches for the best $p$ per root candidate to optimize time smoothness

5. Optimization Step

Once scVelo computes all $ t_{i,o} $ for different candidate root cells $ o $ , it looks for the root and quantile $ p $ that maximize correlation between the latent time series
$$ \{ t_{i,o} \} $$
and its convolution over a local neighborhood graph (e.g., KNN graph of cells).

Mathematical Formulation

The optimization step seeks:

$$ \arg\max_{o,p} \text{corr} \left( t_{i,o}, \sum_{j \in N(i)} w_{ij} t_{j,o} \right) $$

where:

$ N(i) $ is the neighborhood of cell $ i $ .
$ w_{ij} $ represents graph weights, determining connectivity between cells.

D) Neighborhood-convolved Latent Time in scVelo

Step-by-step explanation

1. Latent time

scVelo computes a latent time $t_n$ for each cell $n$ , indicating how far along a differentiation trajectory it is, starting from the inferred root cells.
This time is estimated based on how consistent a cell’s transcriptional state is with the direction of gene expression changes (RNA velocity).

2. Gene-shared latent time course

This is the original latent time vector computed across all cells.
It’s “gene-shared” because it’s not just one gene’s trajectory – it aggregates information across genes to produce a common timeline.

3. Neighborhood convolution

scVelo smooths the latent time values using k-nearest neighbor (KNN) information (e.g., based on transcriptomic similarity).
The convolved time for a cell is the average of the latent times of its neighboring cells.

Think of it as:

$$\tilde{t}_n = \frac{1}{|N(n)|} \sum_{m \in N(n)} t_m \tag{1}$$

where $N(n)$ is the neighborhood of cell $n$ .

4. Regression & consistency check

scVelo performs linear regression: it tries to predict each cell’s latent time using its convolved value.
If a cell’s latent time significantly deviates from the trend of its neighbors, it is likely noisy or unreliable.

5. Correction / Replacement

For these outlier cells, scVelo replaces their latent time with the convolved (smoothed) version.
This acts like a denoising step – ensuring latent time values are locally smooth and robust to noise.

Numeric Explanation of Neighborhood Convolution

Let’s say we have 5 cells, and we already computed latent times for them as:

Cell	Latent Time $t_n$
A	0.10
B	0.15
C	0.85 ❗
D	0.20
E	0.25

Clearly, cell C looks suspicious – it’s way ahead in time compared to its neighbors.

Step 1: Define neighbors

Let’s say the neighborhoods (e.g., via KNN) are:

A → neighbors: [B, D]
B → neighbors: [A, D]
C → neighbors: [B, D, E]
D → neighbors: [A, B, E]
E → neighbors: [C, D]

Step 2: Compute convolved (smoothed) latent times

$$\tilde{t}_n = \text{average of neighbor latent times} \tag{2}$$

Cell	Neighbors	Neighbor Times	Convolved Time $\tilde{t}_n$
A	B, D	0.15, 0.20	0.175
B	A, D	0.10, 0.20	0.15
C	B, D, E	0.15, 0.20, 0.25	0.20
D	A, B, E	0.10, 0.15, 0.25	0.167
E	C, D	0.85, 0.20	0.525 ❗

Step 3: Compare original vs convolved times

Cell	$t_n$	$\tilde{t}_n$	Difference
A	0.10	0.175	0.075
B	0.15	0.15	0.000
C	0.85 ❗	0.20	0.65 ❗
D	0.20	0.167	0.033
E	0.25	0.525 ❗	0.275 ❗

Step 4: Detect inconsistent cells

scVelo runs a regression of $t_n \sim \tilde{t}_n$ and checks for large residuals.
Cells C and E have large mismatches between original and smoothed values → likely noisy or incorrect.

Step 5: Replace inconsistent values

Let’s say we define a cutoff: if difference > 0.3, replace with $\tilde{t}_n$ .

Cell	New $t_n$
A	0.10 (unchanged)
B	0.15 (unchanged)
C	0.20 (replaced!)
D	0.20 (unchanged)
E	0.525 (replaced!)

Final corrected latent time:

$$\mathbf{t} = [0.10, 0.15, \mathbf{0.20}, 0.20, \mathbf{0.525}] \tag{3}$$

Cell C is now brought back in line with its local trajectory.