Fundamental limits on dynamic inference from single-cell snapshots

A First-Principles Relationship of Cell Dynamics to Static Observations.

When reconstructing a dynamic process from single-cell snapshot data, cells are typically observed in a continuous spectrum of states owing to asynchrony in their dynamics. The goal is to reconstruct a set of rules governing possible dynamic trajectories in high-dimensional gene expression space that are compatible with the observed distribution of cell states. The inferred rules could represent a single curve or branching process in gene expression space, or they could reflect a more probabilistic view of gene expression dynamics. In some cases, multiple time points can be collected to add clarity to the temporal ordering of events. In other cases, a single time point could capture all stages of a dynamic process, such as in steady-state adult tissue turnover.

To develop a framework for dynamic reconstruction from first principles, we want to identify a general, model-independent, mathematical formulation linking cell dynamics to static observations. One possible starting point is the “population balance equation” [also known as the flux balance law (19)], which has the following form:

Eq. 1 states that, in each small region of gene expression space, the rate of change in the number of cells (left-hand side of the equation) equals the net cell flux into and out of the region (right-hand side) (Fig. 1A). The equation introduces the cell density, c(x,t), which is the distribution of cell states from which we sample a static snapshot of cells in an experiment. This density depends on the net average velocity, v(x), of the cells at point x, a feature of the dynamics that we want to infer. Notably, being an average quantity, v is not necessarily a description of the dynamics of any individual cell, but it alone governs the form of the sampled cell density c. Eq. 1 also introduces a third variable: $R$(x) is a rate of cell accumulation and loss at point x caused by the discrete phenomena of cell proliferation and cell death, and by entrance and exit from the tissue being isolated for analysis. Although Eq. 1 is likely a good starting point for analyzing many biological systems, it nonetheless introduces some specific assumptions about the nature of cell state space. First, it approximates cell state attributes as continuous variables, although they may in fact represent discrete counts of molecules such as mRNAs or proteins. Second, it assumes that changes in cell state attributes are continuous in time. This means, for example, that the sudden appearance or disappearance of many biomolecules at once cannot be described in this framework.

Multiple Dynamic Trajectories Can Generate the Same High-Dimensional Population Snapshots.

Given knowledge of the cell population density, c(x,t), we hope to infer the underlying dynamics of cells by solving for the average velocity field v in Eq. 1. This approach falls short, however, because v is not fully determined by Eq. 1, and even if it were, knowing the average velocity of cells still leaves some ambiguity in the specific trajectories of individual cells. This raises the question: Does there exist a set of reasonable assumptions that constrain the dynamics to a unique solution? To explore this question, we enumerate the causes of nonuniqueness in cell state dynamics.

First, assumed cell entry and exit points strongly influence inferred dynamics: For the same data, different assumptions about the rates and location of cell entry and exit lead to fundamentally different inferences of the direction of cell progression in gene expression space, as illustrated in Fig. 1C. Cells can enter a system by proliferation, by physically migrating into the tissue that is being analyzed, or by up-regulating selection markers used for sample purification (e.g., cell surface marker expression). Similarly, cells exit observation by cell death, physical migration out of the tissue being studied, or by down-regulation of cell selection markers. Referring to Eq. 1, this discussion is formally reflected in the need to assume a particular form for the rate field R(x) when inferring dynamics v from the observed cell density c.

Second, net velocity does not equal actual velocity: A second unknown is the stochasticity in cell state dynamics, reflected in the degree to which cells in the same molecular state will follow different paths going forward. A net flow in gene expression space could result from imbalanced flows in many directions or from a single coherent flow in one direction (Fig. 1D). If the goal of trajectory analysis is to go beyond a description of what states exist and make predictions about the future behavior of cells (e.g., fate biases) given their current state, then it is necessary to account for the degree of such incoherence of dynamics.

Third, rotations and oscillations in state space do not alter cell density: Static snapshot data cannot distinguish periodic oscillations of cell state from simple fluctuations that do not have a consistent direction and periodicity (Fig. 1E). As with incoherent motion above, predictive models may need to explicitly consider oscillatory behaviors. The inability to detect oscillations from snapshot data are formally reflected in Eq. 1 by invariance of the concentration c to the addition to v of arbitrary rotational velocity fields u satisfying $\nabla \left(c\mathit{u}\right)=0$.

Fourth, hidden features of cell state can lead to a superposition of different dynamic processes: Stable properties of cell state that are invisible to single-cell expression measurements, such as chromatin state or tissue location, could nonetheless impact cell fate over multiple cell state transitions (Fig. 1F). The existence of such long-term “hidden variables” would clearly compromise attempts to predict the future fate of a cell from its current gene expression state.

Because of these phenomena, no unique solution exists for dynamic inference. However, sensible predictions about dynamics can still be made by making certain assumptions. Our framework for cell trajectory analysis (see below) is based on explicit, reasonable assumptions that together are necessary and sufficient to constrain a unique solution (Fig. 2). These assumptions may nevertheless be inaccurate in certain situations.

Construction of the PBA Framework.

To infer cell dynamics from an observed cell density c, we make the following assumptions.

Reducing to Practice: Solving the Population Balance Equation with Spectral Graph Theory.

Equipped with single-cell measurements and estimates for each fitting parameter, we now face two practical problems in using of Eq. 3 to infer cell dynamics: the first is that Eq. 3 is generally high dimensional (reflecting the number of independent gene programs acting in a cell), but numerical solvers cannot solve diffusion equations on more than perhaps 10 dimensions (23). Indeed, until now, studies that used diffusion-drift equations such as Eq. 3 to model trajectories (10, 13, 20) were limited to one or two dimensions, far below the intrinsic dimensionality of typical scRNA-seq data (24). The second practical problem is that we do not in fact measure the cell density c: we only sample a finite number of cells from this density in an experiment.

Overcoming these problems represents the main technical contribution of this paper. We drew on a recent theorem by Ting et al. (25) in spectral graph theory to extend diffusion-drift modeling to arbitrarily high dimension. The core technical insight is that an asymptotically exact solution to Eq. 3 can be calculated on a nearest-neighbor graph constructed with sampled cells as nodes. Our approach, which we call PBA, actually improves in accuracy as dimensionality increases, plateauing at the underlying manifold dimension of the data itself. We thus avoid conclusions based on low-dimensional simplifications of data, which may introduce distortions into the analysis. In practice, some intermediate degree of dimensionality reduction could still be useful (say, to tens or hundreds of dimensions), a point elaborated in Discussion. SI Appendix provides technical proofs and an efficient framework for PBA in any high-dimensional system.

The inputs to PBA are a list of sampled cell states $\mathit{x}=(x_1,\dots x_N)$, an estimate $R=(R_1,\dots R_N)$ for the net rate of cell accumulation or loss at each state $x_i$, and an estimate for the diffusion parameter $D$. We are assuming steady state, so $\partial c/\partial t=0$. The output of PBA is a discrete probabilistic process, that is, a Markov chain that describes the transition probabilities between the states $x_i$. The analysis is asymptotically exact in the sense that—if a potential exists and the estimates for $R$ and $D$ are correct—the inferred Markov chain will converge to the underlying continuous dynamical process in the limit of sampling many cells ($N\to \infty$) and high manifold dimension of the data (SI Appendix, Theory Supplement, Theorem 4). We note that the requirement for many cells is to reduce the variance of PBA estimates, and the requirement for high dimension is to reduce the bias.

PBA computes the transition probabilities of the Markov chain using a simple algorithm, which at its core involves a single matrix inversion. Briefly:

• i)

  Construct a k-nearest-neighbor (knn) graph G, with one node at each position ${\mathit{x}}_i$ extending edges to the k nearest nodes in its local neighborhood. Calculate the graph Laplacian of G, denoted L.

• ii)

  Compute a potential $V=1/2L^{+}R$, where $L^{+}$ is the pseudoinverse of $L$.

• iii)

  To each edge $(x_i\to x_j)$, assign the transition probability

With the Markov chain generating Eq. 3 available, it is possible to calculate the temporal ordering of states, and the fate biases of progenitor cells in a differentiation process, by integrating across many trajectories (Fig. 2). These calculations are simple, generally requiring a single matrix inversion. Temporal orderings can be calculated from mean first-passage times; fate biases can be calculated as absorbing fate probabilities of the Markov chain, starting from each observed state. Specific formulas are provided in SI Appendix, Theory Supplement, section 3. Confidence intervals on inferred mean first-passage times and fate biases, and their sensitivity to parameter choices, can be estimated by bootstrapping. Code for implementing these and other aspects of PBA is available online at https://github.com/AllonKleinLab/PBA.

PBA Accurately Reconstructs Dynamics of Simulated Differentiation Processes.

We tested PBA on a sequence of simulations, first using an explicit model of diffusion-drift process, and then moving on to direct simulations of gene regulatory networks (GRNs). In the first simulation (Fig. 3 and SI Appendix, Fig. S1), cells drift down a bifurcating potential landscape into two output lineages, defined formally as absorbing states in the dynamical system. Cell trajectories span a 50-dimensional gene expression space (two of which are shown in Fig. 3A). With 200 cells sampled from this simulated system (Fig. 3B), PBA predicted dynamical properties of the measured cells, including (i) their fate probabilities, defined as the probabilities of reaching either of the absorbing states; and (ii) their temporal ordering, defined as the mean first-passage time to reach them from the simulation starting point. PBA made very accurate predictions (Pearson correlation, ρ > 0.96; Fig. 1 C and D) if provided with correct estimates of proliferation, loss, and stochasticity (parameters R and D). Estimates of temporal ordering remained accurate with even fivefold error in these parameters (ρ > 0.93), but predictions of fate bias degraded (ρ > 0.77; SI Appendix, Fig. S2 A–D). Thus, even very rough knowledge of the entry and exit points in gene expression space is sufficient to generate a reasonable and quantitative description of the dynamics. Interestingly, PBA also remained predictive in the presence of implanted oscillations (SI Appendix, Fig. S3; fate probability ρ > 0.9; temporal ordering ρ > 0.8). In addition, the simulations confirmed the theoretical prediction that inference quality improves as the number of noisy genes (dimensions) increases, and as more cells are sampled: maximum accuracy in this simple case was reached after ∼100 cells and 20 dimensions (SI Appendix, Fig. S2 E–G). These simulations showcase the ability of PBA not just to describe continuum trajectories, but additionally to predict cell dynamics and by extension cell fate. At the same time, they show the fragility of dynamic inference to information not available from static snapshots alone.

Having demonstrated the accuracy of PBA on an explicit model of a diffusion-drift process, we next tested its performance on gene expression dynamics arising from GRNs (Fig. 4). As before, we simulated cell trajectories, obtained a static snapshot of cell states, and supplied PBA with this static snapshot as well as the parameter $R$ encoding the location of entry and exit points. We began with a simple GRN representing a bistable switch, in which two genes repress each other and activate themselves (Fig. 4A). Simulated trajectories from this GRN begin with both genes at an intermediate expression level, but quickly progress to a state where one gene dominates the other (Fig. 4B). In addition to the two genes of the GRN, we included 48 uncorrelated noisy dimensions. With 500 cells sampled from this process, PBA predicted cell fate bias (absorbing state probabilities) and temporal ordering (mean first passage time from the start point) very well (r > 0.98 for fate bias and r > 0.89 for ordering; Fig. 4C), although the precise accuracy depended on the assumed level of diffusion $D$ (SI Appendix, Fig. S4; the values that gave the best results were used for Fig. 4).

PBA assumes the absence of oscillations in gene expression space. Therefore, it is unclear how well PBA can infer cell trajectories that result from GRNs with oscillatory dynamics. We simulated an oscillatory GRN in the form of a “repressilator” circuit (26) with the addition of positive-feedback loops that create two “escape routes” leading to alternative stable fixed points of the dynamics (Fig. 4D). Simulated trajectories from this GRN begin with all genes oscillating, followed by a stochastic exit from the oscillation when one of the genes surpasses a threshold level (Fig. 4E). With 500 cells sampled from this process, PBA was significantly less accurate than for the previous simulations (Fig. 4F). Although PBA correctly identified which cells were fully committed to the two “escape routes,” it was entirely unable to resolve the fate biases of cells in the uncommitted oscillatory state. PBA also made poor predictions of mean first-passage time, underestimating the amount of time that cells spent in the oscillatory state. Unsurprisingly, when the assumptions of PBA are strongly violated, its prediction accuracy suffers.

PBA Predictions of Fate Bias in Hematopoiesis Reconcile Past Experiments.

To test PBA on experimental data from real biological systems, we made use of single-cell gene expression measurements of 3,803 adult mouse HPCs from another study by our groups (ref. 16; data at Gene Expression Omnibus, accession no. GSE89754).

HPCs reside in the bone marrow and participate in the steady-state production of blood and immune cells through a balance of self-renewal and multilineage differentiation. Descriptions of HPC differentiation invoke a tree structure, with gradual lineage-restriction at branch points. However, the precise tree remains controversial (27, 28), since existing measurements of fate potential reflect a patchwork of defined HPC subsets that may have internal heterogeneity (29) and provide only incomplete coverage of the full HPC pool. We asked whether PBA applied to single-cell RNA profiling of HPCs could generate predictions consistent with experimental data, and possibly help resolve these controversies by providing a global map of approximate cell fate biases of HPCs.

The single-cell expression measurements—derived from mouse bone marrow cells expressing the progenitor marker Kit—represent a mixture of multipotent progenitors as well as cells expressing lineage commitment markers at various stages of maturity. Since PBA depends on analyzing a k-nearest-neighbor (knn) graph of the cells, we developed an interactive knn visualization tool for single-cell data exploration, called SPRING (https://kleintools.hms.harvard.edu/tools/spring.html; ref. 30). The SPRING plot (Fig. 5A) revealed a continuum of gene expression states that pinches off at different points to form several downstream lineages. Known marker genes (SI Appendix, Table S1) identified the graph endpoints as monocytic (Mo), granulocytic (G), dendritic (D), lymphoid (Ly), megakaryocytic (Mk), erythroid (Er), and basophil or mast (Ba/Mast) cell progenitors (SI Appendix, Fig. S5); we also identified cells in the graph expressing hematopoietic stem cell markers. The lengths of the branches reflect the timing of Kit down-regulation and the abundance of each lineage.

For steady-state systems, PBA requires as fitting parameters an estimate of the diffusion strength $D$ and the net rates of cell entry and exit at each gene expression state ($R)$. We estimated R using prior literature (Materials and Methods) and tested a range of values of D, with the range chosen to ensure that the number of cells predicted to be multipotent by PBA would lie within bounds established in the literature. All results that follow hold over the physiological range of PBA parameter values (SI Appendix, Fig. S6).

Applied to single-cell measurements, for the range of fitting parameters, PBA estimated seven fate probabilities for each cell, defined formally as the probabilities that a trajectory initiated at each cell would terminate among the most mature cells in each of the seven lineages (Materials and Methods). We compared PBA results to fate probabilities reported for subsets of the HPC hierarchy, which have been previously defined by cell surface marker expression, and transcriptionally profiled using microarrays. To carry out the comparison, we identified the cells in our data that were most similar to each subset using their published microarray profiles (Fig. 5B, red dots), and we then computed the average PBA-predicted fate probabilities across cells in each subset. Remarkably, for a panel of 12 progenitor cell populations from six previous papers (31–36) (SI Appendix, Table S2) the PBA-predicted fate outcomes (Fig. 5B, bar charts) closely matched fate probabilities measured in functional assays (defined as the proportion of clonogenic colonies containing a given terminal cell type; SI Appendix, Fig. S5). The main qualitative disagreement between PBA predictions and experiment was in the behavior of Lin⁻Sca1⁻Kit⁺IL7R⁻FcgR^lowCD34⁻ HPCs, previously defined as megakaryocyte-erythroid precursors (31). Our prediction was that these cells should lack megakaryocyte potential, which is indeed consistent with recent studies (27, 29, 37). Excluding these cells, our predicted fate probabilities matched experimental data with correlation ρ = 0.91 (Fig. 5C). In another study (16), we test several novel predictions in hematopoiesis emerging from PBA.

Experimental Design for Trajectory Reconstruction from Static Snapshot Measurements.

We have shown that accurate dynamic inference requires knowledge of the density of cells in high-dimensional state space, as well as the rates of cell entry and exit across the density. These requirements immediately suggest a set of principles for experimental design to optimize dynamic inference. First, to minimize distortions in the cell density in gene expression space, it is useful to profile a single, broad population, and to avoid merging data from multiple subpopulations fractionated in advance. Where possible, one should consider testing for uniform cell-sampling probability across states. Second, if cells of interest are sorted before analysis, it is best to minimize the number of sorting gates and enrichment steps, since each introduces an additional term to the entry/exit rates and subsequently a risk of distortion to the inferred dynamics. The HPC dataset analyzed in this paper was well suited for trajectory reconstruction because it consisted of a single population, enriched using a single marker (Kit). This contrasts with previous scRNA-seq datasets of hematopoietic progenitors that included a composite of many subpopulations (38) or used complex FACS gates to exclude early progenitors (29).

How PBA Could Go Wrong.

To constrain a unique solution for trajectory reconstruction, PBA makes several strong assumptions, such as memoryless dynamics with respect to measured states, the absence of oscillations in gene expression space, and an adequately “large” number of sampled cells.

Normalization, Principal-Components Analysis, and Other Coordinate Transformations.

In this study, we described a framework for modeling the movement of cells in a space of gene expression, the units of which might be considered to be (dimensionless) counts of individual molecules. How then should one think about routine transformations of gene expression coordinates performed during practical low-level processing of single-cell expression data, such as transformation into logarithmic space, or dimensionality reduction by principal-components analysis? Here, the asymptotic analysis of PBA makes clear that coordinate transformations may not be important when cells are densely sampled, as they should leave the empirical single-cell graph topology unchanged. The equations of PBA are indeed invariant to coordinate transformations, with the exception of the diffusion operator, which is isotropic and spatially homogeneous but may not remain so upon coordinate system transformation. Since our assumption of isotropic and invariant diffusion is already an approximation, it does not support a priori one coordinate system over another. For small and noisy datasets, the choice of coordinate system could affect conclusions, however, and it is probably best to use the coordinate system that provides the richest view of single-cell population structure, or that agrees most with known biology.

A related question is whether technical noise in gene expression measurements can be distinguished from biological effects. For scRNA-seq methods that use techniques such as linear amplification and unique molecular identifiers, technical noise can mostly be attributed to undersampling of the mRNAs in a cell, and therefore follows a predictable Poisson distribution determined by the underlying gene expression and technical sampling rate. Distinguishing technical noise using known features of this distribution would be an area for improvement.

Fundamental Limits on the Inference of Gene Regulatory Networks.

One promise of single-cell expression measurements is their possible use for reconstructing GRNs (2, 11). However, since any GRN model entails specific hypotheses about the gene expression trajectories of cells, efforts to infer GRNs from single-cell data must also confront the limits of knowledge identified in our framework. In particular, GRN inference may benefit from an explicit consideration of cell entry and exit rates (embodied by R) and the rate of change in the cell density ($\partial c/\partial t$), as well as acknowledging the inability to distinguish oscillations from fluctuations.

Indeed, the inability to detect oscillations in single-cell data, embodied in our framework by the use of a potential landscape, suggests severe limits on the types of underlying gene regulatory relationships that can be modeled. In fact, potential landscapes can only emerge from GRNs with strictly symmetric interactions, meaning every “arrow” between genes has an equal and opposite partner. This result follows from observing that the arrows in a GRN describe the influence of gene $i$ on gene $j$, which is given by $\partial {\mathit{J}}_i/\partial x_j$ (Fig. 6A), where $\mathit{J}$ is the deterministic component of average cell velocities (Eq. 2). The assumption of a potential landscape (i.e., $\mathit{J}=-\nabla F)$ then imposes symmetry on the GRN because $\partial J_i/\partial x_j=\partial J_j/\partial x_i$ = $-{\partial }^{2}F/\partial x_i\partial x_j$. Although a few well-known GRN motifs follow this symmetry rule—such as the “bistable switch” resulting from the mutual inhibition of two genes—many others do not, such as negative-feedback loops and oscillators (Fig. 6B). Potential landscapes are frequently invoked to explain gene expression dynamics (10, 41, 42), and we have shown them to be useful for predicting HPC fate outcomes in the context of PBA. It seems paradoxical that a tool that provides realistic phenomenological descriptions of gene expression dynamics reflects an entirely unrealistic picture for the underlying gene regulatory mechanisms. Resolving this paradox is an interesting direction for future work.

How Should We Visualize Single-Cell Data?

At its core, the PBA algorithm performs dynamic inference by solving a diffusion-drift equation in high dimensions. This computation relies on a 2011 result in spectral graph theory by Ting et al. (25) that describes the limiting behavior of k-nearest-neighbor graph Laplacians on sampled point clouds. Interestingly, several recent studies (8, 15, 43) have developed k-nearest neighbor graph-based representations of single-cell data, and others have suggested embedding cells in diffusion maps (24, 44) on the basis of other similarity kernels. It has been unclear, until now, how to evaluate which of these different methods provides the most useful description of cell dynamics. Our technical results (SI Appendix, Theorems 1–4 in SI Appendix, Theory Supplement) confirm that certain graph representations provide an asymptotically exact description of the cell state manifold on which dynamics unfold, suggesting them to be useful techniques for visualizing single-cell datasets. Therefore, PBA formally links dynamical modeling to choices of single-cell data visualization.