Comparison of cell state models derived from single-cell RNA sequencing data: graph versus multi-dimensional space

Single-cell sequencing technologies have revolutionized molecular and cellular biology and stimulated the development of computational tools to analyze the data generated from these technology platforms. However, despite the recent explosion of computational analysis tools, relatively few mathematical models have been developed to utilize these data. Here we compare and contrast two cell state geometries for building mathematical models of cell state-transitions with single-cell RNA-sequencing data with hematopoeisis as a model system; (i) by using partial differential equations on a graph representing intermediate cell states between known cell types, and (ii) by using the equations on a multi-dimensional continuous cell state-space. As an application of our approach, we demonstrate how the calibrated models may be used to mathematically perturb normal hematopoeisis to simulate, predict, and study the emergence of novel cell states during the pathogenesis of acute myeloid leukemia. We particularly focus on comparing the strength and weakness of the graph model and multi-dimensional model.


Introduction
The ability to apply genome sequencing methods to single-cells has revolutionized biology [1]. Technologies enabling single-cell sequencing are advancing rapidly, with datasets as large as hundreds of thousands of cells are common [2]. RNA-sequencing is currently the most prevalent form of single-cell genomic analysis [1]. The sequencing of RNA at the cellular level enables the interrogation of gene transcription, which is used as a high-dimensional fingerprint which characterizes the identity of the cell. For this reason, single-cell RNA-sequencing data (scRNA-seq) has been used as a tool to study cell identity and state-transitions at the cellular level.
The most frequently studied cell state-transition is cellular differentiation; the process of a cell and its progeny to perform specialized tasks through transformation from a less differentiated stem-like state to a more differentiated state. ScRNA-seq is used to identify cells in various states of differentiation primarily through one or both of two primary methods: 1) clustering of cells with similar features [3], or 2) though trajectory inference (TI) [4]. Clustering analysis relies on the definition of a similarity metric, and may rely on a pre-defined number of clusters (e.g., k-means), or may use optimization methods to identify clusters (e.g., Leiden). There are a wide variety of clustering methods and similarity metrics to choose from, which may give drastically different results [5]. Similarly, trajectory inference methods may use pre-defined relationships between cells or may use optimization methods to identify these relationships to construct graphs which are then used to infer paths, or trajectories, between cell states. In addition, various approaches aim to characterize the cell fate landscape, for instance, by a parameterized landscape based on bifurcation analysis [6,7], by using a measure of entropy of cell states: SCENT [8] and scEpath [9], or by mapping cells to a landscape on optimized parameters: HopLand [10] and Topslam [11].
A significant limitation of these approaches is if the graph structure and underlying relationships between the cells is unknown. As shown in a comprehensive review of trajectory inference methods by Saelens et al. (2019) [4], most TI algorithms have difficulty inferring even simple graphs which may include cycles or disconnected subgraphs. Because of the limitations of clustering and trajectory inference in analysing these data, we suggest that a hypothesis-driven and mathematical approach to the analysis of scRNA-seq data to study cell state transitions is warranted.
Moreover, single-cell genomic sequencing suffers from a number of challenges in analysis. Beyond the several choices to be made for even simple analyses such as clustering or visualization, the data may be sparse and incomplete. Gene "drop outs" and background signal (noise) can complicate differential expression and clustering analyses. For this reason, analysis of these data have remained fairly superficial despite the wealth of information contained in these high-dimensional datasets. Moreover, results obtained from analysis of single-cell sequencing datasets may be very sensitive to choices in the method of analysis and algorithm parameters. To date, single-cell sequencing data have not been effectively leveraged as inputs into mathematical models.
Here we compare two cell state geometries of cell state-transitions modeling with scRNAseq data. Building on our prior work [12], we model cell differentiation as a continuous process. To elaborate this concept, when cell type-A becomes cell type-B, the cell states during the transition process are often classified into more steps as type-A½ or types-A¼, , A¾. The continuous cell states can be considered as a limit of these states. We develop phenotype structured cell state models assuming continuous cell states using reaction-diffusion-advection partial differential equations (PDE) solved on: (i) an abstracted graph and (ii) a multi-dimensional continuum space. We compare and contrast these two cell state geometries with hematopoeisis as a model system. This manuscript is structured as follows: first we present the PDE model on a graph and in continuous space, then we apply the model to two datasets, see [13,14]. We examine the impact of various graph construction and trajectory inference methods on the geometry of the cell state space, and solve the model on these geometries. We then use the model the study the effects of perturbing 1) the graph structure 2) expression of select subsets of genes 3) and cell state transition dynamics by perturbing flow on the graph or by modifying the dynamics in the continuous space. We predict novel dynamics of leukemia pathogenesis by perturbing normal hematopoeisis and conclude with a comparative analysis of our approach and description of future directions for mathematical modeling with single-cell genomic sequencing data. A summary of our workflow is shown in Figure 1.

Modeling cell state-transitions in a continuous cell state-space
In this section we develop PDE models of cell dynamics in the continuous phenotype space identified by dimension reduction techniques. For a given single-cell genomic sequencing dataset g i i = 1 N , g i = g 1 i , g 2 i , …, g G i ∈ ℝ G , where N is the number of cells and g i is a G-dimensional vector of gene expression of the i-th cell, the dimension reduction method can be written as an operator P: ℝ G Γ ⊂ ℝ n where the reduced dimensional space is truncated at the n-th dimension and n ≪ G. We denote the reduced space variable as θ = P(g) ∈ Γ ⊂ ℝ n , θ = θ 1 , θ 2 , …, θ n , (2.1) and the i-th single-cell data can be transformed into the reduced space as P g i = θ i = θ 1 i , θ 2 i , …, θ n i . Various dimension reduction techniques exist to construct a mapping P, including principal component analysis (PCA), diffusion mapping, and tdistributed stochastic neighbor embedding (t-SNE). While different techniques provide different shapes and differentiation spaces, we choose diffusion mapping due to its ability to capture non-linear structure of high-dimensional data, and to well reproduce global trajectory of data [15]. We comment that if the reduced dimensional space is not clear to truncate at a low-dimension, one can consider low-dimensional marker genes according to the cell state of interest, and semi-supervised learning approach can be applied to obtain the low-dimensional reduced space. We also comment that it is common to remove the effect of cell cycles from the gene expression data to eliminate the cell state regarding their location along the cell cycle [16].

PDE model of cell state-transitions on a multi-dimensional reduced component space
We first develop a cell state model that describes the dynamics of cell distribution u(t, θ) on the reduced component space Γ, where θ ∈ Γ is the variable that represents continuous cell state. Three highly distinctive dynamic regimes of the cell states are considered, namely, directed cell transition, birth-death process, and random phenotypic instability. Such model can be written as an advection-reaction-diffusion PDE that governs the cell distribution u(t, θ) as ∂ t u(t, θ) = − ∇ ⋅ (V (t, θ)u(t, θ)) + R(θ, u(t, θ)) + ∇ ⋅ (D(θ) ∇u(t, θ)), (2.2) with zero Dirichlet boundary condition. The three terms in our equation that involve parameterized functions V, R, and D, represent cell differentiation, population growth, and phenotypic instability, respectively.
Let us first describe the advection term V ∈ ℝ n that represents directed cell differentiation, where we propose two candidates for modeling V, denoted as v 1 and v 2 . The first candidate v 1 assumes an attractor cell states of homeostasis. Assuming that the magnitude of phenotypic instability is with a magnitude v, that is, D(θ) = ν, one can compute the advection term v 1 as v 1 (θ) = ν ∇ θ U(θ), (2.3) where U(θ) can be computed from the homeostasis distribution u s (θ) that can be regarded as the cell landscape that the hematopoiesis system desires to maintain. As in the Boltzmannlike distribution from equilibrium statistical mechanics [17], we compute U(θ) as the exponent of u s (θ) in the exponential form, in other words, U θ = − ln u s θ . There are multiple methods to compute the cell landscape, so called quasi-potential, that focuses on relative stability of multiple attractors and models cell differentiation as transition between the attractor states [6][7][8][9][10][11]. Here, we compute the cell landscape empirically by assuming that the entire single-cell data is a representative subset of the entire hematopoiesis system, and by using density approximation methods. In particular, we use kernel density estimation [18] from the projected single-cell data θ i ∈ Γ, i.e., u s θ = 1 N i = 1 N K ℎ θ − θ i where we chose the standard normal density function as the kernel function K h with bandwidth parameter h > 0.
The second candidate v 2 models the dynamics of cell state transition. We model this term using a mechanistic approach that describes the symmetric and asymmetric cell division of stem cells to more differentiated cells. In particular, we consider the following form v 2 (t, θ) = c(θ)[2(1 − a(θ))r(θ)]s(t), (2.4) that is parameterized by the proliferation rate r(θ) and the self renewal rate a(θ) [19]. c θ ∈ ℝ n represents the direction and magnitude of differentiation on the phenotype space that we can estimate with either temporal data or pseudotime inference methods [4]. We note that the self renewal rate a(θ) is the proportion of cells that remains in cell state θ, while 1a(θ) cells further differentiate into matured states. This can be counted from symmetric and asymmetric stem cell division. In addition, we assume a signal parameter s(t) that controls the active differentiation term, where s(t) = 1/(1 + km(t)) and m(t) is the number of matured cells. Finally, we comment that the directed cell transition is simulated as V = v 1 + v 2 , that is a sum of cell transition to attain homeostasis and active cell differentiation.
The reaction term represents the growth rate that consists of proliferation and apoptosis. We comment that the calibration of this term requires additional data to scRNA-seq, particularly, the population level growth data, to uniquely calibrate the model. It has been studied that the cell dynamics cannot be uniquely determined without imposing the reaction term [20]. More recently, there has been efforts to estimate the proliferation rate directly from scRNA-seq data by cellular barcoding techniques [21]. In our simulations, we cluster the single-cell data into biologically well known cell types, for instance, in case of hematopoiesis, myeloid progenitors, lymphoid progenitors, macrophages, and obtain the proliferation rate and self renewal rate from the literature. We consider the logistic growth term as following R(θ, u) = r(θ)(1 − d(θ, u))u, (2.5) where r(θ) is the proliferation rate and d(θ) is the apoptosis term assuming a logistic growth as d θ, u = min u u s θ , d , where d models the maximum magnitude of apoptosis rate.
The second-order diffusion term represents the instability on the phenotypic landscape of the cells that should be taken into account when modeling the macroscopic cell density. We simply consider a constant term D(θ) = v. Assuming that the cell state trajectory is subject to Gaussian white noise, the diffusion coefficient can be estimated as the variance of the cell trajectory θ(t) on the reduced space, ν = Var(θ(t))/4. However, since we do not have the data of cell trajectories, one can estimate the value as a limit of random walk as ν = (Δx 2 )/(4Δt), assuming that Δx is the step size of the phenotypic fluctuation in Δt time [22]. See Appendix A for the detail of the model.

PDE model of cell state-transitions solved on a graph
Although the continuum-based multi-dimensional model provides a framework to study cell states, it is not always straightforward to map back locations in the space to novel or otherwise unknown cell states. Moreover, a central feature of contemporary analysis of scRNA-seq data is clustering and inferring relationships between clusters of known cell types [4]. Therefore, we develop a model that can describe cell state-transition dynamics on a graph that represents relationships between known cell types identified with clusters, extended from our previous work in [12]. An immediate advantage of this cell state geometry is that it is convenient to employ biological insights from well-known classical discrete cell states.
The continuum of differentiation cell states is assumed to be on the graph obtained from the scRNA-seq data, for instance, using partition-based graph abstraction (PAGA) algorithm [23]. We project the graph on the reduced component space, and denote the nodes as v k k = 1 n v and the edges as e ij connecting in the direction from the i-th to the j-th node. For convenience of notation, the edges are also denoted as e k k = 1 n e with the index mapping The model follows the dynamics of the cell distribution on the graph, u(x, t), where x ∈ e k is a one-dimensional variable that parameterizes the differentiation continuum space location along the edges. We annotate the cell distribution on each edge e k as u k (x, t) such that u(x, t) = u k (x, t) k = 1 n e , and model the cell density by an advection-reaction-diffusion equation [24] as The three terms are similarly modeled as the multi-dimensional model in Eq (2.2), representing cell differentiation, population growth, and phenotypic instability. To summarize once more, the advection coefficient V k (x) models the cell differentiation and the transition between the nodes, that is, different cell types. We model the advection term in two parts as in the reduced component space model, V k (x) = v k,1 (x) + v k,2 (x), and compute them as v k, 1 (x) = v ∂ x U k (x), v k, 2 (t, x) = 2 1 − a k (x) r k (x) s(t) . (2.7) Here, u s,k (x) = e −U k (x) is the homeostasis cell distribution on the k-th edge, ν is the magnitude of phenotypic instability, r k (x) is the proliferation rate, a k (x) is the self-renewal rate, and s(t) is the signal parameter. Cell proliferation and apoptosis can be modeled by the reaction coefficient R k (x) as Finally, the diffusion term that represents phenotypic fluctuation is taken as D k (x) = ν. i, n ∈ J ℬ I i, n u, b I i, n = n, j ∈ J ℬ I n, j u, a I n, j , u b I i, n = u a I n, j , for all i, n ∈ J, n, j ∈ J, (2.9) where

Quantification of cell state-transition dynamics
Let us define some useful quantities to interpret model predictions in the multi-dimensional cell state-space and on a graph. The total number of cells from the cell distribution on either a graph or a continuous manifold can be computed as (2.10) respectively. In addition, we compute the number of cells of a specific cell type by assigning a weight, w k , that corresponds to cells in the k-th cluster as ρ k (t) ≐ Γ u(t, θ)w k (θ)dθ, (2.11) with ∑ k w k (θ) = 1. We assign weights based on the relative cell density of each clusters estimated with kernel density estimation. In the graph model, we assign the cell states along the edge to be the cell type of the closest node, so that we can compute the number of the k-th node cell type as ρ k (t) ≐ Although we can understand the continuum of cell states by mapping cells in intermediate states back to known discrete cell types, we also desire to interpret the continuous cell states in their location without reference to the canonical cell identities. For such purpose, we characterize cell states by identifying genes that are strongly correlated to a location in the space or moving in a direction toward a cell state. This extends finding the genes that are correlated to specific reduced space components to analyze the reduced cell state space [25]. First, to characterize the cell state θ* in the reduced space, we use a function f θ* (θ) centered at θ* as f θ * (θ) = 1 2πσ exp − θ − θ * 2 /2σ 2 , and compute the correlation between the function values and the j-th gene expression levels as r f, j ≐ corr f, g j , (2.12) where f represents the vector of function evaluation at each single-cell data point θ i , that is, N . An alternate quantity to examine is the genes that are related to a certain direction in the reduced component space. For instance, the correlation between the j-th gene and the k-th reduced component θ k = θ k i i = 1 N and to a certain vector v = v k i = 1 n can be computed as r k, j ≐ corr θ k , g j , r v, j ≐ k v k r k, j , (2.13) respectively. Regarding Eq (2.13) as global quantities, we can also compute the local correlation on the subdomain of the reduced space Ω d by collecting the cell indices that lie within the subdomain Γ d = i θ i ∈ Ω d , that is, r k, j | Γ d ≐ corr((θ k , g j )| i ∈ Γ d ) and Although these metrics may provide candidates of potential genes that are related to the cell state to be analyzed, we emphasize that these need to be verified experimentally by observing the cell state change after perturbing the genes. See Section 4.2 for the limitations of this approach.

Simulation of continuous cell state models on multi-dimensional space versus graph
In this section, we employ the framework developed in the previous section to the mouse hematopoiesis cell data from Nestorowa [14]. We obtain the graph and multi-dimensional space models of hematopoiesis cell state and focus on comparing the strengths and weaknesses of the two models.
The hematopoiesis single-cell data from [13,14] projected on the first two diffusion component space are shown in Figure 2A, where distinct cell types, including lymphoid-primed multipotent progenitors (Lymph); common myeloid progenitors (CMP); megakaryocyte-erythroid progenitors (MEP); granulocyte-macrophage progenitors (GMP); erythrocytes (Ery); neutrophils (Neu); monocytes (Mo); megakaryocytes (Mk); basophils (Baso), classified in the original papers are illustrated with different colors. We truncate the diffusion component at two since the reduced two-dimensional space describes the dynamics of our interest, that is, from strong to weak stemness. The first two diffusion components θ 1 and θ 2 represent cell maturation in both data sets. In Nestorowa data, the first diffusion component separates the stem cells to myeloid lineages, particularly MEP cells and the second diffusion component describes GMP cells and the lymphoid lineages. In Paul data, the first and second reduced component represents MEP and GMP lineage, respectively. We remark that the cells that are the most stem-like in Paul data are CMPs, that is more matured than the ones in Nestorowa data, that includes the long-term and short-term HSCs. In addition to the single-cell data, the Figure 2B presents the abstracted graphs obtained using PAGA [23]. Further refinements of the graph will eventually become the full single-cell data, where each single-cell being counted as distinct cell states, and it depicts the hierarchy of cell states (see Figure A5).
The homeostatic cell distribution u s on the reduced dimensional space is computed by kernel density approximation [18,26]. The computed cell landscapes viewed in different angles are shown in Figure 2C,D. The cell distribution on the graph can be similarly obtained after reallocating the cells to the node, that is, the center of each cluster. The cell distribution on the continuum space provides an intuitive method to compare the relative concentration of different cell lineages, including the intermediate cell states. We observe high concentration of MEP and Ery cells that are localized at the far right (large θ 1 ) in both data. The Nestorowa data has more diverse cells including the common lymphoid progenitors that are visible on the left (small θ 1 , and intermediate θ 2 ), while the Paul data has evenly distributed cell states among the most stem-like cells (CMP) and the two different lineages.
Let us summarize the properties of the graph and multi-dimensional space models before we present simulation results. The graph model has its strength that the cell lineages between the known cell states can be more easily identified as compared to the multidimensional space model. The cell concentration moving toward different edges can be clearly distinguished as the cell lineages to different cell states. Accordingly, counting the number of cells in each discrete cell state is more straightforward, for instance, by integrating the cell distribution along the edges half way. Although the multi-dimensional space model has ambiguity on classifying the cells into discrete cell types, the number of cells in each discrete cell type can be computed by assigning weights to integrate as in Eq (2.11). Moreover, we emphasize that the advantage of clear cell states in the graph model is also its limitation at the same time, since it restricts the model to only study the known cell types and lineages. The advantage of the multi-dimensional space model is its potential of exploring novel cell states that deviates from known cell types. While the graph model cannot explore the cell states that are not already included in the graph structure, the multidimensional space model can immediately study the abnormal trajectories and emergence of cells at any space location. We will show later in our simulation that the hypothesis of genetic alterations can be studied directly in the multi-dimensional space model, without projecting it on the graph structure. Moreover, the multi-dimensional space model is more sensitive to genetic variations than the graph model, although when the variation is large and a considerably distinct cell state arises, the graph model can append another cluster node. See Table 1 for a summary.
In the following sections, we consider mainly two application problems, namely, normal hematopoiesis and abnormal hematopoiesis differentiation, resulting in myeloid leukemia as application examples of our modeling approach.

Calibrating the mathematical models to normal hematopoiesis
We demonstrate that normal hematopoiesis can be visualized by both models on the graph and on the space of two-dimensional diffusion components, (see Figures A1 and A2 for the advection and reaction terms used in the multi-dimensional space model). Figure 3 shows a cluster of stem cells differentiating into the entire cell states on the graph and reduced space using Nestorowa data [13] and Paul data [14]. The initial condition is imposed as approximately 10% of cell capacity in normal condition mostly composed with stem cells.
On both graph and multi-dimensional space models, nontrivial amount of most matured cell states, particularly, Ery and Neu/Mo cells arise around pesudotime t = 30, and further recovers the full landscape after approximately t = 50. In particular, the observation that the matured cells quickly proliferate to fill in the space agrees in both simulations from Nestorowa and Paul, while the effect is more significant in Paul's data due to shorter distances of the matured cell states from the initially administered cells.
The advantage of the graph model is apparent that we can observe distinct cell states as a mass of cells shifting from a node to distinct edges toward different cell states. For instance, the cells differentiating from the MPP cluster to either Neu/Mo lineage and Ery lineage can be clearly separated in the graph models, while those can be ambiguous in the two-dimensional distribution. Still, we can compute the number of cells in each cell types in both models as shown in Figure 3B. We observe that the number of cells reaches the maximal capacity at later times around t = 100, with the intended ratio of cell numbers in each discrete cell type approximating the given data [13] in Figure 3C. The recovery is more rapid for larger values of ν and larger number of initial stem cells ρ(0) (see Figure A6). We remark that the continuous cell states of hematopoeisis is also depicted in conventional flow cytometry which is typically used to identify distinct cell populations based on expression of cell surface markers. We performed Fluorescence Activated Cell Sorting (FACS) analysis of bone marrow cells isolated from normal C57Bl/6 mice (age 6-8 weeks). Distinct myeloid progenitor types (CMP, GMP and MEP) are typically differentiated by the expression of CD16/32 and CD34 markers within the myeloid lineage progenitor cell compartment. Figure  3D shows representative FACS data with respect to CD16/32 and CD34 expression that is used to identify the CMP, GMP, and MEP cell types within the normal myeloid progenitor compartment. Although the FACS data is conventionally clustered and gated into three cell types, continuity of CD16/32 and CD34 expression can be observed that agrees with our graph abstraction and multi-dimensional cell state geometries.

Using the model framework to simulate acute myeloid leukemia (AML) pathogenesis and progression
In this section, we once more compare the graph and multi-dimensional space models with an application to abnormal differentiation under leukemia pathogenesis and progression. We first consider AML model in the context similar to the previous section that involves known progenitor cells that exemplifies the advantage of graph models. However, we will show how aberrant differentiation and phenotypic plasticity of leukemia pathogenesis motivates the spatial model.
AML results from aberrant differentiation and proliferation of transformed leukemiainitiating cells and abnormal progenitor cells. We model AML pathogenesis based on known behavior of a genetic Cbfb-MYH11 (CM) knock-in mouse model that recapitulates somatic acquisition of a chromosomal rearrangement, inv(16)(p13q22) [28,29], commonly found in approximately 12 percent of AML cases. Inv(16) rearrangement results in expression of a leukemogenic fusion protein CBFβ-SMMHC, which impairs differentiation of multiple hematopoietic lineages at various stages [30][31][32]. Most notable in such leukemia pathogenesis and progression is the increased in abnormal myeloid progenitors, which has an MEP-like immunophenotype and a CMP-like differentiation potential [31]. Experimental studies [27,33] show that such MEP attains a predominant increase in pre-megakaryocyte/ erythroid (Pre-Meg/Ery) population (ranging from 5 to 12 fold) accompanied by impaired erythroid lineage differentiation as shown in Figure 4E. This refined phenotypic Pre-Meg/Ery population consists partly of the CMP and MEP populations which are identified using conventional markers [13,34].
In our model, abnormal leukemic progenitors are regarded as intermediate cell states along the edges connecting CMP (or MPP) and MEP (and Ery) in the graph model, and the corresponding locations in the multi-dimensional space model. We assume a 10-fold increase on average in those population by lowering d(θ) and d k (x) in Eqs (2.5) and (2.8) that controls the local cell capacity. In addition to over-proliferation, another aspect of the leukemia pathogenesis of our interest is the impaired differentiation of erythroid lineage differentiation, where it can be modeled by blocking the cell differentiation V(θ) in Eq (2.6) and V k (x) in Eq (2.2) toward Ery.
The corresponding results are shown in Figure 4, where we modify the model to leukemia progression at t = 10. The cell distribution changes from the normal cell landscape at t = 10 to an increased population of Ery (MEP) and nearby cells at t = 20 in both graph and continuum models. We observe a 10-fold increase in the Ery (MEP) population, which includes the abnormal myeloid progenitors, in both graph and multi-dimensional space model across the data sets as shown in Figure 4D,E. The rapid emergence of AML occurs within two week period, corresponding to the expansion of the leukemic cell phenotype. We observe a rise in the MPP or CMP cluster as shown in the results from Nestorowa data, that is similar in Paul data as well. The total proportion of leukemic cells comprise 50-60% of the total population.
In the leukemia pathogenesis simulation in this section, we focus on studying the leukemic cells as a variation of cell states classified using conventional markers. In this case, the graph model can interpret and include the dynamics of such cells, as well as the multi-dimensional space model. While the simulation outcome between the graph and multi-dimensional space model is similar, the graph model is computationally more efficient due to the fewer number of unknowns as compared to the two-dimensional space model. However, to study the abnormal cell states that may appear far away from the known or existing landscape, we will show in the following section that the multi-dimensional space model has more freedom to include those new cell states and disrupted trajectories. We will study the impact of perturbation of genes in the graph and multi-dimensional space model, particularly focusing on alterations of genes known to be involved in leukemia pathogenesis.

In silico experiments of gene expression perturbation
In this section we investigate the sensitivity of altering specific genes in a prescribed manner and the impact of this perturbation in the graph and multi-dimensional space models. We keep our focus on leukemia pathogenesis and progression and consider alterations of 38 genes that are reported to be related to leukemia stem cells [35,36], although we emphasize that these genes serve simply as examples, and are not intended to model the precise biological process of AML pathogenesis. The j-th gene expression level of i cell, g j i , is modified as g j i = 2 γ j g j i , 0 ≤ log 2 g j i + 1 ≤ 16, (3.1) where γ j is the log 2 -fold change compared to normal cells. The full list of altered genes and magnitudes γ j from [35,36] are shown in A1. Details of the model equation and parameters, and the log 2 -fold change is in the range of γ j ∈ [−3.5, 2.7]. In addition, we consider the extreme case of gene alteration as the maximum level log 2 g j i + 1 = 16 for up-regulated genes and log 2 g j i + 1 = 0 for down-regulated genes. Figure 5 shows examples of the gene expression levels in log scale that we modify including the up-regulated genes, GPR56, GATA2, and MZB1, and the down-regulated genes, LGALS3, LY86, and ANXA5. The given single-cell data in normal condition is plotted together with the hypothetically altered levels of gene expressions in regular leukemia pathogenesis and extreme levels of alteration. Although the case of extreme alteration is unrealistic, we consider such case to illustrate an example where the graph abstraction and dimension reduction algorithm clearly distinguishes the leukemic cells.

Effects of gene perturbation on graph abstraction and multi-dimensional reduced component space
We first study the sensitivity of the reduced component space using diffusion mapping [15]. Figure 6A,B compares the altered leukemic single-cell data g i projected on the normal reduced space (θ 1 ,θ 2 ). The left-most column shows the projected leukemic single-cell data P g i in the normal reduced space, where the leukemic cells are located toward the left-bottom compared to the normal cells in Nestorowa data, and upwards in Paul data. The effect of gene modification is shown more clearly in the presented vector field P g i − P g i .
Similarly, we study the impact of leukemia-associated gene perturbation in graph abstraction using PAGA [23]. Figure 6C,D shows the clustered cell types and the corresponding graph using perturbed leukemic scRNA-seq data. The presented results are computed with Nestorowa data. The clustered cell types and leukemic cells are annotated to depict the cluster properties. In Figure 6C, which is the case of leukemia progression with single-cell data altered in regular magnitude, we observe that there is no cluster that separates the leukemic cells. Thus, the information of leukemic gene alteration is lost within the clustering algorithm, and the model on such abstracted graph is not capable of studying the perturbed cells. On the other hand, when the gene levels are modified to their extreme levels, the perturbed leukemic single-cell data are clustered into separate nodes as shown in Figure 6D. In this case, although the graph model is able to study the perturbed cells as separate nodes, we comment that this level of perturbation is an unrealistic scenario due to the extreme levels of gene expression.
A strength of the multi-dimensional cell state model is its capability of interpreting the perturbation of gene expression levels or new incoming cell data regardless of its relation to the primary data ( Figure 6). As shown in our results, the leukemic alteration is successfully projected in the reduced space, while the abstracted graph lost the information. Although the projected directions in the reduced space can be once more projected on the graph, it does not make sense to do so when the direction is orthogonal to the edges. The multidimensional space model has its advantage especially in this case, where the projected direction of cell states can be directly implemented.

Simulating AML pathogenesis by perturbing known leukemia-associated genes
In this section, we incorporate the perturbed leukemia-associated gene data in the AML simulation using the multi-dimensional space model. In particular, we are interested in studying the impact of leukemia-associated gene alteration on the cell distribution during AML progression. We compute the abnormal differentiation of leukemic cells by projecting the altered single-cell data of MEP and Ery cells to the normal diffusion component space P g i as it is shown in Figure 7. The aberrant differentiation vector v aml 1 = P g i − P g i shows the shifts of cells toward the location where no cell data occupies. Therefore, in addition to modifying the advection term according to the altered gene data, we assume an emergence of new abnormal cell state. In particular, we take the cell state location at θ* = (0.610, 0.215) in Nestorowa data, and at θ* = (0.6, 1.0) in Paul data, and use Gaussian functions centered at θ* to obtain v aml 2 . The corresponding vector fields are also shown in  Figure 8B, where the effects of the parameters, c aml and r aml , are shown more clearly. The total number of cells increases more than 10 times the initial size after t = 30 when c aml = 10 and r aml = 0. When the over-proliferation term is appended as r aml = 1, the total number of cells increases more rapidly, for example, up to 100 times the initial size and the number of cells in most of the myeloid lineage increases. Our simulation results agree with the experimental data, where unconventional cell states emerge during leukemia progression and eventually overtakes the entire progenitor population as observed by FACS analysis of bone marrow progenitor cells isolated from CM knock-in preleukemic and leukemic mice ( Figure 8C). The predominant population observed in leukemic bone marrow does not fall within the typical gates in conventional cell clustering based on data from normal control mice ( Figure 4E). Although this novel population would had been classified as MEP, pre-Meg/E, Pre-GM, and GMP cells in the graph model ( Figure 4A,B), we emphasize that they are distinct population and the multi-dimensional model is capable of incorporating novel cell states. Although we comment that, in the graph model, a new cell type can be included by adding a new node to the original graph.

Interpretation of new cell states in the multi-dimensional model
The remaining question is how to interpret the new cell states in the multi-dimensional space model that may arise far away from the cell states identified by conventional markers. Hence, we propose some measures in Eqs (2.12),(2.13) to guide the interpretation. Figure 9 shows an example of the rescaled correlation quantities r f,j and r v,j computed with Nestorowa data. The first row show results of the correlation r v,j to the average leukemic directional

Discussion
We have shown how to construct mathematical models of cell state-transitions using scRNAseq data. We compare two cell state geometries: solving equations on graphs and solving equations on a multi-dimensional cell state-space. Each cell state geometry has its strengths and limitations. Selecting a model for a given application or dataset will depend on the type of biological data and the nature of the scientific question.
When the modeling application and quantity of interest includes well-known cell lineages and relation between the conventional cell states, the graph model is more appropriate due to its ability of distinguishing distinct cell lineages more clearly compared to the multi-dimensional space model. Dynamics of cell numbers in specific cell states, alteration of proliferation and apoptosis in particular cell state, differentiation block, and emergence of intermediate cell states can be quantified and studied in a straightforward manner.
However, to explore cell states beyond known cell lineages, the continuum space model is more advantageous since it includes all intermediate and pathological cell states, rather than confining the model into presumed cell lineages. Moreover, the continuum model can incorporate a relatively small genetic and epigenetic alteration that the graph abstraction may not recognize, and study abnormal trajectories that yield unconventional cell states.
We selected and perturbed genes to simulate AML based on genes known to be associated with leukemia pathogenesis. We do not intend for this to be an accurate model of the biological process, rather, as an illustration of how one may select sets of genes and perturb them in a prescribed fashion in order to study the effect on cell state-transition dynamics. This approach assumes that AML pathogenesis originates from changes in gene expression in specific cell subsets, which is limited by our identification of these genes based on published literature. We acknowledge this is a limitation of the modeling approach, although we also note that our model predictions are consistent with known features of AML progression.

Comparison to other approaches
Although at the time of this work there are relatively few mathematical models published which utilize single-cell sequencing data, there are a few notable exceptions. Of particular note are works which use modeling and simulation to generate synthetic in silico gene expression datasets [37]. These important approaches to mechanism-based mathematical modeling may also be used to study and predict the effects of perturbations on cell state distributions. They may also be used as computational controls to benchmark analysis tools and potentially to benchmark and compare mathematical models, although using a model to benchmark other models can lead to consistent but incorrect circular reasoning and caution is warranted. Another example is Ferrall-Fairbanks and Papalexi et al, who use mathematical analysis to generate novel quantifications of cell heterogeneity in cancer or immune cell subsets respectively [38,39]. These methods may be used to map and interpret novel cell states predicted by mathematical models or similarly as a method to interpret model-predicted changes in cell heterogeneity following a perturbation.
Schiebinger et al compute and predict differentiation trajectories in cell development using optimal transport (OT) [40,41]. This approach considers the optimal transport of cells as a mass flowing along differentiation trajectories, and is conceptually the most similar to our approach. As presented, Schiebinger et al do not use the OT framework to examine perturbations of cell states or genes along the differentiation trajectory, although this is possible with an OT model. Setty et al present a method to compute cell fate probabilities [42], which may also be achieved by inferring cell state-transition dynamics with lineage trees [43]. Fischer et al have demonstrated a method for inferring population dynamics from single-cell sequencing data [44], where the model equation is identical to our graph based model developed in [12]. Jiang et al develops a dynamic inference approach to derive a Fokker-Planck type PDE on a graph considering an energy landscape and optimal transport [45]. Sharma et al use longitudinal sequencing to study drug-induced infidelity in the stem cell hierarchy [46], and Karaayvaz et al show how to use single-cell sequencing to examine drug resistance in breast cancer [47]. These approaches and analysis methods may be used to inform and potentially calibrate mathematical models of cell population dynamics or response to treatment-induced perturbations.
Recently, vector fields derived from RNA velocity [48] have been used to infer potential energy or fitness landscapes for cell state-transitions. These approaches may be used to inform the computational domain for mathematical models that we present here, however, an important limitation of the RNA-velocity approach is extrapolation of the vector field outside of the data range. This underscores the need for hypothesis-based and model-guided approaches to inform the shape of these fields. This limitation also applies to the rapidly growing field of deep learning approaches [49] to analyze single-cell sequencing data, namely, whether the learning algorithm can effectively make predictions to datasets which are not sufficiently similar to those upon which it has been trained. We believe that the future likely involves a merger of mathematical modeling with machine learning, in which mathematical models are used to inform learning approaches and impute sparse data as has been recently shown by Gaw and Rockne et al [50,51]. Among the recent works that align with this direction, PRESCIENT algorithm aims to learn the underlying differentiation landscape from time-series scRNA-seq data [21]. Moreover, dynamo framework improves RNA velocity using kinetic models to reconstruct continuous vector fields that predict cell fates [52].

Opportunities and limitations of modeling with single-cell sequencing data
There are pros and cons, opportunities and limitations to mathematical modeling with single-cell sequencing data. The advantages and potential opportunities include: a wealth of available data, richness and complexity of each data set, a focus on the cell level, opportunity to study dynamics in hierarchically structured state-based relationships between cells, and an ability to perturb individual cells and/or genes within cells to predict dynamics of state-change at cellular level. The most significant strength of mathematical modeling is the ability to use and generate hypotheses that may not be directly evident from the data; for example, extrapolation of RNA velocity fields beyond the dataset boundaries or to interpret and predict novel cell states which may not otherwise be clearly identified with known canonical cell state markers. Another advantage of our approach is the use of pseudo-time analysis of data collected at a single timepoint to calibrate the models, however, the models can also be calibrated directly to time-sequential single-cell datasets, which we expect to become more commonly available as single-cell sequencing continues to be used as a tool to study cell dynamics.
The disadvantages and limitations include: the potential for misleading or incorrect inference due to poor data quality including drop-outs, small non-representative samples of large heterogeneous populations, batch effects, no physical or micro-environmental context, no direct or physical interactions between cells, and the possibility of model predictions to be sensitive to methods of dimension reduction, graph abstraction, state-space construction, and potentially sequencing platform. Sensitivity of the modeling to experimental and computational methods may be directly studied and mitigated as we have shown in this work, however this remains potentially a significant source of uncertainty and variability in the modeling calibration and predictions. Studying the sensitivity of our modeling framework regarding different noise scenarios and applying noise reduction methods is our future work [53].
In terms of computational cost, the graph model is more efficient since it is a multiple of one-dimensional cost, while the cost of implementing the space model increases exponentially as the dimension of reduced space increases. In our simulation, the computational cost to simulate up to time t = 50 with step size Δt = 10 −3 and O(100 2 ) degrees of freedom in one-dimension is around 25 seconds in the graph model with 8 nodes, while it takes around 230 seconds in the continuum model with two dimensions. In short, the continuum model runs approximately 10 times longer than the graph model with 8 nodes in our example. Therefore, the multi-dimensional cell state geometry will be reasonable only when the reduced component can be truncated at two-to three-dimension, unless the numerical method is carefully implemented, and we emphasize that the graph model will be more advantageous in terms of computational cost than the continuum model especially when higher dimensional reduced space is necessary.

Future work and applications
Future applications of this approach is to explore hypothesis in the resolution of single-cell genomics and study altered and novel cell states with genetic and epigenetic alterations in various biological systems and pathogenesis. We look forward to compare the model prediction to sampling/sequencing of perturbed biological system, for instance, to examine scRNA-seq data from leukemic progenitor cells. Moreover, we anticipate to incorporate effects of external perturbations such as therapy in future studies.
There are opportunities for further enhancements in our model in improving the model of cell landscape dynamics to accurately estimate cell transition pathways in the reduced component space, for instance, minimum action paths [6] and bifurcation [7,54]. The model can be improved by obtaining parameter functions or mappings of biological quantities directly from single-cell sequencing data, for example, more precisely infer the proliferation rate function. Also, developing methodologies to obtain reduced component space that captures desired characteristic of cell states [55] will help us explore our approach for other biological settings where cell states are less clearly characterized. Moreover, we propose to develop quantities, such as index of critical state transitions [54,56], in the phenotype space that could be used to predict forthcoming major alterations in development and diseases. We also expect to be able to infer the potential landscape directly from the RNA velocity vector field [48,52].

Conclusions
In summary, despite the explosion of computational tools to analyze single-cell sequencing data, there have been relatively few mathematical models developed which utilize this data. Here we begin to explore the possibilities-and limitations-of dynamical modeling with single-cell RNA-seq data. We hope this work paves the way for development of mathematical models to guide the interpretation of these complicated datasets as they begin to be collected after biological perturbations (eg., cancer, treatment, altered developmental processes), sequentially over time, or sampled spatially within biological tissues.
where we take θ = 0.04. The self-renewal rate functions a k (x) and a(θ) are computed similarly. See Cell proliferation rate r(θ) and self-renewal rate a(θ) computed from the single-cell data. The black dots are the rates of data. for r(θ) and a(θ) computed for Nestorowa data.
To compute the multi-dimensional function on the continuum space from the single-cell data, we employ the kernel density method [18,26], that is a non-parametric way to estimate the density function based on a finite data sample. Using the single-cell samples in the reduced component space, θ i i = 1 N , the method approximates the density function as where K is the kernel smoothing function that we take it as a Gaussian function and h is the bandwidth. The optimal bandwidth to estimate normal densities can be computed by , where σ is the standard deviation and N is the sample size, and the optimal bandwidth for our data is computed as 0.0383 to 0.0456, however in our simulation, we choose a slightly smaller value, h = 0.03, to reveal more features of multiple modes. Figure  2(C,D) shows the results computing u s (θ) using Nestorowa and Paul data, and A2(a) shows the corresponding ν 1 (θ).
In the diffusion term, we explore the parameter ν so that the phenotypic instability does not dominate the cell maturation. We compute the parameter in the range of v ≤ L/T d 2 /4, where the distance in the diffusion space is L = 1 and the time that HSC differentiates to the progenitors is T d = 5 ~ 30 (day), that is, ν < 0.0027 ~ 0.01, and we consider ν = 0.001.
Quantifying the local phenotypic instability in the reduced component space, and justifying this term is our future work.
To compute the reduced component space using dimension reduction approaches, we employ diffusion mapping. See [15,62] for the detail of the algorithm. We take the cosine distance, k(x i , x j ) = 1 − corr(x i , x j ) for the Nestorowa data and the gaussian distance For the pseudotime inference, we use the algorithm developed in [63]. The diffusion distance between two cells are computed as  Cell proliferation rate r(θ) and self-renewal rate a(θ) computed from the single-cell data.
The black dots are the rates of data. Figure A2.
Pseudotime dynamics. The homeostasis cell differentiation vector v 1 (a), and the direction of active cell differentiation obtained from diffusion pseudotime analysis (b), and that interpolated at the grid points v 2 (c) are presented. We remark that, v 2 corresponds to the cell differentiation along the edges in the graph model.      Step-by-step illustration of our modeling process. 1. Processed single-cell sequencing expression matrices are represented in a reduced dimension space through one of many dimension reduction techniques such as diffusion mapping, PCA, t-SNE, or UMAP. 2. Cell clusters are inferred to construct the cell state geometry either the graph or multidimensional continuum of cell states. 3. From these representations, models are calibrated to the transport of cell distribution along the graph or in the cell state space. 4. The models can then be utilized to perturb genes and cell states. The calibrated models predict cell state-transitions and the emergence of novel cell states.   show an increase within ten days in both graph and multi-dimensional space models. E) Perturbing genes associated with leukemia stem cells. Examples of expression levels of genes in log 2 scale, that are associated with leukemia stem cells and pathogenesis, including up-regulated GPR56, GATA2, and MZB1, and down-regulated LGALS3, LY86, and ANXA5. We show these subset of genes simply to illustrate the process. The normal single-cell data log 2 g j i + 1 (blue circle) and modified gene expression log 2 g j i + 1 (red square) computed as Eq (3.1) are shown together, with the case of extreme levels of either 16 or 0 (purple diamond).    Table 1.
Comparison of the cell state model on graph versus multi-dimensional space in n dimensions. The computational cost is estimated by denoting M as the number of discretized grid points in one-dimension. We comment that the computational cost of a PDE solver can be up to a third power of the degree of freedom.