We employed two different approaches to model the cell differentiation processes. In the first approach we used the experimental data to visualize the Waddington landscape of early mouse embryogenesis and identified the clusters of the genes differentially expressed in each developmental bifurcation. In the second approach, we theoretically compared the dynamical systems generated by the smaller (two-TF) and the extended (two-cluster) regulatory circuitries, using ODE based models.
Waddington landscape: preprocessing of the experimental data
We obtained the expression profiles of 48 genes in 442 single mouse embryonic cells from zygote to 64 cells stage, that were generated by the TaqMan qRT-PCR assay [8]. These genes were selected after analyzing the expression levels of 802 TFs, due to differential expression in blastomeres or known function in early development. The initial Ct values ranged from 10 to 28, and the expression values were assigned by subtracting the Ct values from the baseline value of 28 (see the Additional file 3). PCA was performed using the mean-subtracted expression values. Correlation heatmap of the genes was generated based on pairwise Spearman correlations of the expression profiles of the cells in the 64-cell stage.
Axes of Waddington landscape
In order to visualize the Waddington landscape of the preimplantation development, we needed to define each dimension and compute it. There are three axes (dimensions) in the epigenetic landscape, as illustrated in Fig. 1a: (i) The x-axis (left-right) through which distinct cell fates are shown as different attractors (valleys). (i) The y-axis (back-front) that shows time of development, as early and late developmental stages are located in backward forward of the landscape, respectively. (iii) The z-axis (down-up) that represents a potential function, which integrates both differentiation potency and stability [39]. The totipotent cells (zygotes) are posed at the highest valley. As the cells undergo more differentiation into pluripotent, multipotent and then unipotent cells, they go towards the deeper and lower valleys. Furthermore the stable cell states (attractors) are distinguished as valleys from the instable and transient cell states that form hills.
It was straightforward to assign the y-axis of the cells since the time of development was available for each cell. To establish the x-axis of the epigenetic landscape, we computed the principal components PC1 and PC2 of the gene expression profiles (Fig. 1b). The coordinates origin was slightly moved into a cell-free region (PC1 = −0.5, PC2 = 0) to ensure all the cells of the same fate are located in the same side of the origin and have close angular coordinates. Then the x-dimension of each cell was computed as its angular coordinate around the origin. Through this dimension reduction – from the initial gene expression profiles consisting of 48 dimensions into a single axis – we aimed to preserve the similarities and differences of the cells.
Pseudo-potential function of the Waddington landscape
We needed to define a form of potential function from the experimental data. The closed form of a potential function is restricted to the gradient systems with stringent mathematical conditions that usually do not hold in biological systems [40]. As a result, most of the previous studies have defined pseudo- or quasi-potential functions based on many different methods: the ODEs with path integration [15, 40], Fokker-Planck equation [41, 42], Langevin dynamics [43], Hamilton-Jacobi equation [44], drift-diffusion models [45], Boltzmann distribution [46] and stochastic simulation [47]. Signaling network entropy, as a measure of promiscuity or undetermined lineage, is the other framework used to define a pseudo-potential function based on the experimental data [48].
In this study we employed the Boltzmann (Gibbs) distribution, which models the probability distribution of the particles in a system over various states with different energy levels [49]. It makes a connection between the energy levels and the probabilities of the particles being in each state. The Boltzmann distribution is expressed as the following equation:
$$ \frac{P(A)}{P(B)}={e}^{\raisebox{1ex}{$-\Delta E$}\!\left/ \!\raisebox{-1ex}{${k}_BT$}\right.} $$
where A and B are two different states, P(x) is the probability of a particle to be in state x, ΔE is the energy difference that a particle should absorb/release to change its state from A to B, k
B
is the Boltzmann constant, and T is the system temperature. By taking the logarithm of two sides we have:
$$ \ln \frac{P(A)}{P(B)}=\frac{-\Delta E}{k_BT}\ \to\ E(A)-E(B)=-{k}_BT\left( \ln \left(P(A)\right)- \ln \left(P(B)\right)\right) $$
in which E(x) is the energy of a particle in state x. Taking the state B as the pseudo-potential reference results in:
$$ U\left(\mathrm{A}\right)=-\uprho\ \ln \left(P(A)\right)+\omega $$
where U is the pseudo-potential function. Both ρ and ω are constant values that scale the landscape and can be omitted in visualization. To compute the pseudo-potential function we should determine the probability of the cells to be in each state, as follows.
Probability distribution of the cell states
At each developmental stage we assumed the expression profiles of the cells of the same type were normally distributed along the x-axis after the dimension reduction. To check this assumption we produced the Q-Q plots of each developmental stage for the angular coordinates of the cells in the PC1-PC2 plane (see the Additional file 4). Up to the 16-cell stage the points are almost fitting a single trend line. In the 32-cell stage there are two distinguished segments, discriminating ICM and TE cells. Each of three segments in 64-cell stage fit a different trend line, which shows this stage is a mixture of three normal distributions, representing EPI, PE and TE lineages. Furthermore we performed the Shapiro-Wilk normality test [50], that confirmed the normality of several segments of different stages.
As a result we considered a cell population including m different cell types would have a mixture of m normal distributions. By assuming τ
k
as the probability of a cell belonging to the k -th cell type (\( 1\le k\le m,\ {\tau}_k\ge 0,\ {\displaystyle \sum_{k=1}^m}{\tau}_k=1 \)), the mixed probability distribution function is:
$$ f(x)={\displaystyle \sum_{k=1}^m}{\tau}_k{\Phi}_k\left(x\Big|{\mu}_k,{\Sigma}_k\right) $$
where μ
k
and Σ
k
are the mean and covariance matrix of all the cells of the k -th cell type, and Φ
k
is a Gaussian function defined as:
$$ {\Phi}_k\left({\boldsymbol{x}}_{\boldsymbol{i}}\Big|{\mu}_k,{\Sigma}_k\right)=\frac{1}{\sqrt{2\pi \left|{\Sigma}_k\right|}}{e}^{-\frac{1}{2}{\left({\boldsymbol{x}}_i-{\mu}_k\right)}^T{\Sigma}_k^{-1}\left({\boldsymbol{x}}_i-{\mu}_k\right)} $$
From the above equations we could calculate the pseudo-potential function:
$$ U(x) \propto - \ln \left(f(x)\right) $$
where x is any point on the x-axis of the epigenetic landscape (projection of the gene expression profiles) at some particular developmental time. The mixed distribution and the pseudo-potential function were recalculated for each developmental stage with the available experimental data. A linear interpolation was used to fill the gaps between consecutive developmental stages. The landscape was visually tilted to show the reduced differentiation potency during development.
Selecting the number of different cell types and assigning each cell to one of cell types can be done either manually (supervised) or computationally (unsupervised). To have an objective and automated framework, we used the unsupervised approach, using the “mclust” package [51] of R statistical language. The projected expression profiles were given to the package to compute the probabilistic model parameters, including the number of cell types (clusters), the mean and covariance values, based on a maximum likelihood criterion. For additional details, one may refer to the “mclust” package reference manual.
Dynamical modeling: Phase space representation of the regulatory circuitry of two transcription factors (TFs)
For the dynamical system analysis of the two-TF regulatory circuitry, we employed the following set of ODEs [3, 39]:
$$ \frac{du}{dt}=\alpha \frac{u^n}{s^n+{u}^n}+\beta \frac{s^n}{s^n+{v}^n}-\gamma u $$
$$ \frac{dv}{dt}=\alpha \frac{v^n}{s^n+{v}^n}+\beta \frac{s^n}{s^n+{u}^n}-\gamma v $$
where u and v are concentrations of the pair of opposite TFs, and α and β are the strengths of the positive and negative regulatory interactions, respectively. For simplicity we used the same protein degradation rate γ for both TFs. The term \( \frac{x^n}{s^n+{x}^n}\ \left(x\in \left\{u,v\right\}\right) \) is a sigmoid function that has 0 value at x = 0, increases to 0.5 at x = s, and asymptotically approaches 1 at the large values of x. It resembles the positive auto-activation regulatory effect of each TF. The steepness of the sigmoid function is defined by the power n. On the other hand, \( \frac{s^n}{s^n+{x}^n} \) is a decreasing sigmoid function that starts from 1 at x = 0 and approaches 0 at large x values, which resemble the mutual inhibitory effects.
In order to model the perturbation in the form of increased decay rate of a particular protein, we increased the degradation rate of the TF u by 50 %, as denoted by γ*:
$$ \frac{d{u}^{*}}{dt}=\alpha \frac{u^n}{s^n+{u}^n}+\beta \frac{s^n}{s^n+{v}^n}-{\gamma}^{*}u $$
Phase space representation of the regulatory circuitry of two clusters of transcription factors (TFs)
To analyze the dynamical system of a gene regulatory circuitry consisting in two clusters of TFs we generalized the previous two-TF model by using the following equations:
$$ \frac{dx}{dt}=\frac{d{u}_1}{dt}+\frac{d{u}_2}{dt}=\eta \left({u}_1,{u}_2,{v}_1,{v}_2,\gamma \right)+\eta \left({u}_2,{u}_1,{v}_2,{v}_1,\gamma \right) $$
$$ \frac{dy}{dt}=\frac{d{v}_1}{dt}+\frac{d{v}_2}{dt}=\eta \left({v}_1,{v}_2,{u}_1,{u}_2,\gamma \right)+\eta \left({v}_2,{v}_1,{u}_2,{u}_1,\gamma \right) $$
where u
1 and u
2 are the concentrations of two proteins of the first cluster, v
1 and v
2 denote the second cluster protein concentrations, and x = u
1 + u
2 and y = v
1 + v
2 are the total concentrations of the proteins in clusters 1 and 2, respectively. We defined the generic function η(a, b, c, d, γ) to compute the concentration rate of any protein a based on the concentration values of the TFs a and b in one cluster, and c and d in the other cluster, as follows:
$$ \eta \left(a,b,c,d,\gamma \right)=\alpha \frac{a^n+{b}^n}{s^n+{a}^n+{b}^n}+\beta \frac{s^n}{s^n+{c}^n+{d}^n}-\gamma a $$
For the perturbation analysis, we used the increased degradation rate γ* for the TF u
2:
$$ \frac{dx}{dt}=\frac{d{u}_1}{dt}+\frac{d{u}_2^{*}}{dt}=\eta \left({u}_1,{u}_2,{v}_1,{v}_2,\gamma \right)+\eta \left({u}_2,{u}_1,{v}_2,{v}_1,{\gamma}^{*}\right) $$
In both models we used the following parameters: n = 4, s = 0.5, α = 1.5, β = 1, γ = 1 and γ* = 1.5. The sample space of (u, v) = [0, 3]2 was used for analyzing the two-TF model, and (u, v, u
1, v
1) = [0, 3]4 for the two-cluster model.
Simulation of the cell differentiation in absence or presence of perturbation
For each of the four different regulatory circuitries depicted in Fig. 3c, d and Fig. 4a, b we simulated the differentiation of 1000 cells. The initial expression rates of the TFs in each cell were assigned from a normal distribution with μ = 1.5 and sd = 1. With this selection of parameters, the majority of the cells were initially in proximity of the progenitor state (the attractor 2 of Figs. 3 and 4).
Each simulation continued 100 steps, in which the expression rates of the TFs in each cell were slightly changed, based on two factors: the dynamical system forces (differential equations above) and a standard Gaussian noise (μ = 0, sd = 1). The strengths of these factors were tuned by two coefficients: the force field coefficient had a constant value of 0.2 during the simulation; and the noise coefficient that started with 0.5 and gradually reduced during the simulation (multiplied by 0.98 in each step) to ensure the convergence of the experiment.
The auto-activation strength α was 1.5 at the beginning of each simulation, but gradually reduced (multiplied by 0.98 in each step). In this way, we forced the cells to leave the progenitor state 2 and differentiate into the attractor states 1 or 3. During this process, the stability of the attractor 2 was gradually decreased and resulted in a bistable system with only attractors 1 and 3. In each attractor of the bistable system, one TF was silent and the other was expressed at a slightly lower rate than the initial circuit configuration, due to the lower value of α.
Implementation
The code was implemented in R statistical language [52]. We used the packages “mclust” to generate the mixed Gaussian model, “rgl” for 3D visualization of the Waddington landscape, “ggplot2” for 2D visualization of the data [53], and “pheatmap” for visualization of the correlation heatmap. We also used the packages “grid”, “gplots”, “plyr”, “Hmisc”, and “Biobase”.
Advantages and limitations
Our method of visualizing the Waddington landscape enables the application of the experimental data at single cell resolution for this purpose. While we used the gene expression profiles of early embryonic cells, our method can be generalized for analysis of the high-throughput DNA methylation, histone modifications and non-coding RNA expression profiles. It is computationally fast and can be used for whole-genome scale of data and a large number of single cells. By application of time-course data, the same method can be applied for visualizing the landscape of reprogramming, transdifferentiation or stem cell differentiation.
Our method interpolates the developmental time between each pair of successive sampling time points; hence the closer the sampling time points, the more realistic the resulting landscape. The valley depth in this method mainly represents the number of cells assigned to the corresponding attractor state. This requires the data to be generated by random sampling of different cell types. For study of distant cell types, the quantity of cells and the depth of attractors can be influenced by cell division rates. In this case we suggest combining our method with an indicator of differentiation potency or stability, such as the cellular network entropy [48].
Availability of supporting data
The preprocessed single-cell resolution gene expression profiles of mouse preimplantation embryonic cells [8] are provided in the Additional file 3. We have also provided in the same additional file the complete source code of this study in R programming language.
