Mathematical Framework of RNA Velocity

Biological Background

Gene expression involves a series of biochemical processes:

1. Transcription: DNA → pre-mRNA (nascent/unspliced)
2. Splicing: pre-mRNA → mature mRNA (spliced)
3. Degradation: mRNA → degraded products

RNA velocity leverages the distinction between unspliced (nascent) and spliced (mature) mRNA to infer the direction and rate of gene expression changes.

The Kinetic Model

The dynamics of mRNA abundance can be described by a system of ordinary differential equations (ODEs):

\[\frac{du}{dt} = \alpha(t) - \beta u\]

\[\frac{ds}{dt} = \beta u - \gamma s\]

Where:

• \(u\): unspliced mRNA abundance
• \(s\): spliced mRNA abundance
• \(\alpha(t)\): transcription rate (time-dependent)
• \(\beta\): splicing rate (constitutive)
• \(\gamma\): degradation rate (constitutive)

Constant Transcription (Induction Phase)

When \(\alpha(t) = \alpha\) (constant), the system has analytical solutions:

Unspliced mRNA: \[u(t) = u_0 e^{-\beta t} + \frac{\alpha}{\beta}(1 - e^{-\beta t})\]

Spliced mRNA: \[s(t) = s_0 e^{-\gamma t} + \frac{\alpha}{\gamma}(1 - e^{-\gamma t}) + \frac{\alpha - u_0 \beta}{\gamma - \beta}(e^{-\gamma t} - e^{-\beta t})\]

Zero Transcription (Repression Phase)

When transcription stops (\(\alpha = 0\)):

\[u(t) = u_0 e^{-\beta t}\]

\[s(t) = s_0 e^{-\gamma t} + \frac{\beta u_0}{\gamma - \beta}(e^{-\gamma t} - e^{-\beta t})\]

Steady State

At equilibrium (\(du/dt = 0\), \(ds/dt = 0\)):

\[u_{ss} = \frac{\alpha}{\beta}\]

\[s_{ss} = \frac{\alpha}{\gamma}\]

The ratio at steady state defines a fundamental relationship:

\[\frac{u_{ss}}{s_{ss}} = \frac{\gamma}{\beta}\]

Phase portrait showing the trajectory of gene expression during induction (blue) and repression (red) phases.

Model 1: Deterministic (Steady-State)

Assumption: Cells are near transcriptional steady state.

At steady state: \(\gamma \approx \beta u / s\)

Velocity: The deviation from steady state indicates the direction of change:

\[v_s = \beta u - \gamma s\]

If \(v_s > 0\): gene is being upregulated If \(v_s < 0\): gene is being downregulated

Implementation in scVeloR:

# Fit gamma by linear regression of u vs s
# v_s = β·u - γ·s = β(u - γ/β·s)
# Use regression: u ~ s to get slope γ/β

Model 2: Stochastic

Assumption: Account for transcriptional noise through second-order moments.

The stochastic model uses both first moments (mean expression) and second moments (variance, covariance):

\[\text{Var}(u) = E[u^2] - E[u]^2\] \[\text{Cov}(u,s) = E[us] - E[u]E[s]\]

This provides more robust velocity estimates when transcriptional bursting is present.

Model 3: Dynamical

Assumption: Full kinetics model without steady-state assumption.

The dynamical model infers all kinetic parameters (\(\alpha\), \(\beta\), \(\gamma\)) and the switching time (\(t_-\)) using an Expectation-Maximization (EM) algorithm.

E-step: Assign latent time to each cell M-step: Update kinetic parameters

EM algorithm iteratively refines parameter estimates and time assignments.

Cosine Similarity

The velocity graph captures cell-cell transitions based on velocity correlation:

\[\text{cosine}(v_i, \delta_{ij}) = \frac{v_i \cdot \delta_{ij}}{||v_i|| \cdot ||\delta_{ij}||}\]

Where: - \(v_i\): velocity vector for cell \(i\) - \(\delta_{ij}\): expression difference between cell \(i\) and \(j\)

Transition Probability

The transition probability from cell \(i\) to cell \(j\):

\[P_{ij} = \frac{\max(0, \text{cosine}(v_i, \delta_{ij}))}{\sum_k \max(0, \text{cosine}(v_i, \delta_{ik}))}\]

Schematic of velocity-based transition probabilities.

Gene-Shared Latent Time

The latent time represents the position of each cell along the differentiation trajectory:

\[t_i = \frac{1}{|G|} \sum_{g \in G} t_{ig}\]

Where \(t_{ig}\) is the inferred time for cell \(i\) based on gene \(g\).

Root Cell Identification

The root cells are identified as those with: 1. High connectivity in the velocity graph 2. Velocity vectors pointing “outward” (high end-point probability)