Main
Recent single-cell technologies can measure multiple genome-wide readouts from the same individual cells, often pairing the transcriptome with epigenomic measurements such as chromatin accessibility, histone modifications, transcription factor (TF) occupancy or the three-dimensional (3D) chromatin interactome1,2,3,4,5,6,7,8,9,10,11,12,13,14, at a single snapshot of time. However, cells are intrinsically dynamic systems where complex regulatory mechanisms govern the gene expression and epigenetic programs that determine cell fate over time. As real-time in vivo measurements of these modalities are generally not feasible, computational methods have been developed to infer cell dynamics from single-snapshot datasets. Here, cells from the same time point are viewed as observations sampled from a dynamically evolving system, and each cell is assigned an inferred pseudotime or latent time relative to this presumed dynamic process. One prominent idea to emerge from these efforts was RNA velocity, initially proposed by Velocyto15 and later enhanced by scVelo16, which models spliced and unspliced reads from scRNA-seq using ordinary differential equations (ODEs) to estimate the rate of change of gene expression in an implicit dynamic process. RNA velocities are useful for providing pseudotemporal directional information in inferring gene expression trajectories, although velocities are only computable for a subset of genes and may lack robustness17,18. With single-cell multiome data, providing the transcriptome and chromatin accessibility in the same cells, single-cell TF binding motif accessibilities derived from scATAC-seq can serve as a proxy for TF regulatory activities. Therefore, TF motif accessibilities and their underlying velocities should provide information beyond RNA expression and RNA velocities to improve modeling of dynamic cell state changes.
Here, we present a deep learning model called DynaVelo for inferring cell dynamics from single-cell multiome data that exploits RNA and TF motif accessibility modalities and predicts shifts in cell fate trajectories under genetic perturbations. Our main hypothesis is that these two modalities influence each other in dynamic cellular systems; the accessibility of TF motifs, as proxies for TF regulatory activity, affects the expression of target genes, while the RNA expression of TFs, cofactors and epigenetic regulators modulates downstream changes in chromatin state. DynaVelo trains a generative neural ODE to jointly model the dynamics of gene expression and TF motif accessibility at a single-cell level. A previous method called MultiVelo19 also uses both scRNA and scATAC to estimate the RNA velocities of cells, extending the gene-based ODE models of scVelo (spliced and unspliced counts) to incorporate a third ODE for gene-level chromatin accessibility, but does not provide a functional form for the underlying vector field. Therefore, MultiVelo allows us to draw trajectories but not to intervene in the system to ask how perturbations would alter the dynamics or trajectories. Recent deep models for gene expression dynamics include Dynamo20, which infers absolute RNA velocity, predicts cell fates and perturbation outcomes and extracts regulatory relationships using differential geometry. However, training this model requires single-cell RNA-seq paired with metabolic labeling, which is not feasible in many experimental systems. scDiffEq21 uses drift-diffusion modeling of cell dynamics but only uses RNA expression and therefore fails to leverage scATAC-derived TF activity dynamics. DynaVelo thus moves beyond scRNA-based and multiome-based RNA velocity methods (scVelo and MultiVelo) and neural network methods (Dynamo and scDiffEq) that model the dynamics of gene expression alone.
We reasoned that DynaVelo’s use of single-cell TF motif activity data would be most relevant in cellular systems with rapid and complex dynamics where cell state transitions are driven by TFs. Furthermore, to be able to validate our findings, we needed a setting where loss of function of epigenetic regulators is known to alter TF activities and cell fate decisions. Therefore, we applied DynaVelo to murine germinal center (GC) B cells representing three genotypes: the wild-type (WT) setting and GC B cells bearing heterozygous loss-of-function mutations in Arid1a (Arid1aHet), a subunit of the canonical BAF chromatin remodeling complex that is commonly somatically altered in human B cell lymphomas22,23, and Ctcf (CtcfHet), a key factor involved in 3D chromatin organization24,25. We found that DynaVelo can recapitulate the complex in vivo dynamics of the GC reaction, generate cells and trajectories de novo from the model and recover evolving TF activities and a dynamic gene regulatory network (GRN) over latent time. DynaVelo can both predict the impact of loss-of-function genetic perturbations and identify TF perturbations that rescue loss-of-function dynamics and immunological phenotypes. Thus, DynaVelo provides a generative artificial intelligence framework for modeling and perturbing dynamic cellular systems from single-snapshot single-cell multiome data.
Results
DynaVelo learns a latent neural ODE model for single-cell dynamics from multiomic data
DynaVelo is a neural ODE framework for learning cellular dynamics from single-cell multiome data (scRNA + ATAC-seq). As TF motif accessibility is an interpretable cell-level readout derived from scATAC, we first used chromVAR26 to obtain a cell-by-TF-activity matrix. We also obtained a cell-by-gene matrix from the scRNA readout after preprocessing, normalization and filtering out unused genes (Methods). DynaVelo models the joint dynamics of RNA expression and TF motif accessibility in each cell. Because solving these ODEs in the original high-dimensional input spaces is challenging, we embed the cells, using a variational encoder on both the expression and motif accessibility inputs, and learn the neural ODE in a lower-dimensional latent space (Fig. 1a). Single-cell states and velocities learned in the latent space define a vector field over cells, which are mapped back to the observational space by a decoder. DynaVelo also learns latent times for all cells, providing a cell ordering in the joint vector field. Similar to variational autoencoders, DynaVelo uses a loss function to reconstruct the RNA and TF motif matrices and ensure that learned initial (latent time t = 0) latent states and cell-specific latent times do not greatly diverge from their specified prior distributions (Methods). Additional components to the loss function (1) enforce the similarity between learned RNA velocities and those provided for a subset of genes by a method such as scVelo and (2) ensure consistency between RNA and TF motif velocities (Methods).
a, The cell-by-gene gene expression and cell-by-TF motif accessibility matrices are first embedded into a lower-dimensional latent space by encoders. The initial latent state (z(0)) and latent time (t) of each cell are predicted in the latent space by neural networks and passed to a neural ODE model, which solves the ODE from time (t=0) to time (t) and starting from (z(0)). After receiving the latent state (z(t)) and velocity ({v}_{z}(t)), which forms the latent vector field, two decoders map the RNA and motif vector fields to the original RNA and motif space, reconstructing the RNA expression and motif accessibility and predicting the RNA and motif velocities. b, DynaVelo can learn dynamic GRNs using two approaches: Jacobian matrices and in silico perturbations. Dynamic GRNs, unlike static ones, capture the cell-state-specific gene regulatory relationship by measuring velocity changes after RNA expression or motif accessibility changes. Jacobians measure derivatives of changes, while in silico perturbations measure finite changes. c, In silico perturbation of WT samples is used to validate the models by comparing the change in velocities with the velocities of KO samples. Similarly, in silico perturbation of KO samples is used to identify candidate perturbation targets to rescue loss-of-function mutations.
After training (Methods), the model can be used for multiple downstream tasks. DynaVelo can infer dynamic and cell-state-specific GRNs using two different approaches. First, we can estimate how a small change in the expression of one gene locally changes the RNA velocity of another gene by computing Jacobian matrices at individual cells (Methods). Second, we can perform in silico perturbation (knockout, KO) of RNA expression and TF motif accessibility values and calculate the resulting change in RNA and motif velocities. When perturbing TFs, we set both the TF’s RNA and motif accessibility to their minimum possible values across cells; for non-TF genes, we only set the RNA to the minimum value (Fig. 1b). Both approaches yield dynamic GRNs by capturing which individual gene perturbations lead to RNA velocity changes of other genes, conditioned on a specific cell state. We can also study the effects of gene perturbations on TF motif velocities, for example, to reveal how epigenetic regulators alter motif accessibility landscapes across cells (Fig. 1b).
As DynaVelo provides a functional form of the learned vector field, the model can be used with a WT single-cell dataset to predict the global effect of (germline or somatic) loss-of-function alterations in disease or, given cells bearing such mutations, identify the best perturbation targets to rescue the loss-of-function phenotypes and restore homeostatic dynamics (Fig. 1c). For this second task, we perform in silico perturbation of genes in the mutant system and see which ones shift the RNA and motif velocities toward WT values (Fig. 1c).
DynaVelo learns complex dynamics and latent times of GC B cells
We first applied DynaVelo to a single-cell multiome dataset profiling WT murine GC B cells (n = 7,166 cells). To guide the dynamics, we used RNA velocity estimates from scVelo for the subset of genes with sufficient coverage, called velocity genes. DynaVelo then learned RNA velocities of all genes, as well as TF motif velocities, which were visualized by projection to the corresponding RNA expression and TF motif accessibility uniform manifold approximation and projection (UMAP) spaces (Fig. 2a). The clear branching of these vector fields into two types of trajectories recapitulated the main well-known dynamics of GC B cells. In one branch, highly proliferative centroblasts (CBs; shown in green) progress through the cell cycle in the GC dark zone. In the second branch, CB cells transition into the centrocytes (CC; shown in blue) in the light zone. A third important and more difficult trajectory to capture is the recirculation of some CCs from the light zone back into the dark zone, visible as flow lines that pass through the recycling CC population (shown in orange). Recycling CC cells are presented with antigen by follicular dendritic cells and receive survival signals from T follicular helper cells, before migrating back to the dark zone to undergo somatic hypermutation and improve their B cell receptor (BCR) affinity27,28,29. Lastly, some CCs exit the GC reaction toward prememory B cells (shown in purple) or plasmablasts (shown in cyan; Fig. 2a), to eventually differentiate into memory B cells and plasma cells. Figure 2b summarizes the major GC dynamics and trajectories together with cell type annotations.
a, RNA and motif velocities projected to the UMAP embeddings of RNA expression and motif accessibility, respectively, showing RNA and motif vector fields. The most important cell state transitions are uncovered: CB → CC, CB cell cycle, CC → prememory, CC → plasmablast and CC → CB migration. b, The summarized dynamics of GC B cells. c, Latent times of the cells learned by DynaVelo and scVelo shown on both RNA and motif UMAPs. DynaVelo latent times better match the dynamics of GC B cells and identify more biologically reasonable initial cells (that is, cells with latent time t close to 0). d, DynaVelo is used to simulate trajectories capturing the dynamics of two main types: CB cell cycle and CB → CC transition. The generated RNA and motif trajectories along these simulated cells are embedded along with the input cells into the RNA and motif UMAP spaces. e. Dynamics of two critical GC TFs, Spi1 and Nfkb1 in two simulated trajectories, cell cycle and CB → CC. Solid lines and shades show the mean and its 95% confidence interval based on standard errors on 20 generated trajectories. Trends of four modalities are shown for each TF: RNA expression, RNA velocity, motif accessibility and motif velocity. Velocity curves are the time derivatives of the state curves.
We next plotted DynaVelo’s learned latent times, shown in both RNA expression and motif accessibility UMAP spaces and compared to scVelo’s latent times (Fig. 2c). This side-by-side comparison showed that DynaVelo’s latent times better recapitulated GC dynamics than those of scVelo. While both methods identified prememory B cells and plasmablasts as terminal states, scVelo failed to identify the major bifurcation between cells that cycle in the dark zone and those that transition to the light zone (Fig. 2c). This discrepancy was even more evident in the motif space, where only DynaVelo’s latent times provided a meaningful pseudo-ordering (Fig. 2c). These results show that by exploiting RNA expression and TF motif accessibility, DynaVelo outperformed scVelo at capturing latent times in complex cell state dynamics within GCs.
We also compared DynaVelo’s RNA velocities and latent times to MultiVelo, a multiome-based RNA velocity method (Supplementary Fig. 1). We observed that MultiVelo’s RNA velocities did not capture the dominant transition of CBs to CCs and instead had all trajectories moving from CC to CB (Supplementary Fig. 1a). Moreover, unlike DynaVelo, MultiVelo failed to identify the plasmablasts as one of the GC exits and instead incorrectly identified them as a starting population (Supplementary Fig. 1a). This mistake was also evident in MultiVelo’s latent times, positioning plasmablasts as the starting point (Supplementary Fig. 1b). More comparisons of DynaVelo to scVelo and MultiVelo (Supplementary Fig. 2), the robustness and reproducibility analysis (Supplementary Fig. 3), generalization analysis of DynaVelo across samples (Supplementary Fig. 4) and RNA-only DynaVelo predictions (Supplementary Fig. 5) are provided in the Supplementary Note.
As a generative model, DynaVelo can generate trajectories starting from or ending at any cell and generate cells de novo along these trajectories, enabling analysis of cell state transitions along latent time and determination of possible cell fates given a starting cell state. We used DynaVelo to generate two example trajectories starting in the CBs, one progressing through the cell cycle inside the dark zone and the other one transitioning toward CCs in the light zone (Fig. 2d). The generated latent trajectories were mapped by the decoders to the original RNA and motif spaces to generate 100 cells, which were embedded together with the observed cells in UMAPs (Fig. 2d).
Using DynaVelo’s trajectory generation capability, we also examined the dynamics of Spi1 and Nfkb1, two critical TFs in GC B cells, along both CB → CC and CB cell cycle trajectories (Fig. 2e). Here, the velocity (RNA and motif) curves are the time derivatives of their corresponding state (RNA expression and motif accessibility) curves. A recent study showed that Spi1 is required for transitioning of CBs to CCs and acts as a pioneer factor to open chromatin, recruiting chromatin remodeling complexes and subsequently enabling binding of Nfkb factors30. Our results validated these findings by showing that Spi1 and Nfkb1 accessibilities monotonically increase from CB to CC but not in CB cell cycle progression, where Spi1 motif accessibility is stable while Nfkb1 monotonically decreases (Fig. 2e). Motif accessibility demonstrated stronger dynamics than RNA expression for these TFs, confirming the importance of this modality in the GC system. These results show that DynaVelo can determine the temporal sequence of TF activity and binding in silico, with the potential to discover novel TF codependencies.
DynaVelo identifies driver TFs for major GC trajectory types
Having confirmed that DynaVelo identifies major classes of GC B trajectories (CB cycling, CB → CC transition and exits toward prememory B cells and plasmablasts), we next analyzed dynamic patterns of TF expression and motif accessibility along latent time in each trajectory type. To this end, we first partitioned the cells on the basis of their trajectories, sorted them by latent time and plotted expression and motif accessibility temporal patterns of important GC TFs (Fig. 3a). We expected that driver TFs for each cell type would display monotonic temporal patterns within the cell type and the temporal order of TF activity might also have a role in cell fate decisions.
a, The temporal trends of motif accessibility and RNA expression for essential GC TFs across DynaVelo latent time in three main trajectory types: CB → CC → prememory, CB → CC → plasmablast and CB cell cycle. b, RNA expression and motif accessibility of several important TFs in each GC trajectory type shown on RNA and motif UMAPs, respectively. Both RNA expression and motif accessibility of TFs are required to understand their roles in each trajectory.
For the prememory B cell state, in addition to Spi1/Spib and Nfkb1/2, TFs including Runx1 and Irf8 displayed an increasing pattern, suggesting their possible roles in CC → prememory progression (Fig. 3a,b)30,31. UMAP embeddings in Fig. 3b show RNA expression and motif accessibility temporal patterns of Spi1 and Runx1, with accessibility of Spi1 gradually increasing from CB to CC to prememory B cells but not in plasmablasts. The RNA expression of Spi1 is also absent in plasmablasts. Notably, the Runx1 motif is not accessible in the GC but becomes accessible later when committing to the prememory B cell fate (Fig. 3b).
According to a similar DynaVelo analysis, important TFs driving the plasmablast fate include Prdm1, Irf4 and Xbp1 (ref. 31) (Fig. 3a). Although TF motif accessibility is a good proxy for single-cell TF activity, we found that a TF’s RNA expression was also required to accurately infer its role. For example, Prdm1 is accessible in both prememory B cells and plasmablasts but is only expressed in plasmablasts, showing that Prdm1 is only responsible for plasmablast differentiation (Fig. 3b). Similarly, Xbp1 exhibits both higher accessibility and expression in plasmablasts, consistent with its role as a plasma cell TF.
While key TFs for prememory B cells and plasmablasts have been studied, less is known about CB-specific drivers of cell cycle. DynaVelo temporal patterns identified Foxo1, Mef2b, Mef2c and Pou2f2, all well-known GC TFs31, as potentially having a role in CB cell cycle progression (Fig. 3a). These TFs are expressed in both CBs and CCs but their motif accessibility increases with CB cell cycle progression and decreases when transitioning toward CCs (Fig. 3a,b). These observations are consistent with known biology, as these TFs are important for maintaining the GC reaction but have a limited role in differentiation toward prememory B or plasmablast fates.
Heterozygous mutations in Arid1a and Ctcf disrupt the normal GC B cell trajectories
To determine whether DynaVelo can capture changes in cell fate dynamics upon perturbation of key GC regulators, we next turned to GC B cells with reduced abundance of Arid1a, a DNA-binding subunit of the BAF chromatin remodeling complex, or the insulator protein Ctcf, both a TF and mediator of 3D chromatin looping. We applied DynaVelo to two multiome samples, Arid1aHet (n = 6,716 cells) and CtcfHet (n = 4,688 cells), with heterozygous conditional deletion of Arid1a and Ctcf, respectively, in GC B cells30. Similar to WT, learned latent times in these samples identified a subset of CBs as root cells and prememory B cells and plasmablasts as terminal cells (Fig. 4a,b). The learned RNA velocities captured some of the main GC trajectories including CB cell cycle, CB → CC transition and differentiation toward prememory B cells and plasmablasts (Fig. 4a,b). However, one physiologically important trajectory was lost in both samples—the recirculation of CCs to CBs. This trajectory is necessary for BCR affinity maturation before differentiation to memory B cells or plasma cells, and its loss most likely leads to a premature exit from the GC at the prememory B cell state30. Indeed, the percentage of prememory B cells in the mutant samples (20.09% in Arid1aHet, 25.45% in CtcfHet and 44.96% in Arid1aHom) is increased compared to WT (5.42%) (Supplementary Fig. 6). This indicates that mutant GC B cells do not undergo multiple rounds of somatic hypermutation to improve their BCR affinity and consequently exhibit a faster exit as premature prememory B cells.
a, RNA and motif velocities of Arid1aHet GC B cells and DynaVelo’s latent times projected to the RNA and motif UMAPs. b, RNA and motif velocities of CtcfHet GC B cells and DynaVelo’s latent times projected to the RNA and motif UMAPs. In a,b, the critical recycling of CCs from the light back to the dark zone is lost, leading to either apoptosis or GC exit as premature prememory B cells. c, TF motif accessibility trends of WT, Arid1aHet and CtcfHet GC B cells in two main trajectory types: CB → CC → prememory and CB cell cycle. The dynamics of some TFs change more than others in the mutant GCs, highlighted by red and blue arrows for Arid1aHet and CtcfHet, respectively. d, Motif accessibility of TFs whose temporal patterns change more in Arid1aHet and CtcfHet GCs.
To identify key TFs disrupted in mutant GCs, we analyzed their altered dynamics in CB cell cycle progression and the CB → CC → prememory trajectory type compared to WT (Fig. 4c). While many TFs displayed subtle changes, several were especially notable for their motif accessibility changes in the Arid1aHet sample, including Foxp1, Bach2, Spi1 and Ctcf (Fig. 4c). The accessibility of Foxp1 and Bach2 decreases in Arid1aHet in both trajectory types compared to WT (Fig. 4c,d), while the accessibility of Spi1 decreases in Arid1aHet during CB cell cycle progression (Fig. 4c,d). Interestingly, Ctcf accessibility increases significantly in both CBs and CCs in Arid1aHet compared to WT (Fig. 4c,d), suggesting that Arid1a and Ctcf may function antagonistically within the GC. Overall, the CtcfHet displayed a stronger phenotype than Arid1aHet, in that TFs exhibited a stronger change in motif activity relative to WT. Additionally, we identified Sox5 as a potential TF whose accessibility decreases slightly in CtcfHet compared to WT, especially during CB cell cycle (Fig. 4c,d).
DynaVelo uncovers dynamic TF–gene relationships in GCs upon loss of Arid1a and Ctcf
So far, we have shown that DynaVelo can learn the major WT GC trajectories and their potential driver TFs and identify disrupted trajectories of mutant GC B cells. We next used DynaVelo to derive dynamic GRNs, comprising dynamic TF–gene regulatory relationships (‘edges’), through examination of Jacobian matrices at cells of each trajectory type or through in silico perturbation of genes, assessing how RNA expression and TF motif accessibility regulate the RNA and motif velocities. We first examined known gene regulatory relationships from the literature31 that were recovered on the basis of the Jacobian matrix ({J}_{{xx}}), which measures how the change in RNA velocity of one gene depends on the change in expression of another gene, leaving TF motif accessibilities unchanged (Fig. 5a). Positive and negative Jacobian values correspond to upregulation and downregulation, respectively. One limitation of this assessment is that known GRN ‘edges’ are reported as static rather than cell-state-dependent and dynamic31. DynaVelo extracted dynamic gene regulatory relationships at the single-cell level and identified temporal patterns along latent time in CB cell cycle and CB → CC → prememory trajectories in WT, Arid1aHet and CtcfHet samples (Fig. 5a,b). One key regulatory edge in the GC is Spi1 → Bcl-6 upregulation31, which matches the DynaVelo Jacobian in the WT CB → CC → prememory trajectory type. Notably, this regulatory relationship is lost in both Arid1aHet and CtcfHet along this trajectory (Fig. 5a).
a, Dynamic GRN edges for WT, Arid1aHet and CtcfHet GCs in two trajectory types (CB → CC → prememory and CB cell cycle) based on analysis of the Jacobian ({J}_{xx}). The majority of the WT dynamic regulatory edges match the known static GC GRN edges but provide more information about the cell state in which they hold. The WT edge, Spi1 → Bcl-6 upregulation, is missing in both Arid1aHet and CtcfHet GCs. b, Dynamic GRN edges for WT, Arid1aHet and CtcfHet GCs in two trajectory types (CB → CC → prememory and CB cell cycle) by in silico perturbations. Two important edges are lost in mutant GCs. Foxo1 → Cxcr4 upregulation is lost in CCs in Arid1aHet, which is required for the migration of CCs from the light to the dark zone, while Bcl-6 → Bcl-2 downregulation is gained in prememory B cells in Arid1aHet. c, The regulatory edges Foxo1 → Cxcr4 and Bcl-6 → Bcl-2 shown in the RNA UMAPs of WT, Arid1aHet and CtcfHet GCs. d, A minimal static GRN of the WT GC is summarized from the dynamic GRN analysis. Red arrows represent upregulation; blue arrows represent downregulation. The Bcl-6 → Bcl-2 repressive edge is lost in the WT GC to enable differentiation to memory B cells (blue cross). The edges that are lost or gained in mutant GC B cells are shown in purple and green lines and crosses.
As a second approach to identify TF regulatory effects, we performed in silico perturbation by setting both the TF RNA expression and motif accessibility to their minimal values and measuring the change in RNA velocities for all other genes. Note that in silico perturbation analysis captures slightly different regulation patterns from ({J}_{{xx}}) Jacobian analysis, as both TF RNA expression and motif accessibility are perturbed across cells. Figure 5b shows several important and known gene regulatory relationships in GCs31, along with the cell states and trajectory types in which they hold. Similar to the Jacobian analysis, in silico perturbation identified regulatory edges that are lost in Arid1aHet-mutant and CtcfHet-mutant GCs. One critical transcriptional GRN edge in WT GCs is Foxo1 → Cxcr4 upregulation, needed for maintenance of the GC and migration of CCs to the dark zone31. Notably, Foxo1 → Cxcr4 upregulation is lost in CCs in Arid1aHet but instead is observed in CB cell cycle, as compared to WT (Fig. 5b,c). This lost regulatory edge may be a reason why CCs in Arid1aHet lose their capability to migrate back to the dark zone and instead exit as premature prememory B cells, as described in our earlier analysis and experimentally validated30.
The inferred sign of Bcl-6 → Bcl-2 regulatory influence in the transition of CCs to prememory B cells is different in both mutants compared to WT (Fig. 5b,c). Bcl-6 is a master transcriptional repressor that downregulates expression of the antiapoptotic protein Bcl-2 (ref. 31). This repression is relieved in CCs that receive T cell help, allowing increased expression of Bcl-2 to promote B cell survival. In WT, the sign of Bcl-6 → Bcl-2 becomes positive, opposite to what we observe in both Arid1aHet and CtcfHet, suggesting that Arid1a and Ctcf regulate GC dynamics through control of the Bcl-6–Bcl-2 axis.
To summarize our findings, we assembled a minimal GRN governing the transcriptional regulation of critical GC TFs and highlighted regulatory relationships that are compromised in Arid1aHet and CtcfHet GCs (Fig. 5d).
In silico perturbations of Arid1a-mutant GCs identify candidate rescue targets
One powerful application of DynaVelo is to perform in silico loss-of-function perturbations and determine their effects on GC dynamics and phenotype. We assessed how accurately we could predict the impact of Arid1a genetic loss of function on RNA and motif velocities, as well as determine loss-of-function perturbation targets that could rescue the Arid1a loss-of-function phenotype. Notably, the change in the RNA and motif velocities (delta velocities) after in silico perturbation of Arid1a in the WT DynaVelo model suggests a decreased connectivity of CB and CCs, potentially reflecting the loss of the CC → CB recycling trajectory. Instead the trajectory between CCs and prememory B cells became highly prominent (Fig. 6a), suggesting a preferential prememory B cell differentiation and subsequent premature GC exit.
a, RNA velocities of WT, Arid1aHet and delta RNA velocities after in silico perturbation of Arid1a in the WT cells. Delta velocities show the difference between the velocity vectors before and after perturbation. These delta velocities are projected to the RNA UMAP derived from all samples, where the arrows indicate the direction toward which the cells are inclined to move after the perturbation. These trajectories show that the cells tend to transition from CB to CC and CC → CB recirculation is lost. b, Evaluation of predicted Arid1a KO in WT cells. DV overlap and direction agreement show how much the significant predicted RNA and motif velocities overlap with the true ones and how much they agree in sign (direction). c, The complete list of TFs for which Arid1a delta motif velocities match the motif velocity and accessibility changes in all cell types when comparing WT and Arid1aHet. d, Predicted rescue rates of RNA and motif velocities after in silico perturbation of essential TFs in Arid1aHet. These rates are the mean of rescue rates across all cell types. The TFs at the top right are the best perturbation targets for rescuing the loss-of-function phenotype in Arid1aHet. e, RNA velocities of Arid1aHet;CtcfHet-DKO cells and delta RNA velocity after in silico perturbation of Ctcf in the Arid1aHet cells. f, True rescue rates of RNA and motif velocities, RNA expression and motif accessibility in the DKO GC in each cell type. g, Bcl6 and Stat3 heterozygous KO in Arid1aHet GC (leading to two DKO GCs: Arid1aHet;Bcl6Het and Arid1aHet;Stat3Het) rescues the DZ/LZ ratio phenotype of Arid1aHet GCs, which is higher than WT GCs. Overall, 54 samples were used. Number of samples in Bcl-6 experiment: WT, 6; Arid1aHet, 6; Bcl6Het, 6; Arid1aHet;Bcl6Het-DKO, 6. Number of samples in Stat3 experiment: WT, 10; Arid1aHet, 10; Stat3Het, 6; Arid1aHet;Stat3Het-DKO, 4. Each point represents one mouse and data are presented as the mean values ± s.d. Genotypes were compared using a two-sided Welch’s t-test.
For a quantitative assessment, after in silico perturbation of Arid1a, we first performed a statistical test to identify RNA and motif delta velocities in each cell type whose mean was significantly different from zero, either positive (velocity increase) or negative (velocity decrease) (Methods). Then, we tested whether the mean velocities were significantly different between Arid1aHet and WT GCs in each cell type (Methods). We defined two metrics to assess the match between predicted and true significant (({p}_{mathrm{adj}} <0.05)) differential velocities (DV): DV overlap, defined as the Jaccard index (size of intersection divided by the size of union) of predicted and true DV sets, and DV direction agreement, defined as the fraction of velocities with the same sign in the DV intersection set. Figure 6b gives the results of Arid1a in silico KO in WT GCs, showing higher prediction accuracy for motif velocities than RNA velocities across cell types. Figure 6c shows the TFs whose predicted and true motif velocity changes match each other in all the cell types. Of note are changes in important GC TFs such as the upregulation of Foxo1 (mostly known as a dark-zone TF) and downregulation of Spi1, Irf4 and Prdm1.
As a final application of DynaVelo, we performed in silico perturbation of important GC TFs and epigenetic regulators in Arid1aHet GC B cells and identified perturbation targets that could potentially rescue the genetic loss of function (Fig. 6d). A velocity was counted as rescued in a given cell type if the in silico perturbation significantly shifted it toward its WT value. For each candidate perturbation in the Arid1aHet model, we computed the number of rescued RNA and motif velocities across cell types and defined the rescue rate as the DV direction agreement. Figure 6d shows the mean rescue rates (RNA and motif) for all cell types, identifying Tox, Bcor, Ctcf, Stat3 and Bcl6 as potential rescue perturbations for the Arid1aHet phenotype (Fig. 6d). Notably, for the perturbation target Ctcf, recently published work32 has established a mechanistic basis for the predicted antagonism between Ctcf and Arid1a and demonstrated a phenotypic rescue of the Arid1aHet DZ increase in animals with Arid1aHet;CtcfHet double-heterozygous KO (DKO) using flow cytometry, confirming a key result from our DynaVelo rescue analysis.
To quantitatively validate this rescue, we performed single-cell multiome on GC cells from Arid1aHet;CtcfHet-DKO mice and control animals. In the Arid1aHet;CtcfHet-DKO GCs, the DynaVelo-inferred RNA velocity field showed that the rapid CC to prememory exit characteristic of Arid1aHet was partially reversed (Fig. 6e, left). This experimental outcome closely mirrored the delta RNA velocity predicted by in silico Ctcf perturbation of Arid1aHet cells (Fig. 6e, right), confirming that DynaVelo’s prediction recapitulated the dynamics of the genetic double mutant. We quantified rescue across cell types as the fraction of features that shifted from Arid1aHet toward WT values in the DKO (denoted as the true rescue rate) across four modalities: RNA expression, RNA velocity, motif accessibility and motif velocity (Fig. 6f and Methods). Rescue was observed across all cell types but varied in magnitude. In CCs, the cell type whose dynamics were most disrupted in Arid1aHet, the rescue rates were 38% (RNA expression), 64% (RNA velocity), 51% (motif accessibility) and 46% (motif velocity).
A defining feature of Arid1aHet GCs is an increased ratio of dark zone to light zone (DZ/LZ), consistent with the premature exit of light zone cells toward the prememory fate (Fig. 6g). While DynaVelo in silico perturbation analysis nominated Bcl6 and Stat3, alongside Ctcf, as candidate rescue targets for Arid1a loss, an independent in vivo test is required to establish functional genetic interactions. We hypothesized that, if DynaVelo predictions are correct, compound heterozygous loss of either Bcl6 or Stat3 with Arid1a should restore the DZ/LZ ratio toward WT. To test this, we generated two entirely new mouse models carrying double-heterozygous mutations, Arid1aHet;Bcl6Het and Arid1aHet;Stat3Het, and performed immunophenotyping of GC B cells from each strain. The DZ/LZ ratio was restored to WT levels in both compound heterozygous backgrounds (Fig. 6g), directly confirming the rescue predicted by DynaVelo (Supplementary Note).
Together with the independent rescue of the same phenotype by Ctcf heterozygous loss in Arid1aHet mice, these experiments establish that all three rescue targets nominated by DynaVelo (Ctcf, Bcl6 and Stat3) functionally rescue an in vivo phenotype of Arid1a loss. These results demonstrate that DynaVelo’s in silico perturbation framework can identify genetic interactions and nominate experimentally tractable rescue strategies for loss-of-function disease alleles.
Discussion
A fundamental goal in biology is to understand the dynamic GRNs that govern cell state transitions in evolving cellular systems and to predict how cell fates change when the network is perturbed. As most single-cell genomics technologies are destructive, cell dynamics and GRNs are typically inferred from data at a single time snapshot using pseudotime strategies, including those based on RNA velocity. However, velocity estimates are prone to errors and can generally be computed only for a small subset of genes. Another key challenge is to distinguish causal drivers from noncausal ones, as single-cell data provide a high-dimensional space of observables comprising many confounding factors. With the advent of single-cell multiome technology, we can now exploit single-cell TF motif activities derived from scATAC-seq, in addition to RNA expression and velocity, in dynamic models of evolving cellular systems. We developed DynaVelo as a robust generative deep learning model that uses single-cell RNA and motif accessibility data at a single time point to infer dynamic gene regulation patterns, and we applied this model to dissect the complex in vivo dynamics of WT and mutant GC B cells.
DynaVelo uses a neural ODE framework where RNA expression and TF motif accessibility inputs are each mapped to latent spaces in which the cell dynamics are learned. Partial RNA velocities from a tool such as scVelo are supplied as input to provide pseudotemporal guidance, and the model learns RNA velocities for all modeled genes and TF motif velocities. Through Jacobian analysis and in silico perturbation of TFs, DynaVelo can reveal dynamic GRNs, where single-cell TF–gene regulatory edges depend on cell state, moving beyond traditional GRN inference methods that infer a single static network per dataset or cell type. Moreover, DynaVelo provides a generative deep learning framework, enabling the generation of trajectories and analysis of in silico genetic perturbations.
We first used DynaVelo in WT GCs to uncover the major trajectories of the GC reaction. GC B cells experience rapid and extensive shifts in their chromatin and transcriptional programs as they cycle back and forth between rounds of proliferative bursting in the GC dark zone and selection in the light zone and exit the GC reaction as memory B cells or plasma cells33. These rapid cell fate shifts are dictated by constantly evolving chromatin patterns, presenting a challenge to existing trajectory and velocity inference algorithms. DynaVelo recovered expected dynamics of GCs and the driver TFs of CB cell cycle progression, CB → CC transition and prememory and plasmablast exits and identified Spi1 and Nfkb1 as two driving TFs for CBs to transition to CCs and then to prememory B cells, consistent with previous findings30. Next, we applied DynaVelo to mutant GCs and found that both Arid1aHet and CtcfHet lost migration from the light zone back to the dark zone, a critical step required for B cell affinity maturation before differentiation to prememory B cells or plasmablasts. As a result, CCs either undergo apoptosis or transition toward premature prememory B cells, consistent with the observation that the percentage of prememory B cells is higher in the mutant GCs and that Arid1a haploinsufficiency in the GC leads to a fast exit as premature prememory B cells30. Through the analysis of TF motif accessibility along inferred latent time, we identified TFs whose dynamics were disrupted in mutant GCs, determining that Bach2, Foxp1 and Spi1 become less accessible and Ctcf becomes highly accessible in Arid1aHet.
Beyond resolving latent time ordering of cells, DynaVelo learns a functional form of the vector field underlying cell dynamics, enabling the inference of dynamic GRNs through Jacobian and in silico perturbation analysis. DynaVelo recovered TF–gene regulatory relationships in the WT GC that matched known biology and identified which regulatory edges are lost or gained in mutant GCs. In particular, upregulation of Cxcr4 by Foxo1 is essential for migration of CCs in the light zone back to the dark zone and this regulatory edge is lost in Arid1aHet CCs, leading to fast and premature exit from the GC. Moreover, Bcl-6 repression of Bcl-2 is relieved in CCs receiving T cell help before exit toward prememory B cells31 in WT GCs but this negative Bcl-6 → Bcl-2 regulatory relationship is apparently retained in both mutants. Lastly, we used in silico perturbations in the mutant GCs to identify candidate perturbation targets potentially leading to rescue of the loss of function in Arid1aHet. We then experimentally validated that heterozygous loss of predicted targets Ctcf, Stat3 and Bcl6 indeed rescued the dark-zone phenotype in Arid1aHet. DynaVelo, thus, enables the identification of TF drivers of cell state transitions, inference of dynamic GRNs and prediction of loss-of-function and rescue genetic perturbations.
Methods
Mice
All animal work was carried out in accordance with the Guide for the Care and Use of Laboratory Animals and the policies of Weill Cornell Medicine, with oversight consistent with AAALAC International accreditation standards. Protocols were reviewed and approved by the Research Animal Resource Center and the Institutional Animal Care and Use Committees of Weill Cornell Medicine and Cornell University (protocol nos. 2011-0031 and 2025-0008). Animals were housed in a specific-pathogen-free vivarium on a 12-h light–dark cycle, with temperature held at 20–22 °C and relative humidity at 40–60%, and were provided food and water ad libitum. Unless stated otherwise, experiments used age-matched and sex-matched mice between 8 and 16 weeks of age, with both males and females represented in every experimental group; no sex-dependent effects were detected. Arid1a (027717), Cg1Cre (010611) and Stat3 (016923) lines were obtained from The Jackson Laboratory. The CtcfHet (Ctcf+/–;Cγ1Cre/+) line was generously provided by the laboratory of G. Lauvau (Albert Einstein College of Medicine). Bcl-6 mice (Bcl6+/–;Cγ1Cre/+) were generated and published by A. Melnick34. Before use, the Arid1a mouse line was backcrossed to C57BL/6J for a minimum of six generations as the Jackson Laboratory strain was in a mixed genetic background. The Cγ1-Cre allele was used in heterozygous form throughout.
Immunizations
Mice received an intraperitoneal injection of 500 μl of a 2% suspension of sheep red blood cells in sterile 1× Dulbecco’s PBS, freshly prepared from sheep blood in Alsever’s solution (Cocalico Biologicals, 20-1334A).
Flow cytometry
Isolated spleens were mechanically dissociated and the resulting cell suspensions passed through a 40-μm strainer. Mononuclear cells were then enriched by density gradient centrifugation over Ficoll-Paque Premium (Cytiva, 17-5442-02). Single-cell suspensions were resuspended in PBE (PBS supplemented with EDTA) containing rat anti-mouse CD16/CD32 (clone 2.4G2, BD) at 0.5 μg ml−1 and incubated on ice for at least 5 min to block Fc receptors. Samples were stained on ice for 20 min with fluorochrome-conjugated anti-mouse antibody cocktails. Antibody mixes were prepared in PBE, except when the cocktail contained two or more reagents conjugated to BV or BUV, in which case the staining buffer consisted of a 50:50 mixture of brilliant stain buffer (BD, 566349) and PBS containing 0.5% BSA (no EDTA) to avoid polymer dye interactions. Stained cells were washed twice and DAPI (Thermo Fisher, D1306) was added immediately before acquisition at a final concentration of 0.5 μg ml−1 for dead-cell exclusion. Samples were run on a BD Symphony A3 cytometer and data were processed in FlowJo. Antibodies and dilutions used were as follows: B220–APC (RA3-6B2, BD 553092, 1:100), B220–BV786 (RA3-6B2, BD 563894, 1:300), CD138–BUV737 (281-2, BD 564430, 1:500), CCR6–PerCP–Cy5.5 (29-2L17, BioLegend 129809, 1:300), CD11b–APC (M1/70, BioLegend 101212, 1:400), CD21–APC (7E9, BioLegend 123412, 1:500), CD21–PECy7 (7E9, BioLegend 123420, 1:500), CD4–BUV563 (GK1.5, BD 612923, 1:500), CD38–BUV395 (clone 90, BD 740245, 1:500), CD38–BUV563 (90/CD38, BD 741271, 1:400), CD80–BV421 (16-10A1, BioLegend 104725, 1:300), CD86–BV605 (GL-1, BioLegend 105037, 1:200), CD138–BUV737 (281-2, BD 624230, 1:500), CXCR4–biotin (BD 551968, 1:200), CXCR4–PE (2B11, Invitrogen 12-9991-82, 1:100), Ephrin B1–biotin (R&D Systems, BAF473, 1:100), FAS–BUV805 (Jo2, BD 741968, 1:100), FAS–PE–Cy7 (Jo2, BD 557653, 1:500), GL7–FITC (BD 553666, 1:500), GL7–PerCP–Cy5.5 (GL7, BioLegend 144610, 1:500), IgD–BV510 (11-26c.2a, BD 563110, 1:500), IgD–PerCP (11-26c.2a, BioLegend 405709, 1:500), IgG1–BV650 (A85-1, BD 740478, 1:500), IgM–BUV661 (II/41, BD 750660, 1:300), IgM–BUV805 (II/41, BD 749307, 1:500), κ–APC–Cy7 (1:200), λ–BV421 (1:400) and PDL2–BUV395 (TY25, BD 565102, 1:75). Secondary detection reagents were as follows: SA–APC–Cy7 (BioLegend 405208, 1:100) and SA–BV605 (BioLegend 405229, 1:200). NP-PE (N-5070-1) was used at a working concentration of 5 µg ml−1.
Single-cell multiome (RNA + ATAC) sequencing
Paired scRNA-seq and scATAC-seq libraries were generated with the 10x Genomics Chromium single-cell multiome ATAC + gene expression kit (10x Genomics, 1000230, 1000283, 1000494, 1000215 and 1000212), following the manufacturer’s instructions. In brief, GC B cells were purified from splenocytes by fluorescence-activated cell sorting using the antibody panels described above and washed in PBS containing 0.04% BSA; a total of 100,000 sorted cells were processed for nucleus isolation per the manufacturer’s protocol, using a hybrid of the standard and low-input procedures. Approximately 20,000 resulting nuclei were loaded onto a 10x Chromium X instrument and downstream library construction was performed as specified by the manufacturer. Final libraries were assessed for quality on an Agilent Bioanalyzer, quantified with a KAPA library quantification kit and sequenced on an Illumina NovaSeq.
Inclusion and ethics statement
All murine studies described here were performed under protocol nos. 2011-0031 and 2025-0008, approved by the Institutional Animal Care and Use Committee at Weill Cornell Medicine, and were carried out in accordance with institutional policies and the Guide for the Care and Use of Laboratory Animals. Each experiment included both male and female animals in age-matched and sex-matched cohorts; no sex-dependent effects were observed.
Biological materials
Commercially sourced reagents and previously reported materials are cited in the relevant sections of the text. Novel biological resources generated as part of this work, namely, the Arid1aflox/WT/Ctcfflox/WT;Cγ1-Cre double conditional KO mouse model and Arid1aflox/WT/Stat3flox/WT;Cγ1-Cre, can be obtained from the corresponding authors upon reasonable request, potentially subject to a material transfer agreement.
scATAC data processing
We used ArchR35 to analyze our scATAC-seq data, starting from fragment files of our GC samples. ArchR tiles the genome with the resolution of 500 bp and then calculates the highly variable tiles. We used 25,000 highly variable genes and latent semantic indexing dimensionality reduction to bring the dimension to 30. We then used these clusters to create pseudobulk replicates and run the peak-calling algorithm by MACS2 (ref. 36).
Motif matrices from chromVAR
After getting the peak matrix from ArchR, we ran chromVAR26 to get the TF motif accessibility z scores in each single cell. We then smoothed these z scores over the 15 nearest cells on the k-nearest neighbor graph derived from LSI dimensions. We used the Vierstra motif set37 in chromVAR. This yielded the motif matrix for ~2,100 motifs. To filter the motifs with nonexpressed or lowly expressed TFs, we only kept the motifs whose TFs were expressed in at least 500 cells. As the TF motifs could be degenerate and multiple motifs could correspond to one TF, we only kept one motif per TF by choosing the pair that had the maximum Pearson correlation between the motif accessibility scores and RNA expression values. In the end, we ended up with a motif matrix for ~200 motifs, which was used as input for the DynaVelo model.
scRNA data processing
We used Scanpy38 to read the CellRanger outputs for all samples. We filtered the cells by keeping the ones with minimum gene count of 200, minimum unique molecular identifier (UMI) count of 1,000, maximum UMI count of 10,000 and percentage of mitochondrial genes less than 25%. We next performed total count normalization to the median value of all cells and log1p-normalized (log(x + 1)) them. Lastly, we kept 2,000 highly variable genes to be used in the DynaVelo model.
Spliced and unspliced counts
We ran Velocyto15 to read the CellRanger outputs and get the spliced and unspliced counts through loom files for each sample. We used Scanpy to read these loom files and merged the spliced and unspliced matrices to the RNA anndata, ready to be used by scVelo for RNA velocity estimates.
Estimated RNA velocities by scVelo
We used scVelo to get the estimated RNA velocities. As cells in each sample might have different dynamics, we ran scVelo per sample. scVelo needs sufficient spliced and unspliced counts to have a reliable RNA velocity estimation. We followed the scVelo pipeline and kept the genes with more than ten shared spliced and unspliced counts. This filtering step usually removed ~20,000 genes. We normalized the count, spliced and unspliced matrices on the basis of the median values and log1p-normalized the count matrix. When running scVelo, we included all genes, not only the highly variable ones. The reason for this choice was to have access to maximum number of velocity genes (genes for which scVelo provides RNA velocity estimates). Depending on the sample, the number of velocity genes was in the range of ~400 to ~1,000 genes.
Genes used in DynaVelo models
We chose four subsets of genes to be used in the DynaVelo models: highly variable genes (2,000 genes), velocity genes (~400–1000 genes), expressed TFs (~200 genes) and genes of interest (~20 genes). Highly variable genes carry more information about the cell state transitions. Velocity genes provide additional directional information to guide the learned dynamics. We also included the TFs that were expressed in at least 500 cells to be included in the GRNs, as well as the genes of interest in GCs such as Arid1a, which is an important epigenetic regulator often mutated in the lymphomas. We used the union of these four subsets (~3,000–4,000 genes) as our modeled genes.
DynaVelo model details
ChromVAR assigns cell-wise TF motif accessibility scores through genome-wide scanning of accessible peaks for TF motifs. Given the crucial role of TFs in regulating target genes within dynamic systems, we hypothesize that integrating TF motif accessibility scores could enhance the modeling of cellular dynamics. If we denote log-normalized RNA expression values of (m) genes in a cell at the latent time (t) by (xleft(tright)in {R}^{m}) and chromVAR z scores of (k) TFs in the same cell by (yleft(tright)in {R}^{k}), we model the cellular dynamics by a joint ODE as follows:
$$begin{array}{c}displaystylefrac{{rm{d}}xleft(tright)}{{rm{d}}t}={f}_{x}left(xleft(tright),yleft(tright)right), displaystylefrac{{rm{d}}yleft(tright)}{{rm{d}}t}={f}_{y}left(xleft(tright),yleft(tright)right),end{array}$$
(1)
where ({f}_{x}) and ({f}_{y}) are functions that model the dynamics (velocities) of RNA expression and TF motif accessibility, respectively. Here, we have full information of (xleft(tright)) and (yleft(tright)) and partial information of RNA velocities ({rm{d}}x(t)/{rm{d}}t) (for velocity genes) and we want to learn the dynamics functions ({f}_{x}) and ({f}_{y}) and the latent times (t) of the cells. After learning the dynamics functions, we recover the RNA velocities of all genes (including nonvelocity genes) and define the TF motif velocities ({rm{d}}y(t)/{rm{d}}t) to reveal dynamics of the chromatin landscape. Instead of working with ODEs in the original high-dimensional spaces, we define an equivalent joint ODE in low-dimensional latent spaces ({z}_{x}in {R}^{{d}_{x}}) and ({z}_{y}in {R}^{{d}_{y}}) and then map them to the original spaces. We set the default values ({d}_{x}={d}_{y}=50) in the models. The rationale behind this is twofold. First, we would like to have velocity-aware latent spaces that will be informative of dynamical patterns of the data. Second, as the high-dimensional data are noisy and sparse, learning the dynamics functions in the original spaces is computationally costly and less robust. We define the joint ODE in the latent space as follows:
$$begin{array}{c}displaystylefrac{{rm{d}}{z}_{x}left(tright)}{{rm{d}}t}={f}_{{z}_{x}}left({z}_{x}left(tright),{z}_{y}left(tright)right), displaystylefrac{{rm{d}}{z}_{y}left(tright)}{{rm{d}}t}={f}_{{z}_{y}}left({z}_{x}left(tright),{z}_{y}left(tright)right),end{array}$$
(2)
where ({z}_{x}left(tright)) and ({z}_{y}left(tright)) are respectively the latent space of RNA expression and TF motif accessibility of a cell at the latent time (t) and ({f}_{{z}_{x}}) and ({f}_{{z}_{y}}) are the dynamics functions of ({z}_{x}) and ({z}_{y}), respectively. If we denote the concatenations of these vectors by (zleft(tright)=left[{z}_{x}(t),{z}_{y}(t)right]) and ({f}_{z}=left[{f}_{{z}_{x}},{f}_{{z}_{y}}right]), we can write the joint latent ODE as (frac{{dz}left(tright)}{{dt}}={f}_{z}left(zleft(tright)right)), which has the solution (zleft(tright)=zleft(0right)+{int }_{0}^{t}{f}_{z}left(zleft(tright)right){rm{d}}t) and is known as the initial value problem. Because we do not know the functional form of ({f}_{z}), we define it as a neural network, which led to a neural ODE introduced in recent seminal work39. Neural ODEs are very powerful and efficient models for learning dynamics. After learning the joint cellular dynamics in the latent space, we map them back to the original spaces of RNA expression and motif accessibility by conditional distributions:
$$begin{array}{c}xleft(tright)|{z}_{x}left(tright) sim Nleft({mu }_{x}left(tright)={g}_{x}left({z}_{x}left(tright)right),{sigma }_{x}^{2}Iright), yleft(tright)|{z}_{y}left(tright) sim Nleft({mu }_{y}left(tright)={g}_{y}left({z}_{y}left(tright)right),{sigma }_{y}^{2}Iright),end{array}$$
(3)
where ({g}_{x}) and ({g}_{y}) are two neural network decoders mapping the latent space to the original spaces and ({sigma }_{x}^{2}) and ({sigma }_{y}^{2}) are the known variances of the Gaussian distributions (we assume ({sigma }_{x}^{2}={sigma }_{y}^{2}=0.01)). We define the RNA and motif velocities as the time derivatives of ({mu }_{x}left(tright)) and ({mu }_{y}left(tright)):
$$begin{array}{c}{v}_{x}left(tright)=displaystylefrac{d{mu }_{x}left(tright)}{{dt}}={{J}}_{{g}_{x}}times displaystylefrac{d{z}_{x}left(tright)}{{dt}}={{J}}_{{g}_{x}}{{f}}_{{z}_{x}}, {v}_{y}left(tright)=displaystylefrac{d{mu }_{y}left(tright)}{{dt}}={{J}}_{{g}_{y}}times displaystylefrac{d{z}_{y}left(tright)}{{dt}}={{J}}_{{g}_{y}}{{f}}_{{z}_{y}},end{array}$$
(4)
where ({{J}}_{{g}_{x}}=leftin {R}^{mtimes {d}_{x}}) and ({{J}}_{{g}_{y}}=leftin {R}^{ktimes {rm{d}}y}) are the Jacobian matrices of the decoder functions ({g}_{x}) and ({g}_{y}), respectively. We assume that the initial points of cells are unknown to have a general framework for less studied and more complex biological systems such as GC B cells. Therefore, we assume that (zleft(0right)) is a random variable with the prior standard Gaussian distribution (zleft(0right) sim Nleft(0,Iright)). We also assume that the latent times of cells (t) follow a prior beta distribution (tsimbeta (mathrm{2,2})), as the times should be positive. This leads to the expectation (Eleft=0.5) and variance ({Var}left=0.05). As such, the problem is to learn the posterior probability distributions of (zleft(0right)) and (t) for each cell, the dynamics function ({f}_{z}) and the decoder functions ({g}_{x}) and ({g}_{y}). To this end, we use variational inference techniques by first assuming approximate posterior distributions for (zleft(0right)) and (t) and then using inference networks (or encoders) to learn the parameters of those approximate posteriors, as shown in Fig. 1a. We assume the same family of distributions for the approximate posteriors as the priors, that is, ({z}_{x}(0)|x(t) sim N({mu }_{{z}_{x}(0)},{sigma^2}_{{z}_{x}(0)}I)), ({z}_{y}(0)|y(t) sim N({mu }_{{z}_{y}(0)},{sigma^2}_{{z}_{y}(0)}I)) and (t|x(t),y(t)simbeta ({alpha }_{t},{beta }_{t})), where ({mu }_{{z}_{x}left(0right)}), ({sigma }_{{z}_{x}left(0right)}), ({mu }_{{z}_{y}left(0right)}), ({sigma }_{{z}_{y}left(0right)}), ({alpha }_{t}) and ({beta }_{t}) are encoder neural networks that take in the inputs (xleft(tright)) and (yleft(tright)) and output the corresponding parameters of the approximate posterior distributions.
In addition to reconstruction loss, we use a Kullback–Leibler (KL) divergence loss function to ensure that the learned initial (latent time t = 0) latent states and cell-specific latent times do not greatly diverge from a standard normal distribution and a beta distribution, respectively. There are two additional loss functions used in DynaVelo. The first one is a cosine similarity loss function for the RNA velocities, where any reliable RNA velocity method, such as scVelo, is used to provide velocity estimates for a subset of genes and supervise the corresponding learned RNA velocities. The second loss function is a velocity consistency loss to ensure that RNA and TF motif velocities lead to similar trajectories of cells in their respective latent spaces. The total loss function has the following form:
$$L={L}_{mathrm{reconstruction}}+{kappa }_{1}{L}_{mathrm{velocity}}+{kappa }_{2}{L}_{mathrm{velocity},mathrm{consistency}}+{kappa }_{3}{L}_{mathrm{regularization}},$$
(5)
where the first term is the reconstruction loss of RNA expression and TF motif accessibility, which maximizes the log likelihoods:
$${L}_{mathrm{reconstruction}}=-{E}_{{z}_{0} sim {q}_{phi }({z}_{0}|{x}_{t},{y}_{t}),t sim {q}_{phi }({t|}{x}_{t},{y}_{t})}.$$
(6)
The second term is RNA velocity loss that maximizes the cosine similarity of predicted velocities and scVelo velocities for the velocity genes:
$${L}_{mathrm{velocity}}=-cos left({v}_{x}left[{G}_{v}right],{v}_{x}^{mathrm{scvelo}}left[{G}_{v}right]right),$$
(7)
where ({G}_{v}) and ({v}_{x}^{mathrm{scvelo}}in {R}^{m}) denote the indices of velocity genes as reported by scVelo and the estimated RNA velocities by scVelo, respectively. The third term is velocity consistency loss, which ensures that cells in two RNA and motif spaces follow consistent trajectories:
$$begin{array}{c}{L}_{mathrm{velocity},mathrm{consistency}}={||}{K}_{{v}_{x}}-{K}_{{v}_{y}}|{|}^{2}, {K}_{{v}_{x}}left=cos left({v}_{x}^{left(iright)},{v}_{x}^{left(jright)}right), {K}_{{v}_{y}}=cos ({v}_{y}^{left(iright)},{v}_{y}^{(j)}),end{array}$$
(8)
where ({K}_{{v}_{x}}) and ({K}_{{v}_{y}}) store respectively the cosine similarities of RNA and motif velocity vectors of the cells i and j. The velocity consistency loss minimizes the distance between cosine similarity of cells in both RNA and motif spaces, thereby helping to retain consistent trajectories in both spaces. The fourth term is regularization loss that minimizes the KL divergence between the approximate posterior and prior distributions of (zleft(0right)) and (t):
$${L}_{mathrm{regularization}}={KL}({q}_{phi }({rm{z}}left(0right)|{x}_{t},{y}_{t}){||p}({rm{z}}left(0right)))+{KL}({q}_{phi }({t|}{x}_{t},{y}_{t}){||p}(t))$$
(9)
Hyperparameter selection
Three important hyperparameters of DynaVelo are ({kappa }_{1}), ({kappa }_{2}) and ({kappa }_{3}), which control the strength of different components of the loss function. Depending on how much we trust scVelo velocities, we adjust ({kappa }_{1}), the coefficient of the velocity loss term. We advise first making sure that scVelo velocities overall make sense in the biological system to be modeled, which requires domain expertise. We would not recommend supervising DynaVelo velocities too strongly on scVelo if they do not capture the expected trajectories; in the current study, we confirmed that scVelo velocities captured the main trajectories of GC B cells. We used different values of ({kappa }_{1}) (1,000 and 10,000) during training and found that 10,000 was a good default value for this hyperparameter in the GC system.
The second hyperparameter, ({kappa }_{2}), controls how consistent the cell trajectories should be in RNA expression and motif accessibility space. Under a small value for ({kappa }_{2}), the trajectories might deviate too much between the two spaces. Therefore, we set ({kappa }_{2}) to the same value as ({kappa }_{1}) (10,000 as the default value) to make the trajectories consistent.
Lastly, ({kappa }_{3}) controls how close the learned initial point ({rm{z}}left(0right)) and latent time (t) of cells should be to their prior distributions. We found that ({kappa }_{3}) is very important in smoothing the learned latent times. If ({kappa }_{3}) is high, it forces all cells to start from a similar initial point ({rm{z}}left(0right)) and (t) becomes the only source of variation across cells, leading to smooth latent times aligned with the trajectories. Although this is satisfactory for learning the latent times, it might sacrifice the accuracy of predicted delta velocities upon in silico perturbations of genes. However, if ({kappa }_{3}) is low, cells have more freedom to start from any initial point ({rm{z}}left(0right)), which might improve the accuracy of perturbation prediction while leading to less smooth latent times that do not align well with the trajectories. Therefore, we believe that there is no universal ({kappa }_{3}) that leads to the best performance in all downstream tasks; rather, this hyperparameter should be selected as a function of the task. We tested different values for ({kappa }_{3}) (1, 10, 100 and 1,000) and found 1,000 to be a good default for learning velocities and latent times and 1 to be a good value for perturbation prediction tasks.
Generating trajectories using DynaVelo
As DynaVelo is a generative model, it can be used to generate cell trajectories. Given any cell (left{xleft(tright),y(t)right}), DynaVelo first maps it to the initial latent state (zleft(0right)=left[{z}_{x}(0),{z}_{y}(0)right]) and finds its latent time (t), which is simply the mean of the posterior beta distribution: (t=frac{{alpha }_{t}}{{alpha }_{t}+{beta }_{t}}in (mathrm{0,1})). For generating trajectories, we start from (zleft(0right)) at time (t=0) and solve the ODEs up to the time (t=1) with the step size of ({t}_{rm{step}}=0.01):
$$zleft({t}_{k}right)=zleft(0right)+{int }_{0}^{{t}_{k}}{f}_{z}left(zleft(tright)right){rm{d}}t,qquad{t}_{k}=k{t}_{mathrm{step}},qquad{k=1},ldots ,100,$$
(10)
which generates a latent trajectory ({T}_{z}=left) with the length of 101 (including the initial point). We then map this latent trajectory to the observational spaces of RNA expression and motif accessibility using their corresponding decoders ({mu }_{x}left(tright)={g}_{x}left({z}_{x}left(tright)right)) and ({mu }_{y}left(tright)={g}_{y}left({z}_{y}left(tright)right)) to generate the trajectories ({T}_{x}=left[{mu }_{x}left(0right),{mu }_{x}left({t}_{1}right),,ldots ,{mu }_{x}left({t}_{100}right)right]) and ({T}_{y}=left[{mu }_{y}left(0right),{mu }_{y}left({t}_{1}right),ldots ,{mu }_{y}left({t}_{100}right)right]). To quantify the uncertainty of these trajectories, we sample (N=50) initial latent states from the learned posterior distributions ({z}_{x}left(0right)|xleft(tright) sim)(N({mu }_{{z}_{x}left(0right)},{sigma^2}_{{z}_{x}left(0right)}I)) and ({z}_{y}left(0right)|y(t) sim N({mu }_{{z}_{y}left(0right)},{sigma^2}_{{z}_{y}left(0right)}I)) and generate (N) latent trajectories ({T}_{z}^{(n)}=left[{z}^{(n)}left(0right),{z}^{(n)}left({t}_{1}right),ldots ,{z}^{(n)}left({t}_{100}right)right]) and (N) corresponding RNA expression and motif accessibility trajectories ({T}_{x}^{(n)}=left[{mu }_{x}^{(n)}left(0right),{mu }_{x}^{(n)}left({t}_{1}right),,ldots ,{mu }_{x}^{(n)}left({t}_{100}right)right]) and ({T}_{y}^{(n)}=left[{mu}_{y}^{(n)}left(0right)right)[{mu}_{y}^{(n)}left(0right)right)(left. {mu }_{y}^{(n)}left({t}_{1}right),ldots ,{mu }_{y}^{(n)}left({t}_{100}right)right]) for (n=1,ldots ,N.) We then calculate the mean and s.e.m. at each time point of the RNA trajectories as (the calculation is similar for motif accessibility trajectories):
$$begin{array}{c}{overline{T}}_{x}=left[{overline{mu }}_{x}left({t}_{k}right)right],mathrm{where},{overline{mu }}_{x}left({t}_{k}right)=frac{sum _{n=1}^{N}{mu }_{x}^{(n)}left({t}_{k}right)}{N},qquad{k=0},1,ldots ,100, {rm{s}}. {rm{e}}. {rm{m}}.left({T}_{x}right)=frac{sigma left({T}_{x}right)}{sqrt{N}},end{array}$$
(11)
where (sigma left({T}_{x}right)) is the s.d. at each time point:
$$begin{array}{l}sigma left({T}_{x}right)=left,mathrm{where},sigma left({mu }_{x}left({t}_{k}right)right)=sqrt{frac{1}{N}mathopsumlimits _{n=1}^{N}{left({mu }_{x}^{left(nright)}left({t}_{k}right)-{overline{mu }}_{x}left({t}_{k}right)right)}^{2}},qquad{k=0},ldots ,100end{array}$$
(12)
Lastly, we calculate the 95% confidence interval of the mean trajectory as
$$begin{array}{c}{mathrm{Upper},95 % ,mathrm{limit}=overline{T}}_{x}+1.96times {rm{s}}. {rm{e}}. {rm{m}}.left({T}_{x}right) {mathrm{Lower},95 % ,mathrm{limit}=overline{T}}_{x}-1.96times {rm{s}}. {rm{e}}. {rm{m}}.left({T}_{x}right)end{array}$$
(13)
Calculating the Jacobian matrices for dynamic GRNs
We can calculate Jacobian matrices, evaluated at individual cells, that represent the derivatives of the RNA and TF motif velocity outputs with respect to the RNA expression and TF motif accessibility inputs. This results in four Jacobian matrices in each cell showing how (1) RNA expression affects RNA velocity; (2) RNA expression affects motif velocity; (3) TF motif accessibility affects RNA velocity; and (4) TF motif accessibility affects motif velocity. Specifically, once we have learned the velocity functions ({f}_{x}) and ({f}_{y}) in the joint ODE model, we define four Jacobian matrices for the cell (c) as follows:
$$begin{array}{c}begin{array}{c}{J}_{{xx}}^{(c)}left={left}_{{x=x}^{(c)},y={y}^{(c)}},qquad{i=1},ldots ,m,qquad{j=1},ldots ,m, {J}_{{yx}}^{(c)}left={left}_{{x=x}^{(c)},y={y}^{(c)}},qquad{i=1},ldots ,k,qquad{j=1},ldots ,m,end{array} begin{array}{c}{J}_{{xy}}^{(c)}left={left}_{{x=x}^{(c)},y={y}^{(c)}},qquad{i=1},ldots ,m,qquad{j=1},ldots ,k, {J}_{{yy}}^{(c)}left={left}_{{x=x}^{(c)},y={y}^{(c)}},qquad{i=1},ldots ,k,qquad{j=1},ldots ,k,end{array}end{array}$$
(14)
for (c=1,ldots ,{n}_{mathrm{cells}}), where ({J}_{{xx}}) and ({J}_{{yx}}) show how the RNA expression of a gene affects all other RNA and motif velocities, respectively, and ({J}_{{xy}}) and ({J}_{{yy}}) show how the motif accessibility of a TF affects all other RNA and motif velocities. This yields dynamic gene regulation information as the Jacobians are cell specific. Because GRNs reported in the literature are static rather than dynamic, we calculate summary statistics of these Jacobian matrices in specific cell states to obtain static cell-state-specific GRNs. For example, we obtain the transcriptional GRNs in cell state (S) as
$${J}_{{xx}}^{S}=sum _{cin S}{J}_{{xx}}^{(c)}$$
(15)
To calculate the Jacobians, we consider the definition of partial derivatives:
$$begin{array}{c}begin{array}{c}{J}_{{xx}}^{(c)}left[:,jright]=mathop{mathrm{lim}}limits_{epsilon to 0}frac{{f}_{x}left({x}^{(c)}+epsilon {e}_{j},{y}^{(c)}right)-{f}_{x}left({x}^{left(cright)},{y}^{(c)}right)}{epsilon },qquad{j=1},ldots ,m, {J}_{{yx}}^{(c)}left[:,jright]=mathop{mathrm{lim}}limits_{epsilon to 0}frac{{f}_{y}left({x}^{(c)}+epsilon {e}_{j},{y}^{(c)}right)-{f}_{y}left({x}^{left(cright)},{y}^{(c)}right)}{epsilon },qquad{j=1},ldots ,m,end{array} begin{array}{c}{J}_{{xy}}^{(c)}left[:,jright]=mathop{mathrm{lim}}limits_{epsilon to 0}frac{{f}_{x}left({x}^{(c)},{y}^{(c)}+epsilon {e}_{j}right)-{f}_{x}left({x}^{left(cright)},{y}^{(c)}right)}{epsilon },qquad{j=1},ldots ,k, {J}_{{yy}}^{(c)}left[:,jright]=mathop{mathrm{lim}}limits_{epsilon to 0}frac{{f}_{y}left({x}^{(c)},{y}^{(c)}+epsilon {e}_{j}right)-{f}_{y}left({x}^{left(cright)},{y}^{(c)}right)}{epsilon },qquad{j=1},ldots ,k,end{array}end{array}$$
(16)
for (c=1,ldots ,{n}_{mathrm{cells}}), where ({e}_{j}) is the unit vector of (j) th element in either ({x}^{(c)}) or ({y}^{(c)}); we use a small value (epsilon={10}^{-4}) to estimate the values.
In silico perturbations
To perform in silico perturbation of the ith non-TF gene, we set its RNA expression to the minimum value over all cells and calculate the change in both RNA and motif velocity vectors (delta velocity) as follows:
$$begin{array}{c}{x}^{left(cright),(ito {rm{pert}})}=left[{x}^{left(cright)}left[1right],ldots ,mathop{min }limits_{c}{x}^{left(cright)}left,ldots ,{x}^{left(cright)}leftright], Delta {v}_{x}^{left(cright),(ito {rm{pert}})}={f}_{x}left({x}^{left(cright),(ito {rm{pert}})},{y}^{(c)}right)-{f}_{x}left({x}^{left(cright)},{y}^{left(cright)}right), Delta {v}_{y}^{left(cright),(ito {rm{pert}})}={f}_{y}left({x}^{left(cright),(ito {rm{pert}})},{y}^{(c)}right)-{f}_{y}left({x}^{left(cright)},{y}^{left(cright)}right),end{array}$$
(17)
for (c=1,ldots ,{n}_{mathrm{cells}}). To perform in silico perturbation of ith TF, we set both its RNA expression and motif accessibility to the minimum value over all cells and calculate the delta velocities as (if the ith element in both (x) and (y) vectors corresponds to the same TF):
$$begin{array}{c}begin{array}{c}{x}^{left(cright),(ito {rm{pert}})}=left[{x}^{left(cright)}left[1right],ldots ,mathop{min }limits_{c}{x}^{left(cright)}left,ldots ,{x}^{left(cright)}leftright], {y}^{left(cright),(ito {rm{pert}})}=left[{y}^{left(cright)}left[1right],ldots ,mathop{min }limits_{c}{y}^{left(cright)}left,ldots ,{y}^{left(cright)}leftright],end{array} begin{array}{c}Delta {v}_{x}^{left(cright),(ito {rm{pert}})}={f}_{x}left({x}^{left(cright),(ito {rm{pert}})},{y}^{left(cright),(ito {rm{pert}})}right)-{f}_{x}left({x}^{left(cright)},{y}^{left(cright)}right), Delta {v}_{y}^{left(cright),(ito {rm{pert}})}={f}_{y}left({x}^{left(cright),(ito {rm{pert}})},{y}^{left(cright),(ito {rm{pert}})}right)-{f}_{y}left({x}^{left(cright)},{y}^{left(cright)}right),end{array}end{array}$$
(18)
for (c=1,ldots ,{n}_{mathrm{cells}}). These delta velocities can be projected to lower-dimensional spaces (such as UMAPs) to visualize cell behavior and trajectory changes after any perturbation. We note that these delta velocities are similar to the Jacobian matrices in that they measure the change in velocity functions after any perturbations. The main difference is that the perturbations are infinitesimal in Jacobians but finite in delta velocities. As a result, we can build two delta velocity matrices:
$$begin{array}{c}{Delta }{v}_{x}^{left(cright)}=-left[{Delta }{v}_{x}^{left(cright),left(1to mathrm{pert}right)},ldots ,{Delta }{v}_{x}^{left(cright),left(mto mathrm{pert}right)}right]in {R}^{mtimes m}, {Delta }{v}_{y}^{left(cright)}=-left[{Delta }{v}_{y}^{left(cright),left(1to mathrm{pert}right)},ldots ,{Delta }{v}_{y}^{left(cright),left(mto mathrm{pert}right)}right]in {R}^{ktimes m},end{array}$$
(19)
to be used for the dynamic GRN inference. The reason for the negative sign is because the perturbations were KO (setting to minimum values) instead of overexpression (setting to maximum values). These delta velocities can be used to understand both transcriptional ((Delta {v}_{x}^{left(cright)})) and epigenetic regulation ((Delta {v}_{y}^{left(cright)})) in each cell (c=1,ldots ,{n}_{mathrm{cells}}), or any cell state (S_{j}):
$$begin{array}{c}Delta {v}_{x}^{(S_{j})}=sum _{cin S_{j}}Delta {v}_{x}^{left(cright)}, Delta {v}_{y}^{(S_{j})}=sum _{cin S_{j}}Delta {v}_{y}^{left(cright)},end{array}$$
(20)
for (j=1,ldots ,{n}_{rm{states}}), where ({n}_{mathrm{states}}) is the number of cell states or types. We note that the dynamic GRNs derived by Jacobians and delta velocities are not identical because of the nature of perturbations performed in each approach.
Calculating predicted rescue rates of RNA and motif velocities
To calculate the predicted rescue rates of RNA and motif velocities after any TF perturbation in Fig. 6d, we first let ({v}_{x}^{({rm{mut}})}) and ({v}_{y}^{(mathrm{mut})}) be the RNA and motif velocity vectors of the mutant sample (Arid1aHet) and ({v}_{x}^{(mathrm{wt})}) and ({v}_{y}^{(mathrm{wt})}) be the RNA and motif velocity vectors of the WT sample. We perform in silico perturbation of the TF i in the mutant DynaVelo model and perform a one-sample t-test on the delta velocity vectors in each cell type (S_{j}), (j=1,ldots ,{n}_{rm{states}}), (Tleft({Delta }{v}_{x}^{left(mathrm{mut}right),left(S_{j}right),left(ito mathrm{pert}right)}right)) and (Tleft({Delta }{v}_{y}^{left(mathrm{mut}right),left(S_{j}right),left(ito mathrm{pert}right)}right)), to test whether the mean of delta velocities are significantly positive or negative (that is, whether this perturbation increases or decreases the velocities in cell type (S_{j})). To see whether the direction of change in velocities after in silico perturbation matches the direction toward the WT sample, we perform another two-sample t-test between the velocity vectors of WT and mutant samples in cell type (S_{j}), that is, (Tleft({v}_{x}^{left({rm{wt}}right),left(S_{j}right)},{v}_{x}^{left({rm{mut}}right),left(S_{j}right)}right)) and (Tleft({v}_{y}^{left(mathrm{wt}right),left(S_{j}right)},{v}_{y}^{left(mathrm{mut}right),left(S_{j}right)}right)). Then, we get the significant DV genes and motifs from (Tleft({v}_{x}^{left({rm{wt}}right),left(S_{j}right)},{v}_{x}^{left({rm{mut}}right),left(S_{j}right)}right)) and (Tleft({v}_{y}^{left(mathrm{wt}right),left(S_{j}right)},{v}_{y}^{left(mathrm{mut}right),left(S_{j}right)}right)) by keeping those with ({p}_{mathrm{adj}} <0.05). Now, we calculate the rescue rates of RNA and motif velocities as the respective fraction of significant DV genes and motifs where the sign between (Tleft(Delta {v}_{x}^{left({rm{mut}}right),left(S_{j}right),left(ito {rm{pert}}right)}right)) and (Tleft({v}_{x}^{left(mathrm{wt}right),left(S_{j}right)},{v}_{x}^{left(mathrm{mut}right),left(S_{j}right)}right)) and between (Tleft({Delta }{v}_{y}^{left(mathrm{mut}right),left(S_{j}right),left(ito mathrm{pert}right)}right)) and (Tleft({v}_{y}^{left(mathrm{wt}right),left(S_{j}right)},{v}_{y}^{left(mathrm{mut}right),left(S_{j}right)}right)) matches in each cell types (S_{j}). Then, we calculate the mean rescue rates over all the cell types, as plotted in Fig. 6d. This shows how much a TF perturbation in a mutant sample can shift RNA and motif velocities toward those of the WT sample.
Calculating DV overlap and direction agreement
We can perform an in silico perturbation of the gene i in the WT sample, predict the delta velocities and check how much these shift toward the velocities of the mutant sample, as shown in Fig. 6b for the perturbation of Arid1a. To this end, we perform one-sample t-tests (Tleft(Delta {v}_{x}^{left({rm{wt}}right),left(S_{j}right),left(ito {rm{pert}}right)}right)) and (Tleft({Delta }{v}_{y}^{left(mathrm{wt}right),left(S_{j}right),left(ito mathrm{pert}right)}right)) and two-sample t-tests (Tleft({v}_{x}^{left({rm{mut}}right),left(S_{j}right)},{v}_{x}^{left({rm{wt}}right),left(S_{j}right)}right)) and (Tleft({v}_{y}^{left(mathrm{mut}right),left(S_{j}right)},{v}_{y}^{left(mathrm{wt}right),left(S_{j}right)}right)) in cell type (S_{j}), (j=1,ldots ,{n}_{mathrm{states}}). Now, we define the DV overlap as the number of overlapped significant DV genes and motifs (({p}_{mathrm{adj}} <0.05)) between the one-sample and two-sample tests over the number of genes and motifs. Then, we define the DV direction agreement as the number of overlapped significant DV genes for which (Tleft(Delta {v}_{x}^{left({rm{wt}}right),left(S_{j}right),left(ito {rm{pert}}right)}right)) and (Tleft({v}_{x}^{left(mathrm{mut}right),left(S_{j}right)},{v}_{x}^{left(mathrm{wt}right),left(S_{j}right)}right)) have matched signs in each (S_{j}) and overlapped significant DV motifs for which (Tleft(Delta {v}_{y}^{left({rm{wt}}right),left(S_{j}right),left(ito {rm{pert}}right)}right)) and (Tleft({v}_{y}^{left(mathrm{mut}right),left(S_{j}right)},{v}_{y}^{left(mathrm{wt}right),left(S_{j}right)}right)) have matched signs in each (S_{j}) over the number of DV overlap genes and motifs, respectively. Higher values of DV overlap and direction agreement mean more correct velocity predictions after in silico perturbations. True rescue rates reported in Fig. 6f for Arid1aHet;CtcfHet-DKO GC are also based on direction agreements in all four data modalities: RNA expression, RNA velocity, motif accessibility and motif velocity.
Model architecture details
DynaVelo is a multimodal variational neural ODE framework that jointly models RNA expression and TF motif accessibility dynamics from paired single-cell multiome data. The model takes as input RNA expression vectors (xin {{mathbb{R}}}^{m}) and TF motif accessibility vectors (xin {{mathbb{R}}}^{k}), where (m) and (k) denote the numbers of genes and TF motifs, respectively. RNA expression and TF motif accessibility are first encoded independently into modality-specific latent spaces using shallow neural encoders consisting of a single fully connected layer followed by GELU activation. Specifically, the RNA encoder maps (x) through a linear layer of size (left) followed by GELU, while the motif accessibility encoder maps (y) through a linear layer of size (left) followed by GELU. Default latent dimensions were (d_x=d_y=50). The posterior distributions of latent initial states ({z}_{x}left(0right)) and ({z}_{y}left(0right)) are parameterized as diagonal Gaussian distributions. Separate linear layers estimate the latent means and variances for each modality, with Softplus activation applied to the variance branches to ensure positivity. Latent initial states are sampled using the reparameterization trick and concatenated into a shared latent state (zleft(0right)=left[{z}_{x}left(0right),{z}_{y}left(0right)right]in {{{{mathbb{R}}}}}^{{d}_{x}+{d}_{y}}.)
Latent cellular time is inferred probabilistically using a variational beta distribution conditioned jointly on RNA expression and motif accessibility. Separate time encoders, architecturally identical to the expression encoders, map RNA expression and motif accessibility into latent embeddings of dimensions ({d}_{x}) and ({d}_{y}), respectively. These embeddings are concatenated and passed through two independent linear layers with Softplus activations to infer the beta distribution parameters ({alpha }_{t}) and ({beta }_{t}). The latent time posterior is, therefore, defined as (qleft(t,|,x,yright)=beta left({alpha }_{t},{beta }_{t}right)), with a symmetric beta prior (beta left(mathrm{2,2}right)). During training, latent times are sampled using differentiable reparameterized sampling. Sampled times are sorted before ODE integration to ensure monotonically increasing integration intervals and numerical stability.
Latent dynamics are governed by a shared neural ODE defined on the concatenated latent state. The default ODE function consists of two parallel subnetworks independently predicting RNA latent velocities and TF motif latent velocities. Each subnetwork contains a fully connected layer mapping ({d}_{x}+{d}_{y}to h), followed by GELU activation and a second linear layer mapping (hto {d}_{x}) or (hto {d}_{y}), respectively, where the hidden dimension was typically (h=200). The outputs are concatenated to produce the latent velocity field ({f}_{z}left(zleft(tright)right)). All linear layers in the ODE function were initialized with Gaussian weights ((mu=0,sigma=0.01)) and zero biases. Latent trajectories were solved using the Dormand–Prince adaptive ODE solver (dopri5) from torchdiffeq with adjoint sensitivity optimization enabled during training.
Decoded RNA expression and motif accessibility are reconstructed independently from the latent states. The RNA decoder consists of a single linear layer mapping ({d}_{x}to m) followed by ReLU activation to enforce nonnegative reconstructed expression values. The TF accessibility decoder consists of a single linear layer mapping ({d}_{y}to k) without output activation, allowing positive and negative chromVAR z scores. Observed space velocities were computed by projecting latent velocities through decoder Jacobians using the chain rule. Specifically, decoder Jacobians were computed with automatic differentiation using vectorized Jacobian evaluation (vmap(jacrev)), and observed velocities were obtained as ({v}_{x}left(tright)={{J}}_{{g}_{x}}{{f}}_{{z}_{x}}) and ({v}_{y}left(tright)={{J}}_{{g}_{y}}{{f}}_{{z}_{y}}), where ({{J}}_{{g}_{x}}) and ({{J}}_{{g}_{y}}) denote the Jacobians of the RNA and TF motif decoders with respect to their latent representations.
The training objective combined reconstruction likelihoods, variational regularization, RNA velocity supervision and cross-modal velocity consistency. RNA and TF motif accessibility reconstructions were modeled with Gaussian likelihoods having fixed variance (({sigma }_{x}^{2}={sigma }_{y}^{2}=0.01)). KL divergence penalties regularized both the latent initial states toward a standard normal prior and latent times toward the beta prior. RNA velocity supervision was implemented using negative cosine similarity between predicted RNA velocities and estimated scVelo RNA velocities, masked to genes with valid velocity estimates. To encourage coordinated multimodal dynamics, a velocity consistency loss minimized the squared difference between pairwise cosine similarity matrices computed from predicted RNA and motif velocities across cells. The total loss was defined as the sum of reconstruction losses, KL penalties, velocity supervision and consistency regularization, with default weighting coefficients ({k}_{t}={k}_{z0}=1000) and ({k}_{mathrm{velocity}}={k}_{mathrm{consistency}}=10,000). During optimization, cells were additionally weighted inversely proportional to cell type abundance to mitigate imbalance across populations.
Model training and inference details
Model optimization was performed using the Adam optimizer with an initial learning rate of 1 × 10−3. For each genotype or sample, cells were randomly partitioned into training and validation sets using a 90:10 split. The validation set was used for monitoring generalization performance during training. Learning rate scheduling was performed adaptively using a plateau-based scheduler with a multiplicative decay factor of 0.5 and patience of five epochs. Early stopping was applied using the validation loss with a patience of ten epochs and the model parameters corresponding to the lowest validation loss were retained. Training was performed on a single NVIDIA RTX A6000 GPU (49 GB of VRAM). We used batch sizes of 128 and 256 depending on dataset size and memory constraints. Training on datasets containing 5,000–20,000 cells typically required at most 3 h when using a batch size of 256. We used the maximum epoch size of 200 but early stopping was usually triggered before then.
For probabilistic inference of latent trajectories and velocities after training, the model supports stochastic posterior sampling of both latent initial states and latent times. During prediction, RNA and motif velocities were estimated by repeatedly sampling from the posterior distributions of latent initial state and latent times. Specifically, the model was evaluated ({n}_{mathrm{samples}}) times (default ({n}_{mathrm{samples}}=100)), generating independent samples of predicted latent states, latent velocities, RNA and motif velocities. Final predicted quantities were obtained by averaging across posterior samples, while posterior uncertainty estimates were computed using the corresponding sample s.d. This Monte Carlo procedure reduces sensitivity to stochastic latent sampling and provides more robust velocity estimates.
For Jacobian estimation, stochastic sampling was disabled to obtain deterministic dynamical functions. In this setting, the posterior means of latent initial states and latent times were used directly instead of sampled values. Consequently, Jacobians were computed deterministically from the learned mean latent trajectories, enabling stable estimation of local regulatory interactions and perturbation responses.
In silico perturbation could also be performed in both deterministic and stochastic modes. In deterministic mode, perturbation effects were estimated using posterior mean latent variables and latent times, resulting in faster inference. In stochastic mode, perturbation effects were averaged across repeated posterior samples, analogous to velocity estimation. The stochastic perturbation procedure generally provides more robust estimates, although computational cost scales approximately linearly with the number of posterior samples. Users can, therefore, trade computational efficiency for robustness depending on the application and dataset size. We use the deterministic mode as a default.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
