Modeling cell differentiation by continuous birth-death process

Cell differentiation is clearly a non-equilibrium process due to frequent birth and death of cells. Thus, to model the cell differentiation, we use the continuous birth-death process as the underlying dynamics of cell differentiation, which is also named as the source-sink Fokker-Planck equation in mathematics and the population balance equation in Klein at al.’s work [27, 29, 30]. The continuous birth-death process assumes that the probability density function c(x, t) of all sample cells develops as (1) where x is a vector of gene expression, t is the time, F(x) is a potential function, D is the noise amplitude, and R(x) is the net-flow of cells at state x. ∇, ∇⋅, and Δ denote the gradient, divergence, and Laplace operators, respectively. When the system reaches a non-equilibrium steady state, i.e. lim_t→∞ c(x, t) = p(x) or ∂c(x, t)/∂t = 0, the potential can be decomposed as F(x) = U(x) + V(x) and calculated by (2) (3) according to [27]. We name the flow that is differentiation-oriented as “advection”, while we use the term “diffusion” which is caused by random noise [31–33]. U(x) is known as the equilibrium potential caused by diffusion without birth and death, and V(x) is a new potential caused only by advection without diffusion. By definition, noise is known to only influence the diffusion process and generates meta-stable wells in U(x). Therefore, V(x) can be taken as the cell differentiation potential to describe the differentiation direction. The cells will evolve from pluripotent cells with high V(x) to differentiated cells with low V(x). V(x) can represent the Waddington landscape [28] of the cell at state x, and the additive inverse of V(x) can be considered as a reflection of the pseudo-time (see Materials and methods).

Estimating cell potential V(x) and constructing differentiation paths

To obtain V(x) from Eq (3), the net-flow R(x) and backward operator need to be estimated from scRNA-seq data. Due to the limitation of sample size, it is difficult to accurately measure net-flow R(x) for every cell. In contrast to [27], which required additional information, we first clustered samples into different cell types/clusters, and then computed the net-flow for each cluster s from the gene expression matrix by the divergence theorem and marginal decomposition. On the other hand, benefiting from the diffusion map theory [34, 35] and model reduction [36], the backward operator was approximated by , which was the coarse-grained discrete matrix representation of between cell clusters, obtained from the cell-to-cell transition matrix. Thus, with the approximated net-flow for every cluster/cell type, and the approximated transition operator , could be obtained numerically, generating a concrete value for the potency of each cell type. Additional details can be found in Materials and methods. We remark that in [27], the net-flow rate was set as the prior knowledge for each cell. However, we were able to obtain the value from the gene expression matrix if cells were clustered into different metastable states. This is one key point of our work.

To illustrate the entire differentiation process, we constructed its landscape, in which nodes were cell types with potential , paths were determined by the transition matrix between clusters, and directions were from high to low potential. We named this procedure as Landscape of Differentiation Dynamics (LDD). Fig 1 provides a flowchart of LDD, which is described in details in Materials and methods. An algorithmic description is provided in S1 Text section S1 and Fig I in S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Flowchart of the Landscape of Differentiation Dynamics (LDD) method.

A: A pool of single cells, from which we can obtain the gene expression matrix by single cell sequencing. B: After preprocessing, quality control, and dimension reduction, a low-dimensional data matrix is obtained. C: The samples are clustered into different types. Undirected differentiation paths are determined by a transition matrix between clusters. D: After applying the continuous birth-death process to model the whole differentiation process, the potential V(x), differentiation directions, and landscape can be constructed.

https://doi.org/10.1371/journal.pcbi.1007488.g001

Computing differentiation landscape of simulated models

We used three simulated models to construct potential and verify our LDD method. The first is a drift-diffusion process, in which particles evolve by (4) where x is the position of a particle in 50 dimensions, F is a potential function, D is the noise amplitude, and ξ is a normal random vector standing for noise. A bifurcation from one branch to two branches occurred in the system, representing cell differentiation. Particles/cells were generated from a source, and two sinks indicated places for particle/cell removal or death. It imitated the cell’s lifespan from birth to death. Details about the model are shown in S1 Text. In the example, 400 samples/cells were simulated. After reducing to two dimensions using principal component analysis (PCA) and clustering the 400 samples into four groups/clusters by k-means, we computed the potential of each cluster. The landscape in a three-dimensional view is illustrated in Fig 2A. Fig A(a) in S1 Text shows the potential in a two-dimensional space. In Fig 2B, the differentiation paths between the four clusters were constructed. The cluster with the highest potential corresponded to the source region, and the two sinks with low potential values were at the end of two lineages.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Differentiation landscape, differentiation paths, and gene networks for simulated models.

A and B are the LDD potential landscape and differentiation paths using data from the simulated drift-diffusion process, in which samples/cells were clustered into four groups. C and D are the two-gene and six-gene regulatory networks for simulation, respectively. For the two-gene network, the potential landscape and differentiation paths for clusters are shown in E and F, respectively, where seven clusters were detected. Constant potential V for each cluster was computed by LDD, while the landscape is an illustration constructed by the method in Materials and methods.

https://doi.org/10.1371/journal.pcbi.1007488.g002

The next two simulated examples are two-gene and six-gene regulatory networks. Their gene interactions are shown in Fig 2C and 2D, respectively. A region for cell birth is shown, as well as several regions for cell death. The detailed simulation method can be found in S1 Text. For the two-gene network, as the genes inhibit each other, two lineages formed. In each branch, one gene had high expression, while expression of the other gene was very low. The illustrative landscape for the two-gene network is in Fig 2E, and the two-dimensional plot is in Fig A(b) in S1 Text. Fig 2F shows the differentiation paths between seven groups. LDD generated the potential in a correct order, i.e. pluripotent cells with higher values, and differentiated cells with lower values. Similar results of LDD were obtained in the six-gene network, which had four branches for the samples. Only two genes had high expression, while the others were low in each branch, which are shown in Fig A(c) and A(d) in S1 Text.

The three examples demonstrate the ability of LDD to quantitatively characterize the stochastic process of cell differentiation from a dynamic perspective, due to our model based on the continuous birth-death process. In particular, only using the observed data, LDD could quantitatively show that pluripotent cells with higher potential differentiated into downstream cells with lower potential, similar to a ball rolling down a mountain as described by Waddington [28].

Computing differentiation landscape of real datasets

To further verify the efficiency of LDD, we applied it to four real datasets, i.e. Guo’s dataset, Nef’s dataset, Xu’s dataset, and Furlan’s dataset.

Guo’s dataset [37] describes cells developing from zygote to blastocyst, through oocyte, 2-cell, 4-cell, 8-cell, morula, E3.5 blastocyst, and E4.25 blastocyst stages. At the end term of morula, the cells could differentiate into trophectoderm (TE) and inner cell mass (ICM). Nef’s dataset [38] focused on the determination of mouse’s sex. 400 samples were selected from E10.5, E11.5, E12.5, E13.5, and E16.5 stages. Two branches named the interstitial progenitor cell lineage and Sertoli cell lineage appeared during the observation time. Xu’s dataset [39] showed that the progenitor hepatoblasts had bipotency to divide into hepatocytes and cholangiocytes. Furlan’s dataset [40] found that a large number of chromaffin cells in adrenal medulla arose from a kind of peripheral glial stem cell called Schwann cell precursors (SCPs). The description and the preprocessing of these datasets can be found in S1 Text.

Using LDD, we computed the potential and differentiation paths of each dataset. The illustrative landscapes, potential values, and lineages of Xu’s dataset are displayed in Fig 3A–3C, which show that the progenitor hepatoblasts (cluster 1) with high potency differentiated into hepatocytes (cluster 3) and cholangiocytes (cluster 5) with low potency. Fig 3D–3F are the corresponding results of Furlan’s dataset. There was a single path through which the SCPs became chromaffin cells. For Guo’s dataset and Nef’s dataset, the results are given in Fig B in S1 Text. For Nef’s dataset, we also computed Pearson correlations between LDD potential and three character genes Pbx1, Gpc3, and Sfrp1 as 0.9015, 0.8582, and 0.8000. These genes decreased during the differentiation process, which indicated that the pseudo-time conformed to the differentiation direction.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Differentiation landscape and differentiation paths for real datasets.

A, D are the LDD potential landscapes, B, E are the potential values plotted in the two-dimensional reduction space, and C, F are the differentiation paths. A-C use Xu’s dataset, which describes that hepatoblasts (cluster 1) differentiate into hepatocytes (cluster 3) and cholangiocytes (cluster 5). D-F use Furlan’s dataset and show that chromaffin cells (cluster 4) are generated from SCPs (cluster 1).

https://doi.org/10.1371/journal.pcbi.1007488.g003

Comparison with other pseudo-time algorithms

From both the simulated and real datasets, we showed that LDD could obtain the correct differentiation paths by quantifying potential values of cells. The pseudo-time is another important topic in single cell research, which provides a time label to each cell. We set the additive inverse of potential as our pseudo-time by LDD, and compare it with several traditional methods. Because entropy-based methods usually need additional information, such as gene interactions or function annotations, we only compared LDD with six distance-based methods, i.e. TSCAN [41], Monocle2 [15], Diffusion Map [11, 12], DPT [13], SLICER [42], and Slingshot [43]. Their properties are listed in Table 1, and the details are given in S1 Text.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Properties of different pseudo-time methods.

https://doi.org/10.1371/journal.pcbi.1007488.t001

For the simulated datasets, we chose the Pearson correlation between the pseudo-times and the true-time labels as the measurement. The second to the fourth columns in Table 2 show the results for the Simu1 (the drift-diffusion process), Simu2 (the two-gene regulatory network), and Simu3 (the six-gene regulatory network) datasets. The closer the Pearson correlation approximated to 1, the more reliable the algorithm was. LDD and DPT outperformed the others; however, DPT required additional information, i.e. a started root sample to represent the stem cell. Hence, it implies that LDD effectively uses the information in the samples, which leads to correct differentiation paths.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Comparison between seven pseudo-time methods.

https://doi.org/10.1371/journal.pcbi.1007488.t002

For the real datasets, the accurate cell time is unavailable; however, several stage-level time labels, such as E10.5 and E11.5, are given. Under one stage label, there are several cell types. Therefore, instead of the Pearson correlation, we used the one-sided Wilcoxon ranksum test to determine whether the stem cell stage had earlier pseudo-time than the differentiated stage. The p-value was chosen as the measurement [19, 20]. The fifth to the eighth columns in Table 2 list the p-values for the four real datasets (Guo’s, Nef’s, Xu’s, and Furlan’s). The alternative hypothesis was that pluripotent stem cells had earlier pseudo-time than differentiated cells. Thus, the smaller the p-value was, the better the algorithm performed. From Table 2, LDD was always shown to produce good results among the seven methods (the second place in Nef’s dataset and the best in other three datasets).

The conclusion from the comparison between different pseudo-time methods is as follows: LDD and DPT performed well in determining the pseudo-time but DPT required additional information, i.e. a root cell. Thus, LDD is an efficient approach to model cell differentiation and quantitatively characterize the Waddington potential landscape or differentiation paths by exploiting dynamical and stochastic features of the differentiation process from the measured data.

Cell differentiation dynamics as continuous birth-death process

We use a stochastic non-equilibrium process Eq (1) to model the cell differentiation, which is also named the source-sink Fokker-Planck equation or the population balance equation [27, 29, 30]. In Eq (1), the function c(x, t) represents the probability density function (pdf) of cells at x, which will be estimated from scRNA-seq data. The term R(x) is crucial to understand the birth and death of the cells involved in the system. If there is a source and cells are increasing locally, R(x) will be positive. Conversely, if there is a sink and cells are removed or die, R(x) will be negative.

Deriving cell differentiation potential V(x)

Assume that the system reaches a non-equilibrium steady state, which means cells keep being born and dying but the whole population distribution is invariant. Let t → ∞, then we will obtain from Eq (1) (5) where p(x) = lim_t→∞ c(x, t). Only depending on p(x), F(x) could have many solutions, one of which can be written as (6) (7) and (8) The explicit forms for Eqs (6) and (7) are (9) (10) where is the backward Kolmogorov operator (11) U(x) can be considered as the equilibrium potential caused by diffusion without samples’ birth and death, while V(x) can be taken as the potential caused by a birth-death flow without diffusion. If the noise amplitude D approaches zero, the diffusion vanishes and thus we can take V(x) as the potential to describe cell pluripotency, or as the reverse pseudo-time of cell differentiation.

To compute V(x), instead of setting R(x) as the given values at each point x like [27], we cluster samples into groups/types to obtain for each group only from the expression matrix. The backward operator is also approximated by a coarse-grained defined on the cell clusters. The cell differentiation potential is then computed for each cluster based on a discrete version of Eq (10). The following subsections will show both theoretical and numerical details, and the overall algorithm is described in S1 Text.

Constructing Kolmogorov operator by diffusion map from cell samples

To approximate the backward Kolmogorov operator in Eq (11), we utilize the diffusion map theory [34, 35]. Denote the kernel function by (12) where x and y are two samples in m-dimensional space, and ε is a parameter adjusting the kernel width. If there are N samples/cells {x₁, x₂, …, x_N} in total obtained from probability density r(x), we can define (13) (14) (15) and the transition matrix between samples i and j as (16) The discrete backward Kolmogorov operator is constructed as (17) When the sample size N → ∞ and the kernel width ε goes to 0, the discrete operator L_ε,α tends to be an operator in the continuous space. When α = 1/2, is exactly the backward operator in Eq (11) (see S1 Text for details). In our LDD algorithm, we first construct a weighted undirected k-nearest-neighbor (kNN) network by the similarity measure Eq (14) with the samples as nodes. After symmetrization, we obtain a Markov chain with transition matrix P = (p_ij)_N×N through Eq (16) as (18) The stationary distribution of the Markov chain is defined as μ = (μ_i)_1×N, where (19) which satisfies μ = μ P. P and μ are used in the procedures below. Usually for the Gaussian kernel, the bandwidth is set as the median value of the distances between all samples.

Clustering cells and computing the coarse-grained operator between clusters

There are lots of well-known clustering methods, such as k-means, spectral cluster, dynamical reduction [36], and new packages designed for single cell datasets, such as SC3 [44] and Seurat [45]. It is still a challenge to make a choice among the diverse approaches [46]. However, for our model, we have several limits when choosing a suitable clustering method. One important point is that different cell types should be separated into different groups. Each cell type holds a metastable state in the dynamical system. The other important point is that we need to maintain continuity of differentiation paths. Some nonlinear clustering approaches may fail, as they may scatter different clusters far away or eliminate the transition point that connects different branches. A suitable clustering is one fully using information from the data and coincident with the biological background. In our simulated datasets and real datasets, as the model is built on the original euclidean space, we used the linear principal component analysis (PCA) to reduce dimension and k-means to cluster the data, which has orthogonal invariance and keeps both the separation and continuity perfectly. When choosing the number of clusters in k-means, the following constrains should be considered. (1) From the expression of Eq (S7) in the S1 Text, we can conclude that the best number of the clusters is exactly the number of cell subtypes (the metastable states). The metastable wells are separated by the ridges of system potential F(x). Hence, there is an upper limit for the cluster number, which is decided by the number of metastable wells. (2) For the differentiation process with one lineage, the lower limit of cluster number is two, and for two lineages, it is four; one is for stem cells, one is around the bifurcation point, and the other two clusters stand for two branches. (3) If two neighboring clusters, which do not include the cluster near the bifurcation point, merge into one, Eq (S7) will still hold, and the pseudo-time order of the cells will not change. Under these constrains, we can ensure the LDD robustness for different cluster numbers. The results are shown in Fig J in S1 Text.

If we get K clusters of all cells by one of the above methods, we can define a coarse-grained matrix between clusters s and t as (20) where and n_s is the number of samples in Ω_s (see S1 Text for details). Correspondingly, we can get the approximated operator through by Eq (17) as (21)

Theoretical derivation of net-flow for each cluster by divergence theorem

One of the most important advantages in the model is that we can compute the net-flow rate of each cluster only based on scRNA-seq data. From Eq (1), if we only focus on the samples/cells in one cluster (or cell type) Ω_s, we can define the conditional probability density function (cpdf) as (22) where is the indicator function of Ω_s, i.e. when x is in Ω_s and 0 otherwise. In the long time limit, c(x, t) and r_s(x, t) will converge to the steady distribution p(x) and r_s(x), respectively. By applying the divergence theorem to the equation satisfied by r_s(x) and eliminating equal terms, we can derive an equation satisfied by the net-flow rate of cluster s as (23) The derivation details are shown in S1 Text.

Numerical computation of net-flow for each cluster

The net-flow rate formula can be simplified by marginal density functions as (24) where x = (x⁽¹⁾, x⁽²⁾, …, x^(m))^T ∈ Ω_s, is the marginal density of r_s(x) on the x^(j)-axis, and is the interval that lies in. By approximating in one-dimensional space through the kernel method and summation of the boundary derivatives, we can compute conveniently.

When computing for each cluster separately by Eq (24), the steady state Eq (5) of the original system may break down due to finite sample size effect and the numerical error in discrete computation. A post-processing of is needed to ensure that the system is steady in discrete setting. By integrating Eq (1) in the whole space, letting t → ∞ and using Eqs (22) and (23), we obtain the equation for as follows: (25) If there are n_s samples in cluster s and , the constraint Eq (25) can be written as (26) Then we post-process the results obtained from Eq (24) as (27) in order to satisfy Eq (26). In other parts of this paper, for notational convenience, we still use instead of , but note that it is a post-processed value.

Computing potential and constructing differentiation paths of cells

By Eqs (10) and (21), the potential in the cluster level satisfies , i.e. (28) where is the coarse-grained matrix between clusters, vector is the LDD potential of different clusters, and vector represents the net-flow of each cluster. As the matrix on the left hand side of Eq (28) is degenerate, we need to compute its pseudo-inverse and get the least square solution of as (29)

On the other hand, the structure of differentiation branches can be inferred from the weight matrix between clusters, whose element is given by (30) and will become the transition matrix between clusters after row normalization [36]. Nonzero elements in represent differentiation paths, while in some cases we use a threshold to eliminate small values. The direction of differentiation is determined by from high to low potential. Thus, using and , we can construct the whole differentiation landscape or picture only by the gene expression matrix obtained from an scRNA-seq dataset. Note that the noise D is not required for evaluating the potential landscape V of cells.

Drawing the illustrative landscape

For most of the datasets in this paper, we constructed an illustrative 3D landscape. It is plotted by fitting functions (31) where (x_i, y_i) are the position of N samples in the two-dimensional reduced space, equals to if sample i belongs to cluster s, g(x, y) is a function (usually linear) positively correlated with , and a, b, σ are adjustable parameters. Some high values could be set as NaN when plotting.