RNA Velocity: Modeling Gene Expression Dynamics

• RNA velocity is a quantitative model that estimates transcriptional dynamics and cell fate transitions by analyzing unspliced and spliced RNA counts.
• It integrates kinetic modeling, Bayesian inference, and network-based approaches to overcome challenges like latent time resolution and computational scalability.
• Applications include mapping differentiation trajectories in development, disease modeling, and designing targeted interventions using control-theoretic principles.

RNA velocity is a quantitative model for gene expression dynamics that leverages single-cell RNA sequencing (scRNA-seq) data to infer the direction and rate of cell-state transitions. By distinguishing between unspliced (nascent) and spliced (mature) transcript counts at the single-cell level, it reconstructs underlying kinetic parameters and vector flows in gene expression space, serving as a foundational framework for mapping cell fate, lineage transitions, and cycling dynamics. Recent methodological advances address limitations of classical and widely-adopted approaches, incorporating rigorous Bayesian inference, network-theoretic modeling, and mathematically principled landscape construction.

1. Differential Kinetic Models and RNA Velocity Estimation

RNA velocity quantifies transcriptional dynamics by modeling the rates of transcription, splicing, and degradation for each gene. For cell $i$ and gene $g$, spliced ($s_i^g$) and unspliced ($u_i^g$) counts are observed. The canonical kinetic ODEs are: $$\frac{du}{dt} = \alpha - \beta u, \qquad \frac{ds}{dt} = \beta u - \gamma s,$$ with $\alpha$ as the transcription rate, $\beta$ the splicing rate, $\gamma$ the degradation rate, and the instantaneous RNA velocity given by $v_i^g \approx \beta u_i^g - \gamma s_i^g$ (Jia et al., 2023).

Bayesian frameworks, such as BayVel, further refine this estimation by modeling raw count data via a hierarchical Negative-Binomial structure, incorporating latent group and subgroup variables and fixing global scaling invariances (e.g., $\beta_g=1$, average capture efficiency $\lambda_c=1$). Velocity is derived for each posterior sample as $$v_{k,r,g} = \beta_g u_{k,r,g}^\sim - \gamma_g s_{k,r,g}^\sim$$ (Sabbioni et al., 6 May 2025).

2. Vector Field Construction and Natural Helmholtz–Hodge Decomposition

Discrete cell-wise velocities ${v_i}$ are mapped onto a continuous vector field $$F(x): \Omega \subset \mathbb{R}^n \rightarrow \mathbb{R}^n$$ (after projection via dimension reduction such as UMAP). The Natural Helmholtz–Hodge Decomposition (nHHD) theorem guarantees a unique partition for sufficiently smooth $F(x)$ within a bounded domain: $F(x) = F_g(x) + F_c(x) + F_h(x),$ where $F_g = -\nabla U(x)$ (gradient, curl-free), $F_c = \nabla \times A(x)$ (rotational, divergence-free), and $F_h$ harmonic ($\Delta F_h = 0$). This decomposition enables explicit distinction of cell differentiation (gradient descent in $U(x)$) and cell cycling (rotation in $A(x)$) in the Waddington landscape metaphor (Jia et al., 2023). Computationally, $U(x)$ and $A(x)$ are identified by solving Poisson-type PDEs or by line integration and convolution with Green's functions on appropriate grids.

3. Network-Based Extensions: Intracellular GRNs and Intercellular Coupling

Recent developments generalize RNA velocity models to account for regulatory context. In a multi-gene, multi-cell system, transcription rates are modulated by an intracellular gene regulatory network (GRN) specified by activating ($W^+$) and repressing ($W^-$) matrices. The regulatory input for gene $g$ is

$$R_g(s) = \frac{\kappa + \sum_{q=1}^{n_g} W_{gq}^{+} s^q}{\kappa + \sum_{q=1}^{n_g} W_{gq}^{-} s^q},$$

where $\kappa$ is a regularization constant. The full spatially-coupled dynamic system incorporates cell–cell adjacency $A$ and a coupling coefficient $c$: $$\frac{du_i^g}{dt} = \alpha_i^g R_g(s_i) - \beta_i^g u_i^g, \quad \frac{ds_i^g}{dt} = \beta_i^g u_i^g - \gamma_i^g s_i^g + c \sum_{j=1}^{n_c} A_{ij}(s_j^g - s_i^g)$$ Equilibrium existence and global stability can be established via spectral radius criteria and Lyapunov functions. The algebraic connectivity $\lambda_2(L)$ of the Laplacian $L = D - A$ determines consensus strength across the cell population (Hou et al., 3 Jan 2026).

4. Algorithmic and Statistical Identification: Time-Scale Fixation and Uncertainty

RNA velocity is subject to nonidentifiability under time-rescaling ($t \to \kappa t$, $\alpha \to \alpha/\kappa$, etc.), demanding post-hoc resolution of latent time parameters. Rational fixation recovers global latent times $t_c$ and true splicing rates $\beta_g$ via constrained maximum likelihood, using either multiplicative ($t_{cg} \beta_g^{-1} = t_c + \epsilon_{cg}$) or additive ($t_{cg} = t_c \beta_g + \epsilon_{cg}$) noise models. Solutions involve Perron eigenvector computation in gene-time matrices (Li et al., 2023).

Bayesian methods further provide principled uncertainty quantification for kinetic rates, latent times, and velocities, propagating posterior samples to credible intervals and posterior predictive diagnostics (Sabbioni et al., 6 May 2025). EM-based frequentist approaches use observed and missing information matrices, expressed via the SEM (Supplemented EM) decomposition, to obtain asymptotic covariance estimates.

5. Downstream Analysis: Transition Kernels, Bandwidth Selection, and Pseudotemporal Distances

To infer developmental trajectories, velocity-induced random walks are constructed via Gaussian-cosine kernels in spliced-count space. The Markov chain transition kernel $P_{\varepsilon,n}$ is tuned for operator convergence; optimal kernel bandwidth is determined by balancing variance ($O(n^{-1/2} \varepsilon^{-d/4})$) and bias ($O(\sqrt{\varepsilon})$), yielding $\varepsilon_{\rm opt} = n^{-2/(d+2)}$ (Li et al., 2023).

Pseudotemporal ordering employs mean first-hitting times in these Markov chains. For cell $i$ and target cluster $A$,

$$k_i^A = \mathbb{E}[\min \{ n \ge 0: X_n \in A \} | X_0 = i]$$

is computed recursively or by linear solvers. For bifurcating or cyclic systems, taboo sets $H$ enforce trajectory restriction, adjusting hitting time computation for lineage-specific transitions.

6. Applications, Comparative Evaluation, and Control-Theoretic Implications

RNA velocity approaches have been validated on synthetic and real scRNA-seq datasets, including dentate gyrus and hematopoiesis (gradient-dominated flows, minimal cycling), and pancreatic development (simultaneous gradient and rotational flows in cycling progenitors) (Jia et al., 2023, Sabbioni et al., 6 May 2025). Bayesian methods, notably BayVel, outperform deterministic ODE-fitting methods (scVelo) in statistical accuracy and uncertainty quantification, especially under simulation.

Control-theoretic formulations allow for targeted intervention in single-cell GRNs and multicellular systems. Drug-regimen optimization, expressible as minimum-time Hamiltonian problems, yields explicit bang–bang control policies, leveraging system reachability as determined via Lie bracket analysis. These frameworks suggest robust estimation and rational perturbation design for spatial transcriptomics and disease-modeling contexts (Hou et al., 3 Jan 2026).

7. Limitations, Practical Workflow, and Outlook

Although recent RNA velocity methodologies provide rigorous frameworks for kinetic inference and landscape construction, they pose significant computational burden (e.g., MCMC convergence over thousands of genes/cells), identifiability challenges (latent time determination, subgrouping), and dependence on model assumptions (e.g., piecewise-constant ON/OFF transcription, kernel choice). Practical workflows proceed from raw count acquisition, grouping/subgrouping, prior specification, batch inference, and downstream visualization in embedding space. Model selection is performed via WAIC and careful evaluation of clustering granularity.

These advances offer enhanced quantification of developmental potency and cycling strength, accurate mapping of bifurcations and cell-fate transitions, and a testbed for perturbation-optimized control in biological systems. The unified dual-landscape and network-consensus frameworks situate RNA velocity as an essential analytic tool for dynamic, high-dimensional transcriptomic systems.