Cell Fate Decision as High-Dimensional Critical State Transition

Abstract

Cell fate choice and commitment of multipotent progenitor cells to a differentiated lineage requires broad changes of their gene expression profile. But how progenitor cells overcome the stability of their gene expression configuration (attractor) to exit the attractor in one direction remains elusive. Here we show that commitment of blood progenitor cells to the erythroid or myeloid lineage is preceded by the destabilization of their high-dimensional attractor state, such that differentiating cells undergo a critical state transition. Single-cell resolution analysis of gene expression in populations of differentiating cells affords a new quantitative index for predicting critical transitions in a high-dimensional state space based on decrease of correlation between cells and concomitant increase of correlation between genes as cells approach a tipping point. The detection of "rebellious cells" that enter the fate opposite to the one intended corroborates the model of preceding destabilization of a progenitor attractor. Thus, early warning signals associated with critical transitions can be detected in statistical ensembles of high-dimensional systems, offering a formal theory-based approach for analyzing single-cell molecular profiles that goes beyond current computational pattern recognition, does not require knowledge of specific pathways, and could be used to predict impending major shifts in development and disease.

Fig 1. Single-cell analysis of transcript expression during binary fate decision in EML cells.

(A) A progenitor EML cell population was stimulated with the cytokines erythropoietin (EPO) (left), interleukin-3/granulocyte macrophage-colony stimulating factor (IL-3/GM-CSF) (right), or with a combination of EPO and GM-CSF/IL-3 (center). Flow cytometry histograms of Sca1 surface expression were gated into Sca1^LOW (L), Sca1^MEDIUM (M), and Sca1^HIGH (H) fractions or subpopulations (green boxes) during fluorescence-activated cell sorting (FACS) of single cells at the indicated days for use in later analysis (Fig 2). At d3, further division to account for the extreme outliers (Lʹ, Hʹ)* indicates “rebellious cells” (see text). As previously reported [4], myeloid differentiation in EML cells is driven by the Sca1^HIGH fraction of cells and the global decrease of Sca1 expression is delayed, and can bounce back to intermediate state once cell have passed the bimodal (d3) stage. (B) For visualization of individual cells’ transcript expression patterns (of m = 17 genes), cells were projected onto a dimension-reduced state space spanned by the first three principal components (PC) following principal component analysis (PCA, see S1 Appendix). Each sphere represents a cell, colored according to treatment: untreated progenitors (grey); cells treated with EPO (red), cells treated with GM-CSF/IL-3 (blue); and combined-treated cells (purple). (C) To calculate a quasi-potential landscape for the three cell types for visualization of the idea of attractors as potential wells, a Gaussian filter with σ = 2 was applied to PC1 and PC2 coordinates of cells at d0 and d6 treated with EPO and GM-CSF/IL-3, leading to a smooth 2-dimensional distribution p. With the (quasi-)steady state assumption, the attractor landscape was visualized relative to a base level of 0 by − log(p +1). This time-invariant schematic is only a visual guide and does not model the landscape distortion during the bifurcation. Numerical data for the PCA graphs can be found in the supporting file S2 Data.

Fig 2. Critical transition during lineage commitment.

(A) Cell–cell correlation matrices displaying the Pearson correlation coefficient R(S^k, S^l) for all pairs of cells in states S^k and S^l (see S2 Appendix). R calculated for a set of 150 progenitor cells, 500 EPO-treated, 500 GM-CSF/IL-3-treated, and 450 combination-treated (COMB) cells from data used in Fig 1. Black squares (diagonal) emphasize the higher correlation between cells within the nominally same population. Two control genes (GAPDH and TBP) were excluded from this analysis. Lʹ, L, M, H, Hʹ indicate the Sca1 fractions shown in Fig 1: extremely low, low, medium, high, and extremely high level of Sca1 expression, respectively. (B) Average Pearson correlation coefficients for cell–cell pairs (left) and gene–gene pairs (center) as well as the state transition index I_c = 〈|R(g_i, g_j)|〉 / 〈R(S^k, S^l)〉 at various time points. Correlation coefficients were calculated for the central fractions/subpopulations in panel A(*). Error bars indicate SEM. (C) Gene–gene correlation matrices for the 17 genes of interest and the two endogenous control genes for the three treatments at various time points where correlation is indicated either by color (lower matrix triangle) or solid color segment in pie chart. Color values for magnitude of correlation coefficient for both matrices (A, C) are shown in color bar. Numerical data for the graphs in (B) and (C) can be found in S2 Data.

Fig 3. Intermediate stage of myeloid commitment along CD11b dimension exhibits destabilization of progenitor state.

(A) Flow cytometry dot plot of expression of Sca1 and CD11b upon treatment of the progenitor EML cells with GM-CSF/IL-3. Three distinct subpopulations on d3, designated, α, β and γ, in the tri-modal distribution of CD11b flow cytometry histogram underneath (red line, treated; blue line, untreated). (B) Cell–cell correlation for 72 progenitor cells and 48 cells from each of the α, β, and γ subpopulations, and gene–gene correlation for all 17 genes of interest and two endogenous control genes. Pearson correlation coefficient displayed as heatmap; same color scheme as in Fig 2. Bar graphs in box show the drastic increase of I_C for all the three subpopulations at d3 compared to the untreated progenitor cells, P(d0). I_C computed as in Fig 1. (C) Rescue by EPO of the “rebellious” = unintended γ subpopulation (pink curve) during myeloid differentiation. Three subpopulations (α, dark blue; β, light blue; γ, pink) were FACS sorted and antibodies were removed and stimulated with EPO. Total cell number and viability were quantified on day of sorting (d3) and four subsequent days. Viability was determined based on percent of cells excluding trypan blue. Each point represents average +/- STD for two biological replicates.

Fig 4. Critical Slowing down of state relaxation during fate commitment.

Apparent critical slowing down of relaxation and restoring of parental distribution of the sorted Sca1-low outlier fraction in the treated population. Clonal EML progenitor cells were stimulated (top) with GM-CSF/IL-3 or not (bottom), and cells with lowest 15% Sca1 expression were FACS-sorted one day after stimulation.

Fig 5. Whole-population transcriptome analysis reveals transient alternative program (rebellious cells).

(A) Sca1 surface expression population distribution in progenitor and cytokine-treated cells and transcriptomes of sorted subpopulations at indicated treatments/time points, displayed as GEDI self-organizing maps [44]. Progenitor EML cells were stimulated with EPO alone, with GM-CSF/IL-3 alone or with the combination of the two, and the Sca1-Medium (M) fractions (d0 and d6) and/or the Sca1-Low and -High subpopulations (d3) were FACS sorted for microarray analysis. (B) Hierarchical cluster analysis of the microarray-based transcriptomes of samples in A (columns, correspondence indicated by the green numbers) for a subset of the 17 genes analyzed in single-cell qPCR (rows). The entire experiment was performed twice in two laboratories with similar results. Color bar represents transcript expression values in log-scale (see Materials and Methods).

Fig 6. Epigenetic landscape model of symmetry-breaking bifurcation event.

The architecture and the specification of the individual interactions of the gene regulatory network (left) determines the topography of the epigenetic landscape (center) in which the elevation (formally, a quasi-potential U) [29] visualizes the relative stability of individual cell states, represented by the position with respect to the x-axis (state space). Thus, valleys represent stable attractor states. Fate commitment is induced by external differentiation signals, which initiate the deformation of the landscape and can be divided into two phases: first, the destabilization of the (meta)stable attractor of the progenitor cells (grey balls) and generation of a poised unstable state; and second, the opening of access to the destination attractors for both the intended and non-intended fate, biased toward the former by the differentiation signal. This allows the cells (balls) to descend to the new attractors. The disappearance of the progenitor attractor marks the critical transition (tipping point). As explained in S2 Appendix, the increase in cell–cell variation and in gene–gene correlation as the progenitor attractor destabilizes and cells spread and align along the particular direction (reaction coordinate) to exit the old attractor gives rise to a gradual increase of the index I_C (right) as cells approach the critical transition. This event coincides with lineage separation in state space and a major shift in cellular transcriptome. (Note, however, that there is no proof that I_C reaches a peak just at the tipping point.)