Comprehensive RNA velocity by modeling the cascade of gene regulation, transcription, and splicing from single-cell RNA sequencing data with TSvelo

We present TSvelo, a method for jointly estimating RNA velocity across high-dimensional genes by integrating both transcriptional regulation and splicing information. Existing RNA velocity models predominantly focus on modeling the splicing process and typically operate in a gene-wise manner. In contrast, TSvelo explicitly models the full cascade of regulation, transcription, and splicing, while capturing their coordinated dynamics across all genes simultaneously. A systematic comparison between TSvelo and representative RNA velocity methods is provided in Appendix 1—table 1. In TSvelo, the unspliced and spliced RNA abundances are initially preprocessed for velocity gene selection and pseudotime initialization (Figure 1a). See Methods for details about preprocessing. Next, to model the velocity gene $g$’s expression dynamics, we suppose $u_g$ and $s_g$ are the abundance of unspliced and spliced RNA, and ${\alpha }_g\left(t\right),{\beta }_g$, and ${\gamma }_g$ are the transcription, splicing, and degradation rates, respectively. The dynamics is modeled as

Furthermore, we assume that the gene- and cell-specific transcriptional rate ${\alpha }_g\left(t\right)$ is influenced by the expression of TFs.

Figure 1 with 7 supplements see all

Download asset Open asset

Considering the wide usage and interpretability of linear models for gene relations in previous studies (Li et al., 2024a; Tran et al., 2022; Song et al., 2023; Liu et al., 2011; Wang et al., 2023), we model ${\alpha }_g\left(t\right)$ as

$ReLU$ denotes the Rectified Linear Unit activation function, defined as $ReLU\left(x\right)=max(0,x)$. The term $w_{gi}$ represents the regulatory weight of TF $i$ on the target gene $g$. $TFs\left(g\right)$ refers to the TFs that could regulate the target gene $g$, which are selected according to the ChEA and ENCODE TF–target databases. In order to include more TFs in the model, we also reserve those TFs that are not selected as velocity genes and directly model the dynamics between transcription and spliced RNA without using the unspliced abundance. Finally, we combine the dynamic models on all selected genes into the ODE matrix form (Figure 1b), which enables directly inferring a unified latent time for each cell using all its genes. See Methods for details.

TSvelo employs an EM framework to iteratively optimize both latent time $t$ and the parameters in ODE. The global pseudotime $t$ is assigned to each cell through grid search. Due to the difficulty in calculating analytical solutions of $u_g\left(t\right)$ or $s_g\left(t\right)$, TSvelo adopted a Neural ODE for estimating those parameters of transcription rate, splicing rate, and degradation rate. Next, TSvelo can model the gene dynamics of both transcription and splicing processes (Figure 1c), predict cell states using pseudotime and velocity stream analysis (Figure 1d), and is also applicable to multi-lineage tasks for analyzing expression patterns across different lineages in large complex scRNA-seq datasets (Figure 1e). We further validated TSvelo on simulated datasets and observed consistent improvements in velocity estimation and trajectory reconstruction across a range of settings (Figure 1—figure supplements 1 and 2).

To validate the TSvelo model, we applied it to the pancreas scRNA-seq dataset (Bastidas-Ponce et al., 2019), which is widely used in RNA velocity studies to model cell differentiation from ductal to endocrine cells (La Manno et al., 2018). Using the predicted RNA velocity, TSvelo can generate velocity stream plot using the functions provided in scVelo as well. Both the pseudotime $t$ (Figure 2a) and the velocity stream plot (Figure 2b) effectively capture the cell differentiation process. We quantitatively compare TSvelo and baseline approaches, including scVelo (Bergen et al., 2020), dynamo (Qiu et al., 2022), UniTVelo (Gao et al., 2022), cellDancer (Li et al., 2024b), and TFvelo (Li et al., 2024a).

Figure 2 with 3 supplements see all

Download asset Open asset

Figure 2—source data 1

Source data of the comparison between the 2D and 3D phase portrait fitting.

https://cdn.elifesciences.org/articles/108950/elife-108950-fig2-data1-v1.csv

Download elife-108950-fig2-data1-v1.csv

In the conventional 2D unspliced–spliced phase portrait, cells from different clusters often overlap, limiting separability. By introducing the latent variable α, the representation is extended to a 3D space, which helps disentangle these mixed states and reveals a clearer phase structure. As shown in Figure 2c, the 3D representation achieves consistently higher k-nearest neighbor (kNN) classification accuracy for cell-state separation than the 2D u–s embedding (one-sided Mann–Whitney U test, p = 4.37 × 10⁻¹⁰). Details are provided in Methods. These results indicate that the 3D phase portrait provides improved separation of cell states and a stronger foundation for modeling the underlying dynamics along cell trajectories.

Additional quantitative evaluation is provided in Figure 2—figure supplement 1. TSvelo achieves the highest median velocity consistency, which demonstrates that the high-dimensional velocity vectors learned by TSvelo are mostly coherent within neighbor cells. TSvelo also achieves the highest median in-cluster coherence and cross-boundary correctness, which validates that TSvelo best fits the differentiation process within these cell types according to the ground-truth annotation.

Next, we show the dynamics for individual genes and demonstrate how the TSvelo model aids gene dynamics analysis by incorporating both transcriptional and splicing information. Maml3, Anxa4, and Gstz1 are selected as examples because their unspliced–spliced 2D phase portrait exhibits mixed or overlapping patterns that are difficult to model using conventional RNA velocity approaches. Figure 2d presents the results of the TSvelo model on these genes, and Figure 2e presents the results of baseline approaches. Many previous approaches assume that the unspliced–spliced phase portrait exhibits an almond-shaped distribution (Weiler et al., 2022), which may not actually hold in the data. For Maml3, while the unspliced–spliced distribution predominantly follows the almond-shaped pattern, some cells, specifically Alpha cells (blue) and Beta cells (light blue), overlap with Ngn3 high EP cells (yellow), indicating that these cell types cannot be distinctly separated using only the unspliced–spliced 2D phase portrait. Nevertheless, TSvelo can infer transcription representation from the expression of multiple TFs, which enables the model to distinguish these cell types (Figure 2d, f). Anxa4 shows higher expression in Ductal cells (in red) compared to Ngn3 low EP cells (in orange), which means its expression pattern exhibits an initial decrease followed by an increase. Such dynamics are not easily captured in the conventional unspliced–spliced phase portrait used by previous approaches, as many baseline methods implicitly assume a decreasing–then–increasing expression pattern. By comparison, TSvelo can still fit such expression patterns by using additional information from the 3D phase portrait. The visualization of dynamics fitting on more genes is provided at Figure 2—figure supplement 2 and Figure 2—figure supplement 3.

We applied TSvelo to the gastrulation erythroid dataset, which is derived from the transcriptional profile of mouse embryos (Pijuan-Sala et al., 2019) and has been used in previous RNA velocity studies. This dataset primarily describes the differentiation process from blood progenitors to erythroid cells. First, using Gene Ontology (GO) term enrichment analysis (Figure 3a), we find that the velocity genes selected through TSvelo’s preprocessing strategies are mostly enriched in the erythropoiesis-related process, providing a preliminary validation of the biological plausibility of the inferred velocities.

Figure 3 with 2 supplements see all

Download asset Open asset

Figure 3—source data 1

The comparison of velocity consistency across different methods.

https://cdn.elifesciences.org/articles/108950/elife-108950-fig3-data1-v1.zip

Download elife-108950-fig3-data1-v1.zip

Figure 3—source data 2

The comparison of in-cluster coherence across different methods.

https://cdn.elifesciences.org/articles/108950/elife-108950-fig3-data2-v1.zip

Download elife-108950-fig3-data2-v1.zip

Figure 3—source data 3

The comparison of cross-boundary direction correctness across different methods.

https://cdn.elifesciences.org/articles/108950/elife-108950-fig3-data3-v1.zip

Download elife-108950-fig3-data3-v1.zip

After applying TSvelo, both the inferred pseudotime (Figure 3b) and the velocity stream plot (Figure 3c) are compatible with current biological knowledge about the cell differentiation process. We also compare TSvelo with previous approaches. The results, shown in Figure 3d–f, demonstrate that TSvelo achieves the highest velocity consistency, the highest in-cluster coherence, and the highest cross-boundary correctness, showing that the high-dimensional dynamics predicted by TSvelo best capture the underlying biological processes.

We next analyze TSvelo modeling at the gene level in Figure 3g–i. The 3D phase portrait again provides a better basis for dynamic fitting, as it more effectively separates cells from different clusters than the traditional 2D phase portrait (Figure 3—figure supplement 1). For genes exhibiting higher noise or more complex patterns, TSvelo demonstrates its robustness in dynamic fitting by integrating both transcriptional and splicing signals. For Hsp90ab1, which exhibits a counter-clockwise pattern in the unspliced–spliced phase portrait (Figure 3g), in contrast to the clockwise dynamics typically assumed by most baseline approaches, it is difficult for previous methods to capture this behavior (Figure 3h), whereas TSvelo can still faithfully model such patterns. For genes such as Rps26, which have critical roles in the development of blood progenitors to erythroid (Liu and Karlsson, 2024), the unspliced–spliced data is so noisy that cells of different types overlap in phase portrait. TSvelo can still capture the gene dynamics and reveals differences in transcription rates (Alpha) across cell types (Figure 3i). In contrast, methods that rely solely on unspliced–spliced data from individual genes, such as scVelo, fail to accurately fit these dynamics and cannot identify them as velocity genes, which further validates TSvelo’s performance in complex scenarios. The visualization of dynamics fitting on more genes is provided at Figure 3—figure supplement 2.

From the TF–target weight matrices learned by TSvelo, we can extract key regulatory relationships. First, we calculate the mean absolute weight of each TF across its interactions with target genes. The TFs with the highest mean absolute weights are ranked and presented in Figure 3j. Notably, Klf1, a critical TF in blood progenitors (Mukherjee and Bieker, 2022), is assigned the highest weight. We then examine the target genes of KLF1, highlighting those with the highest weights (Figure 3k), among which Hba-x (Perseu et al., 2011), Aals2 (Singleton et al., 2008), and Gypa Kuvardina et al., 2015 have been previously identified as key genes in erythroid differentiation and are known to be regulated by Klf1. A time-delay pattern is also observed between Klf1 and its target genes (Figure 3l), indicating that Klf1 expression increases first, followed by a corresponding upregulation of its target genes.

We apply TSvelo to a 10x multi-omics dataset from the embryonic mouse brain, which includes both Assay for Transposase-Accessible Chromatin with sequencing (ATAC-seq) (Buenrostro et al., 2013) and scRNA-seq data. This dataset was previously used in the Multivelo study, which introduced an RNA velocity model designed for multi-omics datasets to capture the dynamics between chromatin accessibility and unspliced mRNA (Riba et al., 2022). As introduced by a previous study (Li et al., 2023), Radial Glia (RG) cells located in the outer subventricular zone give rise to neurons, astrocytes, and oligodendrocytes. During cortical development, neurons follow an inside-out layering pattern in which earlier-born neurons populate the deep cortical layers, whereas later-born neurons migrate past them to occupy more superficial layers (Nadarajah and Parnavelas, 2002). Additionally, RG cells are capable of producing intermediate progenitor cells (IPCs), which serve as neural stem cells and generate a variety of mature excitatory neurons within the cortical layers.

Multivelo is chosen as the baseline in this study since it connects the regulation signal and RNA velocity in ATAC–unspliced–spliced space. TSvelo uses the same data processed by Multivelo for a fair comparison. The velocity stream obtained by TSvelo and Multivelo is shown in Figure 4a, b, respectively. TSvelo produces velocity patterns that appear more consistent with expected differentiation trends, especially for the RG (in red) to IPCs (in brown) process. The pseudotime inferred by TSvelo is also coherent with the whole cell development process (Figure 4c).

Figure 4 with 3 supplements see all

Download asset Open asset

We next further analyze the details in the phase portrait fitting. On Meis2 (Figure 4d), both Multivelo and TSvelo successfully model its dynamics, which follows a clear pattern on the unspliced–spliced phase portrait. However, due to the high noise and sparsity inherent in ATAC-seq data, Multivelo encounters difficulty on genes that show significant overlap between different cell types in the unspliced–spliced phase portrait. For example, in the case of Basp1 (Figure 4e), Multivelo incorrectly models their expression patterns as monotonically increasing, failing to capture the true dynamics, which involve initial upregulation followed by downregulation. In contrast, TSvelo’s prediction is consistent with known biological processes and identifies the delay from transcription to the unspliced mRNA and also the delay from unspliced to spliced mRNA (The rightmost plots in Figure 4e). Msi2, which has been reported to be highly expressed in neural stem/progenitor cells (Liu et al., 2019), is more accurately modeled by TSvelo, exhibiting a decreasing expression pattern during the differentiation process. By comparison, Multivelo fails to capture the correct trend at the initial stage (Figure 4f). We also show the learned transcriptional rate α as well as the predicted dynamics of unspliced and spliced RNA along pseudotime t, which clearly illustrate the time delays between transcription to unspliced and unspliced to spliced RNAs. The visualization of dynamics fitting on more genes is provided at Figure 4—figure supplement 1. Considering that ATAC- and TF RNA-based signals capture distinct aspects of gene regulation (Figure 4—figure supplement 2), integrating both modalities may be beneficial for future models that aim to more comprehensively characterize transcriptional regulation. Results on an additional mouse pons dataset are shown in Figure 4—figure supplement 3.

Given that multi-lineage differentiation is a common phenomenon in larger and larger scRNA-seq datasets, developing the RNA velocity model that can handle such complexity is essential. Because TSvelo demonstrates robust performance across various tasks, its applicability can be extended to more complex situations, such as scRNA-seq datasets where cells differentiate into multiple fates. We apply TSvelo to a multi-lineage dentate gyrus scRNA-seq dataset, which captures the differentiation process from neural blast cells to various cell types (Hochgerner et al., 2018). During preprocessing, the velocity genes selected by TSvelo are enriched in GO terms related to neural development, such as axonogenesis (GO:0007409), axon guidance (GO:0007411), and axon development (GO:0061564) (Figure 5a), which aligns with the biological processes described by the dataset.

Figure 5 with 2 supplements see all

Download asset Open asset

TSvelo could correctly identify three lineages from this data. The lineage segmentation process is fully automated and does not require prior knowledge about the presence of multiple branches (details are provided at Methods and Figure 5—figure supplement 1). After applying the TSvelo model to each lineage and combining the results, the inferred pseudotime and velocity stream plots (Figure 5b, c) align with the ground-truth differentiation trajectory well, which begins with neural blast cells. We also compare TSvelo to baseline methods, including scVelo and cellDancer. Notably, cellDancer is designed to handle such multi-lineage data by modeling cell- and gene-specific transcriptional rates, splicing rates, and degradation rates using neural networks. As shown in Figure 5d, e, scVelo struggles with this multi-lineage dataset, and cellDancer fails to correctly capture the trajectory of immature granule cells (colored in purple).

Next, we analyze gene expression patterns across different lineages. TSvelo provides inferred dynamics along each lineage, revealing that the expression of Ank3 in the granule and Cornu Ammonis (CA) lineages follows an increasing and then decreasing pattern. In contrast, the expression of Ank3 in the glial lineage decreases to zero (Figure 5f). The similar pattern is also observed for other genes, such as Map1b (Figure 5g). Both Ank3 (Leussis et al., 2012) and Map1b (Meixner et al., 2000) are associated with GO terms axonogenesis and axon guidance. This observation is consistent with the fact that the granule and CA (Douglas, 1967) predominantly consist of neurons, where axons are a critical component. In contrast, glial cells, such as astrocytes, typically lack axons, providing a potential explanation for these expression patterns. Slc1a2 exhibits a distinct expression pattern, showing a significant increase in the glial lineage and a decrease in the granule and CA lineages. This pattern aligns with the observation that astrocytes are the primary cell type expressing Slc1a2 (Sun et al., 2023). Details about dynamics fitting for genes on each lineage are provided at Figure 5—figure supplement 2.

The Lineage and RNA Recovery (LARRY) method utilizes barcoded hematopoietic cells to trace both cell lineage and gene expression over time. LARRY has been successfully employed to track the in vitro differentiation of human blood cells, accurately capturing lineage trajectories and cell fates (Weinreb et al., 2020).

On this LARRY dataset, which encompasses a total of 49,302 cells on multiple lineages, TSvelo effectively detects the initial Leiden clusters based on the unspliced-to-spliced delay (Figure 6a, b) and separates its lineages (details are provided at the Lineage segmentation and pseudotime initialization section in Methods). Furthermore, Figure 6c, d shows that the pseudotime and velocity stream inferred by TSvelo can capture the differentiation process, progressing from undifferentiated cells to distinct cell fates. The velocity genes selected by TSvelo are significantly enriched in processes related to neutrophil biology, such as neutrophil degranulation (GO:0043312), neutrophil activation in immune response (GO:0002283), and neutrophil-mediated immunity (GO:0002446).

Figure 6

Download asset Open asset

We further analyzed the expression patterns of genes associated with neutrophil degranulation, focusing specifically on the neutrophil lineage. Using four representative genes as examples (Figure 6f, g), TSvelo can help observe the significant variation in the expression patterns across these genes. Pygl, which has been widely reported as a key gene in neutrophil (Borella et al., 2022), is almost exclusively expressed in neutrophil cells and exhibits an increasing expression pattern during neutrophil differentiation. In contrast, Ms4a3 was examined across all lineages, and using TSvelo, we observed that its expression initially increases and then decreases in both the neutrophil and monocyte lineages. This pattern is consistent with prior studies identifying Ms4a3 as a gene specifically expressed by granulocyte–monocyte progenitors (Liu et al., 2019). Similarly, Clec12a, another gene associated with neutrophil degranulation, shows a comparable expression trajectory during neutrophil differentiation. Previous research has indicated that Clec12a expression is highest in granulocyte–macrophage progenitors (Bill et al., 2018). Lta4h exhibits a pattern similar to that of Clec12a, suggesting it may play a significant role in the early stages of neutrophil differentiation. These findings further confirm the utility of TSvelo for gene-level analysis in multi-lineage differentiation datasets.

The unspliced and spliced RNA abundances are preprocessed with multiple steps, which include highly variable genes (HVGs) selection, normalization, log transformation, kNN smoothing, and clustering. Following the data preprocessing procedures outlined in scVelo (Bergen et al., 2020), we first select the top 2000 HVGs and normalize their expression profiles by dividing by the total counts in each cell. A nearest-neighbor graph (with 30 neighbors by default) was constructed based on Euclidean distances in principal component analysis space (with 30 principal components by default) on log-transformed gene expression data. Subsequently, we compute the first- and second-order moments (mean and uncentered variance) for each cell across its nearest neighbors. These steps are performed by using scvelo.pp.filter_and_normalize() and scvelo.pp.moments().

To construct a prior gene network for the selected HVGs, we utilize TF–target annotations from the ENCODE (Feingold, 2004) and ChEA (Lachmann et al., 2010) databases. If a regulatory interaction between a TF and its target gene is identified in either of these two databases, the TF is included in the set of regulators for modeling the dynamic behavior of the target gene. During model optimization, the contribution of each prior regulatory edge is adaptively adjusted, allowing unsupported interactions to be effectively down-weighted. As a result, TSvelo is expected to be relatively robust to false-positive TF–target annotations. By contrast, missing true regulatory interactions are not represented in the candidate regulator set and therefore cannot contribute to the learned regulatory dynamics, making the model potentially more sensitive to false negatives in the prior network. Ablation studies on the selection of TF–target resources (Figure 1—figure supplements 4 and 5) verify that integrating multiple TF–target resources can improve performance by increasing regulatory coverage and reducing false negatives in the prior network.

Previous studies Bergen et al., 2020 have demonstrated that only a subset of genes, termed ‘velocity genes’, can be accurately fitted in the unspliced–spliced phase portrait by RNA velocity models. Inspired by the notion that velocity genes provide high-quality data for dynamic modeling in the phase portrait, we propose to select velocity genes during the preprocessing step. The velocity genes are selected based on the premise that they exhibit clear dynamics in the unspliced–spliced phase portrait, which is necessary for precise modeling of splicing dynamics. The selection is based on the similarity between cell–cell neighborhood relationships in UMAP space and those observed in gene-specific unspliced–spliced phase portraits. In detail, we first compute the cell–cell neighborhood graph, Graph, using scvelo.pp.neighbors() on all HVGs. Next, we construct an anndata object adata_g for each gene, which contains only the unspliced and spliced expression data of that gene. We then apply scvelo.pp.neighbors() to each adata_g to calculate the gene-specific neighborhood graph, denoted as Graph_g. Subsequently, we compute the similarity between Graph and each Graph_g. The top 100 genes with the highest similarity are selected as velocity genes. These genes are characterized by phase portraits that exhibit a structure similar to that of the UMAP space, thereby enhancing the separation of cells from different types. As a result, the splicing dynamics of these velocity genes are more likely to be captured by RNA velocity models.

Based on the normalized gene expression data, we perform Leiden clustering using scanpy.tl.leiden() with a low resolution (default resolution = 0.1). These Leiden clusters are utilized to identify lineages and determine the differentiation direction for each lineage. Subsequently, we apply PAGA (scanpy.tl.paga()) to assess the connectivity between Leiden groups.

By setting one Leiden group as the initial state, we infer pseudotime using diffusion pseudotime (DPT) with scanpy.tl.dpt(). To detect lineages, we start from the selected initial Leiden cluster and filter the PAGA graph by applying a threshold (default = 0.02). Edges with weights below this threshold are removed, leaving only the shortest path from the initial state to each group. The remaining paths are considered as the detected lineages.

Given a lineage, TSvelo could detect the orientation of differentiation along it, which is based on the fact that unspliced RNA precedes spliced RNA in the expression pattern (Ge et al., 2025). For each gene, TSvelo computes the Spearman correlation between unspliced and spliced abundance across different moving steps along the DPT. The time number of moving steps along unspliced to spliced expression, which maximizes the Spearman correlation, is considered the U-to-S delay (Figure 1a shows an example of U-to-S delay with gene EML5). The U-to-S delay for each lineage is calculated as the mean U-to-S delay for all genes on it. And the average U-to-S delay across all lineages provides the overall U-to-S delay for the initial Leiden cluster. Finally, the initial Leiden cluster with the highest U-to-S delay is considered the correct initiation. The DPT under this condition is used to initialize pseudotime for TSvelo. If multiple lineages are detected, TSvelo will model each lineage independently and subsequently merge them in the final step. A detailed explanation with an example illustrating how lineages are segmented is provided in Figure 5—figure supplement 1. Results of the initialization strategy on all datasets are shown in Figure 1—figure supplement 3 demonstrating that clusters with the highest U-to-S delay scores can be reliably identified as the initial cell states. We also provide a simulation study to verify this initialization strategy in Figure 1—figure supplement 2.

Suppose $u_g$ and $s_g$ are the abundance of unspliced and spliced RNA, and ${\alpha }_g\left(t\right),{\beta }_g$, and ${\gamma }_g$ are the transcription, splicing, and degradation rates, respectively. The dynamics of a velocity gene $g$ is modeled as

We assume that the gene and cell-specific transcriptional rate ${\alpha }_g\left(t\right)$ is influenced by the expression of TFs, while ${\beta }_g$ and ${\gamma }_g$ are gene-specific constant parameters. Considering the wide usage of linear models for gene relations in previous studies (Li et al., 2024a; Tran et al., 2022; Song et al., 2023; Liu et al., 2011; Wang et al., 2023), we model ${\alpha }_g\left(t\right)$ as

In order to include more TFs in the modeling, we reserve those TFs that are not selected as velocity genes. These TFs are excluded from the velocity genes selection because their unspliced–spliced expression does not provide sufficient information, primarily due to high noise in the unspliced RNA. Consequently, we incorporate these TFs, whose spliced abundance is denoted as $s`$, into the TSvelo by directly modeling the process between transcriptional signal to mature RNA,

As a result, we can get the ODE model with parameters matrix $A$:

$U=[u_1,u_2,\dots ,u_G{]}^T\in R^G$, $S=[s_1,s_2,\dots ,s_G{]}^T\in R^G$, $S=[s_{G+1},s_{G+2},\dots ,s_{G+K}{]}^T\in R^K$, $G$ is the number of velocity genes, and $K$ is the number of additional TFs which are not selected as velocity genes. ${\rm B}\in R^{G\times G}$ and $\Gamma \in R^{G\times G}$ are the diagonal matrices consisting of splicing rates $\left[{\beta }_1,{\beta }_2\dots ,{\beta }_G\right]$ and degradation rates $\left[{\gamma }_1,{\gamma }_2\dots ,{\gamma }_G\right]$ for velocity genes, respectively. $\Gamma `\in R^{K\times K}$ is the diagonal matrix consisting of degradation rates $\left[{\gamma }_{G+1},{\gamma }_{G+2}\dots ,{\gamma }_{G+K}\right]$ for these additional TFs. $W\in R^{G\times G}$ and $W^{`}\in R^{K\times K}$ denote the matrices representing the TF–target relationships for the TFs selected as velocity genes and those not selected as velocity genes, respectively. Using these gene–gene weight matrices, TSvelo can model the gene- and cell-specific transcriptional rate ${\alpha }_g$. The whole parameters matrix $A\in R^{(2G+K)\times (2G+K)}$.

In TSvelo, the dynamics described in Equation 7 are implemented using a Neural ODE model. Given an initial state and the parameters in matrix $A$, the Neural ODE model can compute the values of U and S at any time step. This approach eliminates the need for an analytical solution for unspliced and spliced expression, enhancing the flexibility of TSvelo. By default, the number of time steps in the ODE is set to 1000. The ReLU activation function is applied to the parameters ${\alpha }_g$, ${\beta }_g$, and ${\gamma }_g$, in order to prevent them from taking negative values. Notably, TSvelo does not incorporate deep neural networks or encoders. Since TSvelo models all genes simultaneously, we normalize the unspliced and spliced abundances of each gene by its standard deviation before inputting the data into the Neural ODE model to ensure that the influence of each gene is balanced.

TSvelo optimizes the pseudotime $t$ and parameters matrix $A$ in ODE iteratively using an EM approach. The maximum number of iterations is set to 30 by default, with an early stopping criterion applied. Given the parameters in the ODE system, we can compute the predicted values of unspliced (for velocity genes) and spliced (for velocity genes and additional TFs) at all time steps, denoted as $\hat{U}\left(t_{step};A\right)$, $\hat{S}\left(t_{step};A\right)$, and $\widehat{S^{\prime }}\left(t_{step};A\right)$, where $t_{step}\in \left\{\mathrm{0,1},2\dots ,999\right\}$. Given the time steps assignment $t_c$ for each cell $c$, we can get the expected values for velocity genes at all cells as $\hat{U}\left(t_c;A\right),\hat{S}\left(t_c;A\right)$, and $\widehat{S^{\prime }}\left(t_c;A\right)$. The loss function used in the EM algorithm is the mean squared error between the expected value and the observed data.

The EM algorithm will iteratively update the time assignment and parameters in the Neural ODE by minimizing the above loss function.

In the E-step, following the strategy adopted in the scVelo model (Bergen et al., 2020), TSvelo updates the time step assigned to each cell using grid search with time steps range $t_{step}\in \{\mathrm{0,1},2\dots ,999\}$, which is achieved by finding the minimum distance between the observed unspliced–spliced expression and the model prediction ($\hat{U}$ and $\hat{S}$) at all time steps.

In the M-step, the parameter matrix $A$ is updated using gradient descent within the Neural ODE model to minimize the loss function. Here we use the prior gene regulation knowledge to constrain the weight matrix $W$ and $W^{`}$ in $A$. Only if gene $i$ is annotated as a TF for gene $j$ in the ENCODE or ChEA databases, the corresponding $w_{ij}$ could be learnable. Otherwise, $w_{ij}$ is fixed at zero. The usage of prior knowledge could keep the matrix sparse and avoid overfitting. The Neural ODE module is implemented using the torchdiffeq package and trained with the Adam optimizer, a learning rate of 0.05, and a maximum of 500 epochs, incorporating an early stopping criterion to prevent overfitting.

By performing the EM optimization, TSvelo could fit the high-dimensional unspliced–spliced data across multiple genes well, and learn the global pseudotime for each cell. After the training procedure, the RNA velocity could be calculated using those parameters in matrix $A,$ which could be further used for cell fate prediction.

(1) Velocity consistency (VCon). We used the scvelo.velocity_confidence() function from scVelo to evaluate velocity consistency, interpreting the results as a measure of how consistent velocities are within neighboring cells. Velocity consistency is especially suitable for evaluating the RNA velocity modeling on single lineage. For each cell $c$, the velocity consistency is calculated as follows:

where $N\left(c\right)$ represents the neighboring cells of a given cell $c$. $v_c$ and $v_{c^{`}}$ denote the low-dimensional velocity vectors of cell $c$ and its neighboring cell $c^{`}$.

(2) Cross-boundary direction correctness (CBDir). Cross-boundary direction correctness is initially introduced by VeloAE (Qiao and Huang, 2021), which assesses the accuracy of transitions from a source cluster to a target cluster by examining the boundary cells, and requires ground-truth annotations. We directly run the function unitvelo.evaluate() provided in UniTVelo Gao et al., 2022 to obtain the Cross-boundary direction correctness. In detail, the CBDir is calculated as follows:

where $C_A$ denotes the set of cells in the target cluster A, and $N\left(c\right)$ represents the neighboring cells of a given cell $c$. ${\mathit{v}}_c$ and ${\mathit{x}}_c$ denote the low-dimensional velocity and state vectors of cell $c$, respectively, and ${\mathit{x}}_{c^{\prime }}$ denotes the state vector of its neighboring cell.

(3) Within-cluster velocity coherence (ICCoh). Within-cluster velocity coherence is initially introduced by VeloAE (Qiao and Huang, 2021), which measures the coherence of velocities within a single cluster using a cosine similarity score between cell velocities. We applied the function unitvelo.evaluate() provided by UniTVelo (Gao et al., 2022) to directly compute the within-cluster velocity coherence. Using the same notation as defined above, the CBDir is calculated as follows:

To quantitatively assess whether TSvelo can distinguish cell types, we evaluated the separability of cell-type labels in both the 2D (unspliced–spliced) phase portrait adopted by previous RNA velocity approaches, and the 3D (α–unspliced–spliced, α denotes the transcriptional rate) phase portrait introduced by TSvelo.

Specifically, we evaluated how well the embedding preserves cell-type information using a kNN classification accuracy with fivefold cross-validation. Given an embedding matrix in 2D or 3D space ($X\in {\mathbb{R}}^{n\ast d}$, where $n$ is the number of cells and $d$ is 2 or 3) and corresponding cell-type labels ($y\in \{1,\dots ,C\}$), we partition the data into fivefolds. For each fold ($k$), a kNN classifier with $K=5$, denoted as $f^{\left(k\right)}$, is trained on the training subset and evaluated on the held-out test subset. The classification accuracy for the $k$ th fold is defined as

where $n_k$ is the number of samples in the test set and 1(·) is the indicator function. The final score is obtained by averaging across all folds:

This metric directly assesses whether cells of the same type are positioned close to each other in the embedding space and is widely used to quantify representation quality.

GO terms enrichment analysis is a widely used method for identifying biological processes, molecular functions, or cellular components that are overrepresented in a given set of genes compared to a background set. Here we use the biological processes enrichment analysis in the GO terms database to explore the associated biological roles of those selected velocity genes. To perform the GO term enrichment analysis, we utilized the ‘gseapy.enrichr()’ function in Python using Fisher’s exact test by default.

For datasets containing multiple lineages, TSvelo applies its model independently to each lineage. The resulting lineage-specific outputs are integrated into a unified representation by aggregating both expression layers and cell-level annotations across branches.

All objects are aligned to a common gene space by restricting to the shared set of genes. Let $b=1,\dots ,B$ index the lineage-specific results. For each cell $c$ and each inferred quantity (e.g., velocity and pseudotime), TSvelo collects the corresponding results $x_{c,b}$ from all lineages in which the cell is present. Layer values are merged on a per-cell basis using a weighted average, where the weight for each branch is given by its number of cells, $n_b$. The merged value for cell $c$ is computed as

Continuous variables are averaged, while discrete variables are cast to integer values after aggregation.

Finally, the aggregated layer values and annotations are assigned back according to the global cell index, producing a unified dataset that integrates lineage-specific inferences.

TSvelo is implemented in Python. The source code can be downloaded from the GitHub repository, https://github.com/lijc0804/TSvelo (copy archived at Li, 2026).

We generated synthetic single-cell RNA velocity datasets using a mechanistic transcriptional dynamics model with one or multiple developmental branches. The system included 200 genes, among which 30 were designated as transcription factors (TFs).

For each branch, we independently sampled a TF–target regulatory matrix $W\in {\mathbb{R}}^{30\times 200}$ from a standard normal distribution to simulate distinct GRN structures. Gene expression dynamics were modeled using a coupled ordinary differential equation (ODE) system describing unspliced and spliced RNA abundances:

where $u$ and $s$ denote unspliced and spliced RNA levels, respectively. The transcription rate α was computed as a nonlinear function of TF expression, defined as a weighted sum of spliced TF abundance, followed by clipping to ensure bounded activation.

Each branch is initialized from the same randomly sampled initial condition drawn from a gamma distribution, allowing controlled divergence of trajectories driven solely by branch-specific regulatory programs.

To simulate observed sequencing counts, we introduced technical noise by scaling latent expression levels with cell-specific library sizes drawn from a log-normal distribution. The resulting expression counts were generated using a negative binomial sampling model:

where $\theta$ controls overdispersion, with smaller values corresponding to higher noise levels. The final datasets consist of paired unspliced (U) and spliced (S) count matrices with realistic transcriptional stochasticity and branching gene regulatory dynamics. For each branch, cells were further divided into three developmental stages for downstream analysis.

We evaluated TSvelo on multiple simulated datasets with varying numbers of branches and noise levels. There are two or three branches starting from the same root cell groups in these datasets (Branch 1: stage 0–stage 1–stage 2. Branch 2: stage 0–stage 3–stage 4. Branch 3: stage 0–stage 5–stage 6). The results of initial state identification based on the unspliced-to-spliced (U-to-S) delay, along with the corresponding 2D velocity stream visualizations, are presented in Figure 1—figure supplement 1. These results demonstrate that the U-to-S delay-based initialization is robust and consistently identifies cells corresponding to the earliest developmental stage (stage 0) across different simulation settings.

We also evaluated TSvelo and those splicing-based RNA velocity approaches on multiple simulated datasets with varying numbers of branches and noise levels. There are one, two, or three branches starting from the same cell group in these datasets (Branch 1: stage 0–stage 1–stage 2. Branch 2: stage 0–stage 3–stage 4. Branch 3: stage 0–stage 5–stage 6). We primarily assessed performance using the cross-boundary direction correctness (CBDir) metric, as it directly evaluates inferred trajectories against ground-truth cell-stage annotations. As shown in Figure 1—figure supplement 2, TSvelo consistently achieves the highest accuracy across all simulation settings, particularly in scenarios with complex branching structures, which pose significant challenges for baseline methods.

Along a correctly oriented trajectory, unspliced (U) expression is expected to precede spliced (S) expression due to transcriptional dynamics. Ideally, this U-to-S delay would be observable at the level of individual genes. However, due to the high noise inherent in scRNA-seq data, such delays are often not consistently detectable on a per-gene basis. To address this, we aggregate U-to-S delay signals across all genes and determine the lineage orientation by maximizing a global delay score. Under this criterion, the cluster from which all outgoing lineages exhibit the highest aggregated U-to-S delay is inferred to correspond to the initial state. The results on all datasets (Figure 1—figure supplement 3) suggest that the highest U-to-S delay scores can be used to detect the initial cluster.

we additionally incorporated the DoRothEA regulon database as an alternative prior with confidence-level filtering. We further performed ablation studies on the pancreas dataset and the gastrulation erythroid dataset using different TF–target resources, including ChEA, ENCODE, and their combinations with DoRothEA.

The results on the pancreas dataset and the gastrulation erythroid dataset are shown in Figure 1—figure supplements 4 and 5, respectively, which come up with the same conclusion. We observed highly consistent results across most TF–target prior combinations, including ChEA, ENCODE, ChEA + ENCODE, ChEA + DoRothEA, ENCODE + DoRothEA, and ChEA + ENCODE + DoRothEA. Using the pancreas dataset as an example, the mean velocity consistency ranged from 0.985 to 0.995, the mean in-cluster coherence ranged from 0.983 to 0.992, and the mean cross-boundary direction correctness ranged from 0.719 to 0.740 across all settings. These consistently high and tightly bounded metrics indicate that TSvelo is largely insensitive to the specific choice of TF–target prior.

The only configuration showing reduced stability was the use of DoRothEA alone, particularly in terms of cross-boundary direction correctness. This is likely due to its comparatively limited coverage of TF–target interactions. For instance, in the pancreas dataset, only 81 out of 2000 highly variable genes (HVGs) could be associated with TFs based on DoRothEA, corresponding to 102 TF–target links in total, which may restrict downstream regulatory modeling. In contrast, ChEA covered 1793 genes with 13,976 TF–target links, and ENCODE covered 1854 genes with 33,076 links.

These results further suggest that integrating multiple TF–target resources can improve performance by increasing regulatory coverage and reducing false negatives in the prior network. Given that unsupported edges can be down-weighted during training, TSvelo appears to be more sensitive to missing true regulatory interactions than to the inclusion of spurious ones.

we have added a systematic comparison of runtime and GPU memory usage across TSvelo and ba methods using simulated datasets of increasing scale (600, 1200, and 1800 cells) on our NVIDIA GeForce RTX 3090 device with 24 GB memory.

Appendix 1—table 2 shows differences in computational efficiency and resource requirements among methods. Specifically, classical methods such as scVelo and Dynamo exhibit very fast runtimes (10–24 s) and do not rely on GPU acceleration, reflecting their relatively lightweight modeling strategies. In contrast, deep learning-based approaches, including UniTVelo, cellDancer, and TSvelo, have higher computational costs due to their increased model complexity.

TSvelo exhibits a stable GPU memory footprint (~1.26 GB) across different dataset sizes, indicating that its memory usage is primarily determined by model architecture rather than the number of cells. This level of memory consumption is well within the capacity of modern GPUs and does not pose practical limitations. In terms of runtime, TSvelo scales approximately linearly with dataset size. The higher computational cost of TSvelo is mainly due to its EM-style optimization procedure, where each M-step also involves multiple optimization updates to infer gene regulatory effects in a global model. This design enables TSvelo to explicitly incorporate regulatory priors and jointly model gene interactions, which is not supported by these baseline methods.

To further improve runtime efficiency, TSvelo allows flexible control of the number of EM iterations. As shown in Figure 1—figure supplement 6 and Appendix 1—table 3, we evaluated performance under different iteration settings on the simulation dataset. The early stopping strategy employed in the EM framework of TSvelo, which will stop modeling if the loss is not further reduced in the last three iterations. Results show that convergence is typically achieved within three iterations for this dataset, and increasing the maximum number of iterations beyond this does not further change the results. Notably, even a single iteration already yields competitive performance, likely benefiting from the strong initialization based on unspliced-to-spliced temporal delay.

Overall, these results highlight a trade-off between computational efficiency and modeling expressiveness. While TSvelo is more computationally demanding than classical approaches, it provides a more flexible framework for incorporating regulatory information and capturing complex gene interactions, which we believe justifies the additional computational cost in scenarios requiring accurate dynamical inference.

We have added null simulations to evaluate TSvelo and baseline approaches on data without underlying dynamic structure. Specifically, we generated a null dataset including 200 genes and 600 cells by independently sampling spliced (S) and unspliced (U) counts, thereby removing any coherent transcriptional relationship between them.

When applying scVelo and UniTVelo to this data, no genes passed the velocity gene selection step under the default likelihood-based filtering, and no velocity field could be obtained. We further tested TSvelo, Dynamo, and cellDancer on the same null data and observed that all three methods still produce trajectory-like patterns despite the absence of true dynamics (Figure 1—figure supplement 7).

Including TSvelo, many RNA velocity and trajectory inference approaches assume that they are applied to datasets reflecting underlying dynamic biological processes. Incorporating additional checks during preprocessing could help prevent applying velocity analysis to non-dynamic datasets.

We conducted comparisons across all methods on the pancreas dataset, with quantitative evaluations shown in Figure 2—figure supplement 1. In each plot, methods are ranked in descending order of their mean values. Numbers at the bottom indicate the sample size for each metric. Statistical significance is assessed using a one-sided Mann–Whitney U test, where *****, ***, **, and * denote p < 0.00001, 0.0001 ≤ p < 0.001, 0.001 ≤ p < 0.01, and 0.01 ≤ p < 0.05, respectively.

TSvelo significantly outperforms all baseline methods in terms of velocity consistency. For in-cluster coherence, TSvelo performs comparably to the best-performing baselines (UniTVelo and TFvelo) and significantly outperforms several competing methods, including CellDancer, Dynamo, and scVelo. For cross-boundary direction correctness, TSvelo shows consistent improvements in mean performance. The pairwise comparisons on cross-boundary direction correctness do not reach statistical significance. It is likely influenced by the limited number of independent samples (n = 7), which reduces statistical power for detecting differences. Importantly, TSvelo still achieves the best average performance among all methods, indicating a consistent overall trend in favor of TSvelo.

TSvelo is capable of accurately modeling the dynamics of spliced mRNA abundance for those TFs which are not selected as velocity genes during preprocessing. These TFs cannot be well modeled in the unspliced–spliced phase portrait. However, TSvelo overcomes this limitation by directly modeling the process from transcription to spliced mRNA. Figure 2—figure supplement 3 shows the modeling on those TFs on pancreas dataset, where TSvelo directly models the dynamic between transcription and spliced abundance.

In the conventional 2D u–s phase portrait, cells from different transcriptional states may overlap, leading to reduced separability. In contrast, introducing the latent variable α expands the representation to a 3D space, which helps disentangle these mixed states and yields a clearer phase structure. We provide quantitative evidence on this gastrulation erythroid dataset in Figure 3—figure supplement 1, showing that the 3D representation achieves consistently higher kNN classification accuracy for cell-state separation compared to the 2D u–s embedding (one-sided Mann–Whitney U test, p-value = 0.002).

Figure 3—figure supplement 2 shows the dynamics fitting on additional genes with TSvelo on the gastrulation erythroid dataset.

Figure 4—figure supplement 1 shows the dynamics fitting on additional genes with TSvelo on the mouse brain dataset.

We have conducted the requested analysis by computing gene-wise chrome accessibility rate used in MultiVelo and the learned transcription rate from TSvelo, and evaluated their correlation across genes. As shown in Figure 4—figure supplement 2, the two estimates exhibit almost no global correlation across genes, indicating that they capture substantially different aspects of regulatory information.

This discrepancy is not unexpected and reflects the fundamental differences between these modalities. scATAC-seq measures chromatin accessibility, which provides a proxy for cis-regulatory potential of genomic regions. However, ATAC signals are inherently sparse and often exhibit a near-binary structure, limiting their ability to directly capture fine-grained temporal regulatory dynamics. In contrast, TF RNA expression reflects downstream transcriptional output, which is shaped by multiple regulatory layers, including post-transcriptional regulation, protein activity, temporal delays, and indirect regulation through intermediate transcriptional or signaling pathways. As a result, these two modalities are expected to capture complementary but not directly comparable aspects of gene regulation.

Overall, this result suggests that ATAC- and TF RNA-based signals capture distinct aspects of gene regulation. This further implies that integrating both modalities may be beneficial for future models that aim to more comprehensively characterize transcriptional regulation.

The pons dataset originates from the hindbrain of adolescent mice and captures the differentiation pathway from oligodendrocyte precursor cells (OPCs) to committed oligodendrocyte precursor cells (COPs), progressing to newly formed oligodendrocytes, and finally culminating in myelin-forming oligodendrocytes. The velocity genes selected through TSvelo’s preprocessing strategies are mostly enriched in Gene Ontology (GO) terms related to the development of the nervous system, including regulation of axon extension involved in axon guidance (GO:0048841), axonogenesis (GO:0007409), and nervous system development (GO:0007399) (Figure 4—figure supplement 3a).

Using the unspliced-to-spliced delay, TSvelo can correctly detect the initial Leiden cluster on the pons dataset (Figure 4—figure supplement 3b, c). After performing the model optimization for the neural ODE system, the pseudotime and velocity learned by TSvelo can accurately fit the differentiation process in the data (Figure 4—figure supplement 3d, e). Figure 4—figure supplement 3f illustrates the dynamics of multiple genes modeled by TSvelo. A challenge in this dataset is the abrupt shift in gene expression between OPCs (red) and COPs (orange). For example, as shown in Figure 4—figure supplement 3f for CHN2 and DPYSL2, there is a significant gap in expression levels between OPCs and COPs. Despite this, TSvelo successfully captures their dynamic patterns.

Using the dentate gyrus dataset, we demonstrate how lineages are segmented during preprocessing, and subsequently combined during postprocessing after being modeled individually. Notably, the lineage segmentation process is fully automated and does not require prior knowledge about the presence of multiple branches.

First, Leiden clustering is performed on the preprocessed scRNA-seq data, with the resolution set to 0.1 by default (Figure 5—figure supplement 1a). A cluster–cluster graph is then constructed using PAGA (Figure 5—figure supplement 1b), followed by a filtering strategy that removes edges with weights below a threshold of 0.02 by default. Next, each Leiden cluster is treated as an initial cluster, and the appropriate initialization is determined using the unspliced to spliced delay.

For instance, using Leiden cluster 3 (colored in red at Figure 5—figure supplement 1a) as the initial cluster, diffusion pseudotime (DPT) is initialized with the selected cluster (Figure 5—figure supplement 1c). Starting from this initial cluster, the shortest paths to all other clusters are computed. Any path that is a subset of another is discarded. The remaining paths correspond to the detected lineages. In the dentate gyrus dataset, initializing with Leiden cluster 3 yields three distinct lineages (Figure 5—figure supplement 1d). By applying DPT to each lineage, we can calculate the U-to-S delay for each lineage and determine the overall U-to-S delay for this initialization with the selected cluster.

After calculating the U-to-S delay for each Leiden cluster initialization, the best initialization cluster is chosen based on the highest U-to-S delay (Figure 5—figure supplement 1e). The corresponding DPT is then used to initialize pseudotime in the downstream neural ODE model. Subsequently, each lineage is processed independently in the TSvelo for optimization, yielding both the optimized pseudotime and velocity model. The results for each lineage are shown in Figure 5—figure supplement 1f.

Finally, on the dataset where multiple branches are detected, to combine the models for all lineages and obtain the overall pseudotime and velocity stream, we compute the average velocity (Figure 5—figure supplement 1g) and pseudotime (Figure 5—figure supplement 1h) across the different lineages. For example, in the dentate gyrus dataset, since Leiden clusters 0, 1, 2, and 4 each belong to a single lineage, the velocity and pseudotime values calculated for these lineages are directly assigned to the combined file. For Leiden cluster 3, which is included in all three branches, the mean velocity and mean pseudotime across the three lineages are calculated for cells in this cluster.

Figure 5—figure supplement 2 shows the detailed dynamics fitting on each lineage for Ank3, Map1b, and Slc1a2 on the dentate gyrus dataset.

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

This work was supported by the National Key R&D Program of China (2023YFF1204500 to YY) and the National Natural Science Foundation of China (No. 62503452 to JL).

1. Preprint posted: August 19, 2025
2. Sent for peer review: September 3, 2025
3. Reviewed Preprint version 1: November 25, 2025
4. Reviewed Preprint version 2: June 22, 2026
5. Reviewed Preprint version 3: July 27, 2026
6. Version of Record published: September 15, 2026

You can cite all versions using the DOI https://doi.org/10.7554/eLife.108950. This DOI represents all versions, and will always resolve to the latest one.

© 2025, Li et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.