Quantifying Landscape‐Flux via Single‐Cell Transcriptomics Uncovers the Underlying Mechanism of Cell Cycle

Abstract Recent developments in single‐cell sequencing technology enable the acquisition of entire transcriptome data. Understanding the underlying mechanism and identifying the driving force of transcriptional regulation governing cell function directly from these data remains challenging. This study reconstructs a continuous vector field of the cell cycle based on discrete single‐cell RNA velocity to quantify the single‐cell global nonequilibrium dynamic landscape‐flux. It reveals that large fluctuations disrupt the global landscape and genetic perturbations alter landscape‐flux, thus identifying key genes in maintaining cell cycle dynamics and predicting associated functional effects. Additionally, it quantifies the fundamental energy cost of the cell cycle initiation and unveils that sustaining the cell cycle requires curl flux and dissipation to maintain the oscillatory phase coherence. This study enables the inference of the cell cycle gene regulatory networks directly from the single‐cell transcriptomic data, including the feedback mechanisms and interaction intensity. This provides a golden opportunity to experimentally verify the landscape‐flux theory and also obtain its associated quantifications. It also offers a unique framework for combining the landscape‐flux theory and single‐cell high‐through sequencing experiments for understanding the underlying mechanisms of the cell cycle and can be extended to other nonequilibrium biological processes, such as differentiation development and disease pathogenesis.


Supplementary Tables
Table S3.The hill function and parameters of gene interaction.

ODEs Interaction term Parameters
The interaction between CCNE2 and CENPE The interaction between CCNE2 and KIF23

Figure S2 .
Figure S2.The differential geometric analysis of cell cycle dynamics of U2OS cells.(A) Ddhodge potential of the reconstructed vector field among all cell cycle phases.(B) Same as in (A) but for the speed.(C) Same as in (A) but for the divergence.(D) Same as in (A) but for the curl.(E) Same as in (A) but for the acceleration.(F) Same as in (A) but for the curvature.

Figure S3 .
Figure S3.Constructing cell cycle landscape-flux by scEU-seq data of RPE1 cells.(A) Cell cycle phase clusters by UMAP.(B) RNA velocity of cell cycle dynamics.(C) Reconstructed vector field of cell cycle dynamics.The red digit 0 to represent the limit cycle attractor, reflecting cell cycle oscillations, and the black digit 1 to represent the absorption fixed-point attractor, reflecting the stationary phase of the cell at G0. (D) Potential landscape of cell cycle dynamics in UMAP.(E) Curl flux (black arrow) of cell cycle dynamics landscape in UMAP.(F) The gradient force (black arrow) of cell cycle dynamics landscape in UMAP.

Figure S4 .
Figure S4.The differential geometric analysis of cell cycle dynamics of RPE1 cell.(A) Divergence of the reconstructed vector field among all cell cycle phases.(B) Same as in (A) but for the curl.(C) Same as in (A) but for the acceleration.(D) Same as in (A) but for the curvature.

Figure S5 .
Figure S5.Constructing cell cycle landscape-flux by scRNA-seq data of human fibroblasts.(A) Cell cycle phase clusters by UMAP.(B) RNA velocity of cell cycle dynamics.(C) Reconstructed vector field of cell cycle dynamics.The red digit 0 to represent the limit cycle attractor, reflecting cell cycle oscillations, and the black digit 1 to represent the absorption fixed-point attractor, reflecting the stationary phase of the cell at G1. (D) Potential landscape of cell cycle dynamics in UMAP.(E) Curl flux (black arrow) of cell cycle dynamics landscape in UMAP.(F) The gradient force (black arrow) of cell cycle dynamics landscape in UMAP.

Figure S6 .
Figure S6.The differential geometric analysis of cell cycle dynamics of human fibroblasts.(A) Divergence of the reconstructed vector field among all cell cycle phases.(B) Same as in (A) but for the curl.(C) Same as in (A) but for the acceleration.(D) Same as in (A) but for the curvature.

Figure S7 .
Figure S7.The trajectory, autocorrelation and power spectral density of cell cycle oscillation dynamics with different strength of noise.(A) The trajectories along time of the cell cycle for different diffusion coefficient D=0.(B) The trajectories in phase space corresponding to (A) for different diffusion coefficient D=0.(C) The autocorrelation corresponding to (A) for different diffusion coefficient D=0.(D) The power spectral density corresponding to (A) for different diffusion coefficient D=0.(E) Same as in (A) but for different diffusion coefficient D=0.00001.(F) Same as in (B) but for different diffusion coefficient D=0.00001.(G) Same as in (C) but for different diffusion coefficient D=0.00001.(H) Same as in (D) but for different diffusion coefficient D=0.00001.(I) Same as in (A) but for different diffusion coefficient D=0.00005.(J) Same as in (B) but for different diffusion coefficient D=0.00005.(K) Same as in (C) but for different diffusion coefficient D=0.00005.(L) Same as in (D) but for different diffusion coefficient D=0.00005.(M) Same as in (A) but for different diffusion coefficient D=0.0001.(N) Same as in (B) but for different diffusion coefficient D=0.0001.(O) Same as in (C) but for different diffusion coefficient D=0.0001.(P) Same as in (D) but for different diffusion coefficient D=0.0001.

Figure S8 .
Figure S8.The period and amplitude distribution of cell cycle oscillation dynamics with different strength of noise.(A) The distribution of cell cycle period for different diffusion coefficient D=0.(B) The distribution of amplitude in UMAP1 for different diffusion coefficient D=0.(C) The distribution of amplitude in UMAP2 for different diffusion coefficient D=0.(D) Same as in (A) but for different diffusion coefficient D=0.00001.(E) Same as in (B) but for different diffusion coefficient D=0.00001.(F) Same as in (C) but for different diffusion coefficient D=0.00001.(G) Same as in (A) but for different diffusion coefficient D=0.00005.(H) Same as in (B) but for different diffusion coefficient D=0.00005.(I) Same as in (C) but for different diffusion coefficient D=0.00005.(J) Same as in (A) but for different diffusion coefficient D=0.0001.(K) Same as in (B) but for different diffusion coefficient D=0.0001.(L) Same as in (C) but for different diffusion coefficient D=0.0001.

Figure S10 .
Figure S10.The genetic perturbation alters of U2OS cell cycle global dynamics.(A) Cell cycle global dynamics landscape of U2OS cells in UMAP with different genetic perturbation.(i) Suppression of both CCNE2 and CENPE.(ii) Suppression of CENPE only.(iii) Activation of CCNE2 and suppression of CENPE.(iv) Suppression of CCNE2 only.(v) WT. (vi) Activation of CCNE2 only.(vii) Suppression of CCNE2 and activation of CENPE.(viii) Activation of CENPE only.(ix) Activation of both CCNE2 and CENPE.

Figure S11 .
Figure S11.The genetic perturbation alters of RPE1 cell cycle landscape-flux.(A) Cell cycle global dynamics landscape of REP1 cells in UMAP with different genetic perturbation.(i) Suppression of both CCNE2 and KPNA2.(ii) Suppression of KPNA2 only.(iii) Activation of CCNE2 and suppression of KPNA2.(iv) Suppression of CCNE2 only.(v) WT. (vi) Activation of CCNE2 only.(vii) Suppression of CCNE2 and activation of KPNA2.(viii) Activation of KPNA2 only.(ix) Activation of both CCNE2 and KPNA2.(B) EPR of cell cycle nonequilibrium thermodynamics with different genetic perturbation.(C) The average Flux of cell cycle nonequilibrium dynamics with different genetic perturbation.

Figure S12 .
Figure S12.LAPs and MFPT in cell cycle initiation and termination.(A) The correlation between barrier height and least action when diffusion coefficients are changed.The line is a fitting correlation line, and the shaded part is the 95% fitting confidence interval.(B) The correlation between the logarithm of flux and least action when diffusion coefficients are changed.(C) The correlation between the logarithm of EPR and least action when diffusion coefficients are changed.(D) LAPs between different cell cycle phases in the velocity vector field of RPE1 cells.The color of digits in each node reflects the type of fixed point: red, emitting fixed point; black, absorbing fixed point.The color of the numbered nodes corresponds to the confidence of the fixed points.The color of the dots along the paths corresponds to the direction.(E) LAPs between different cell cycle phases in 2D landscape of RPE1 cells.The black arrows present the curl flux and the white arrows present the gradient force in the landscape.(F) LAPs between different cell cycle phases in the velocity vector field of human fibroblasts.The color of digits in each node reflects the type of fixed point: red, emitting fixed point; black, absorbing fixed point.The color of the numbered nodes corresponds to the confidence of the fixed points.The color of the dots along the paths corresponds to the direction.(G) LAPs between different cell cycle phases in 2D landscape of human fibroblasts.The black arrows present the curl flux and the white arrows present the gradient force in the landscape.

Figure S13 .
Figure S13.Phase diffusion and phase coherence of cell cycle oscillation with different noise.(A-C) Raster plot of the peak times for 500 different trajectories starting with the same initial condition at the different diffusion coefficient, D=0 (A), D=0.00005 (B) and D=0.0001 (C).(D-F) Peak time variance σ 2 goes linearly with the average peak time for different diffusion coefficient, with the linear coefficient defined as the peak-time diffusion constant, D=0 (D), D=0.00005 (E) and D=0.0001 (F).(G) The logarithm of phase diffusion and the logarithm of diffusion coefficient D. The line is a fitting correlation line, and the shaded part is the 95% fitting confidence interval.(H) The change of coherence time and phase coherence when diffusion coefficients (fluctuations) are changed.(I) Sketch map for the definition of phase coherence.

Figure S14 .
Figure S14.Jacobian analyses of the cell cycle regulatory interactions in the gene expression space.(A-D) Jacobian analyses of the U2OS cell cycle regulatory interactions in the CCNE2 and CENPE expression space.Activation from CENPE to CCNE2 in the CENPE and CCNE2 expression space (A), repression from CCNE2 to CENPE in the CCNE2 and CENPE expression space (B), (C) self-activation of CENPE in the CENPE and CCNE2 expression space (C), and self-activation of CCNE2 in the CCNE2 and CENPE expression space (D).(E-H) Jacobian analyses of the U2OS cell cycle regulatory interactions in the CCNE2 and KIF23 expression space.Repression from KIF23 to CCNE2 in the KIF23 and CCNE2 expression space (E), activation from CCNE2 to KIF23 in the CCNE2 and KIF23 expression space (F), self-activation of KIF23 in the KIF23 and CCNE2 expression space (G), and self-activation of CCNE2 in the CCNE2 and KIF23 expression space (H).

Figure S15 .
Figure S15.Inference of cell-cycle-phase specific gene regulatory networks for the U2OS cell cycle by single-cell transcriptomics data.(A) The cell-cycle-phase specific gene-gene interaction jacobian matrices of each cell cycle phase.The positive coefficients (red) imply an activation from the regulator gene (x-axis) to the target gene (y-axis), while negative coefficients (blue) imply inhibition.(B) The cell-cycle-phase specific gene regulatory networks.The node size is scaled based on gene expression in the specific cell type and the color of nodes represent the node centrality in the network, the black arrows represent the activation and the red arrows represent the inhibition, and the width of connection arrows is scaled based on the strength of gene interactions.

Figure S16 .
Figure S16.Inference of gene regulatory networks for the RPE1 cell cycle by single-cell transcriptomics data.(A) The interaction network of genes of each cell cycle phase of RPE1 cells.All differentially expressed genes in each cell cycle phase are used to construct network; nodes (genes) with more than one edge are shown.Node colors represent different cycle phases (red: G1-S, yellow: S, blue: G2-M, green: M, purple: M-G1), the black arrows represent the activation and the red arrows represent the inhibition.Genes with black representing TF and red representing non-TFs.(B) The diagram for the cell cycle model with key genes of each cell cycle phase.(C) CCNE2 has high expression in G1-S phase and KIF23 has high expression in G2-M phase.(D) Molecular mechanisms underlying the maintenance of the cell cycle.(i) Repression of CCNE2 by KIF23.(ii) Self-activation of KIF23.(iii) Self-activation of CCNE2.(iv) CCNE2 activates KIF23.(E) Fitting the function of Jacobian versus gene expression with derivatives of a simplistic inhibitory or activation Hill equation.White dashed line corresponds to the zero Jacobian value.The blue stars at each x axis grid point correspond to the weighted mean of the Jacobian values for that point.The blue solid lines are the resultant fittings for the Jacobian.(F) Schematic summarizing the interactions involving CCNE2 and KIF23.(G) The velocity kinetic curves over gene expression changes of the corresponding fitted Hill equations of (E).

Figure S17 .
Figure S17.Inference of cell-cycle-phase specific gene regulatory networks for the RPE1 cell cycle by single-cell transcriptomics data.(A) The cell-cycle-phase specific gene-gene interaction jacobian matrices of each cell cycle phase.The positive coefficients (red) imply an activation from the regulator gene (x-axis) to the target gene (y-axis), while negative coefficients (blue) imply inhibition.(B) The cell-cycle-phase specific gene regulatory networks.The node size is scaled based on gene expression in the specific cell type and the color of nodes represent the node centrality in the network, the black arrows represent the activation and the red arrows represent the inhibition, and the width of connection arrows is scaled based on the strength of gene interactions.

Table S1 .
Differential expression genes for U2OS cell cycle regulation.

Table S2 .
Differential expression genes for RPE1 cell cycle regulation.