Decoding the mechanisms underlying cell-fate decision-making during stem cell differentiation by random circuit perturbation

Abstract

Stem cells can precisely and robustly undergo cellular differentiation and lineage commitment, referred to as stemness. However, how the gene network underlying stemness regulation reliably specifies cell fates is not well understood. To address this question, we applied a recently developed computational method, random circuit perturbation (RACIPE), to a nine-component gene regulatory network (GRN) governing stemness, from which we identified robust gene states. Among them, four out of the five most probable gene states exhibit gene expression patterns observed in single mouse embryonic cells at 32-cell and 64-cell stages. These gene states can be robustly predicted by the stemness GRN but not by randomized versions of the stemness GRN. Strikingly, we found a hierarchical structure of the GRN with the Oct4/Cdx2 motif functioning as the first decision-making module followed by Gata6/Nanog. We propose that stem cell populations, instead of being viewed as all having a specific cellular state, can be regarded as a heterogeneous mixture including cells in various states. Upon perturbations by external signals, stem cells lose the capacity to access certain cellular states, thereby becoming differentiated. The new gene states and key parameters regulating transitions among gene states proposed by RACIPE can be used to guide experimental strategies to better understand differentiation and design reprogramming. The findings demonstrate that the functions of the stemness GRN is mainly determined by its well-evolved network topology rather than by detailed kinetic parameters.

Keywords: gene regulatory circuit; hierarchical structure; random circuit perturbation; stem cell; systems biology.

Figure 1.

Schematic illustration of random circuit perturbation (RACIPE). The gene regulatory network governing a specific cellular function can be divided into two parts—a core decision-making module and the rest functioning as input signals to the core. Through randomization, RACIPE generates an ensemble of mathematical models, each of which is simulated by the same set of chemical rate equations but with randomly sampled parameters. The simulation results of the model ensemble are subject to statistical analysis, such as hierarchical clustering analysis (HCA), and in silico gene/link perturbation analysis.

Figure 2.

The RACIPE method uncovers robust gene states allowed by the stemness GRN. (a) Diagram of the core gene regulatory network governing stem cell differentiation. Red arrows represent excitatory regulation; blue bar-headed arrows represent inhibitory regulation. (b) Probability distribution of the number of stable steady states generated by 10 000 RACIPE models. Different colours represent different cases characterized by different numbers of initial conditions (blue: 1000 times, red: 1500 times and green: 2000 times) that are used to simulate each RACIPE model. Each case was repeated 10 times to estimate the mean and the standard deviation of the distribution. (c) 2D probability density map of the RACIPE-predicted gene expression profiles projected onto the 1st and 2nd principal component (PC1 and PC2) axes. (d) Contribution of each gene to PC1 and PC2. The PCs were obtained by performing the principal component analysis (PCA) using the gene expression profiles from all 10 000 RACIPE models.

Figure 3.

Comparison of the RACIPE-generated gene expression profiles and single-cell gene expression data of mouse embryo. (a) Robust clusters (gene states, coloured hierarchical trees) were identified for both datasets by unsupervised HCA. Four RACIPE-predicted gene states match those from the late stage single-cell gene expression data. The histogram of the predicted expression levels for each gene is shown at the bottom (blue, 50 bins in each histogram). In both heat maps, each column represents a gene; each row represents the gene expression profile of a stable steady state of a RACIPE model (left) or that for a single cell (right). (b) A total of 54 gene clusters (only show 30 here) were identified by HCA. With a minimum probability cut-off of 0.005, we identified 15 clusters, referred to as major gene states. The colouring scheme for these 15 clusters is consistent with that used in (a), and the other clusters are shown in grey. (c) The characteristic gene expression of each gene state ranked by the likelihood in the RACIPE models. The four gene states that match the experimental data are highlighted by blue asterisks and are shown with their likelihoods. (The method to classify the gene states in the presence of external signal can be found in electronic supplementary material, §S8.)

Figure 4.

Comparison between the stemness GRN and the randomized networks (10 for Type I (a) and 10 for Type II (b)). Percentage of the RACIPE-predicted gene expression data matching each experimental gene state shown in figure 3a (right) for the stemness GRN and random networks. The details of the 10 Type I randomized networks and the 10 Type II networks can be found in electronic supplementary material, figure S9.

Figure 5.

Key parameters that are involved in the transitions among certain gene states. (a) A summary of the results depicted on top of the probability density map (figure 2c) of the RACIPE-generated gene expression data. In (a), along with each transition, the key parameters that have shifted the most have been marked. Red up arrows represent upregulation and blue down arrows represent downregulation. The mean of the normalized values for each parameter for the two corresponding gene states (x-axis for the first gene state, and y-axis for the second gene state) are shown in (b) (states 1 and 2) and (c) (states 1 and 6). The change of the parameters values between any of the transitions shown in (a) can be found in electronic supplementary material, figure S14.

Figure 6.

Hierarchical structure of the stemness GRN inferred from the perturbation analysis. (a) The Kullback–Leibler (KL) divergence between the probability distributions of the number of stable states for each RACIPE model computed before and after the knockout (KO) of each gene. ‘Oct4–Sox2*’ represents the removal of the protein complex Oct4–Sox2. (b) Similar to (a), but the KL divergences are between the distributions before and after removal of each regulatory link. (c) Schematic diagram of the stemness GRN highlighting the important genes and regulatory links. The larger the gene element and the thicker the regulatory link, the more important the component to the network behaviour, as inferred from the analyses in (a) and (b). (d) The hierarchical structure of the stemness GRN (left) is consistent with the two-step decision-making of mouse embryonic development (right). (e) The roadmap of stem cell differentiation inferred from the RACIPE simulations. All the original RACIPE models (WT) were treated by activating (↑) or inhibiting (↓) the maximum production rate of the corresponding genes by 50-fold. The probability of different gene states is proportional to the area in the pie chart.

Figure 7.

Analysis results elucidated by applying the extended RACIPE (RACIPE-wb) to the stemness GRN reveals consistent characterization. (a) Top panel: Diagram of the core stemness GRN highlighting the binding interactions between Oct4, Sox2, and the OCT4–SOX2 complex. Bottom panel: The main changes in the mathematical equations simulating the dynamics of Sox2 (X9), Oct4 (X7), OCT4–SOX2 (X8) to capture their binding/unbinding interactions. The full equations for RACIPE-wb are listed in electronic supplementary material, §S2. (b) The 2D probability density map of the results for the RACIPE-wb model projected onto the first two principal components. (c) Unsupervised HCA for RACIPE and RACIPE-wb, left and right, respectively. Clusters were identified using a probability cut-off of 0.005. The lines between clusters show the majority of clusters from the original RACIPE are also present in RACIPE-wb. Additionally, some of the clusters obtained by using RACIPE-wb were seen to be over/underrepresented as compared to results using the original RACIPE framework; the number of solutions belonging to that cluster is shown by the coloured vertical bars on the left and right of the middle which correspond the dendrograms of HCA for RACIPE and RACIPE-wb, respectively. (d) The KL divergence of RACIPE-wb distributions before and after knocking out a gene. ‘Oct4–Sox2*’ represents the removal of the protein complex Oct4–Sox2. (e) The KL divergence of RACIPE-wb distributions before and after the removal of a regulatory link. Also included are the blocking of binding between Oct4 and Sox2 and the blocking of the unbinding of OCT4–SOX2. (f) A schematic diagram depicting the relative importance of each gene and link as inferred by the analysis in (d–e).

Figure 8.

Schematic illustration of the revised Waddington's epigenetic landscape for stem cell differentiation. For each cell potency, the accessible cell types are shown by the attractors. Stem cells are induced by external signals toward differentiation along the valleys (highlighted by arrow lines with different colours) with the decrease of cell potency.