Student or educator? Sign up with your school email and get 6 months of Pro free.
Cited by 350

The DOI matches its Crossref record; 10 of the 10 references it cites were found in a registry or an index.

Cell Fate Decision as High-Dimensional Critical State Transition

Mojtahedi, Mitra, Skupin, Alexander, Zhou, Joseph, Castaño, Ivan G., Leong-Quong, Rebecca Y. Y., Chang, Hannah, Trachana, Kalliopi, Giuliani, Alessandro, Huang, Sui·2016·PLOS Biology

Published version of a preprint

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.

Overview

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Integrity signals

Registry

Crossref · 10.1371/journal.pbio.2000640

Cites

10 of 10 references verified0 citation errors · 0 not checkable

Statistics

Figures

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 Sca1LOW (L), Sca1MEDIUM (M), and Sca1HIGH (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 Sca1HIGH 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 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 Sca1LOW (L), Sca1MEDIUM (M), and Sca1HIGH (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 Sca1HIGH 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.CC BY · PubMed Central
Fig 2 Critical transition during lineage commitment. (A) Cell–cell correlation matrices displaying the Pearson correlation coefficient R(Sk, Sl) for all pairs of cells in states Sk and Sl (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 Ic = 〈|R(gi, gj)|〉 / 〈R(Sk, Sl)〉 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 2 Critical transition during lineage commitment. (A) Cell–cell correlation matrices displaying the Pearson correlation coefficient R(Sk, Sl) for all pairs of cells in states Sk and Sl (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 Ic = 〈|R(gi, gj)|〉 / 〈R(Sk, Sl)〉 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.CC BY · PubMed Central
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 IC for all the three subpopulations at d3 compared to the untreated progenitor cells, P(d0). IC 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 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 IC for all the three subpopulations at d3 compared to the untreated progenitor cells, P(d0). IC 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.CC BY · PubMed Central
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 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.CC BY · PubMed Central
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 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).CC BY · PubMed Central
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 IC (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 IC reaches a peak just at the tipping point.)
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 IC (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 IC reaches a peak just at the tipping point.)CC BY · PubMed Central

Verification

Claims the paper makes

The claims and statistics of this paper haven’t been checked yet.

Run the check — free, no account needed. It reads the open-access text and checks the sources it cites.

A few free checks a day · papers already verified open instantly · results are public

Claims other papers make about this one

Each is checked against this paper’s own text: does it say what the citing paper says it does?

Self-check

Percentages, sums and reported differences re-computed from the paper’s own numbers.

The self-check needs the paper’s full text.

Citations

Cites

Show all 10 references

Research works printed but not linked: 38

  1. BD Macarthur, 2009, Systems biology of stem cell fate and cellular reprogramming, Nat Rev Mol Cell Biol, 10, 672, 10.1038/nrm2766 · DOI · Paper page
  2. S Huang, 2010, Cell lineage determination in state space: a systems view brings flexibility to dogmatic canonical rules, PLoS Biol, 8, 1
  3. J Wray, 2010, The ground state of pluripotency, Biochem Soc Trans, 38, 1027, 10.1042/BST0381027 · DOI
  4. HH Chang, 2008, Transcriptome-wide noise controls lineage choice in mammalian progenitor cells, Nature, 453, 544, 10.1038/nature06965 · DOI
  5. SR Hough, 2009, A continuum of cell states spans pluripotency and lineage commitment in human embryonic stem cells, PLoS ONE, 4, e7708, 10.1371/journal.pone.0007708 · DOI
  6. S Kauffman, 1969, Homeostasis and differentiation in random genetic control networks, Nature, 224, 177, 10.1038/224177a0 · DOI · Paper page
  7. S Huang, 2007, Bifurcation dynamics of cell fate decision in bipotent progenitor cells, Dev Biol, 305, 695, 10.1016/j.ydbio.2007.02.036 · DOI · Paper page
  8. D Kaplan, 1995, Understanding Nonlinear Dynamics, 10.1007/978-1-4612-0823-5 · DOI
  9. C Trefois, 2015, Critical transitions in chronic disease: transferring concepts from ecology to systems medicine, Curr Opin Biotechnol, 34, 48, 10.1016/j.copbio.2014.11.020 · DOI
  10. V Chickarmane, 2008, A computational model for understanding stem cell, trophectoderm and endoderm lineage determination, PLoS ONE, 3, e3478, 10.1371/journal.pone.0003478 · DOI
  11. C Duff, 2011, Mathematical modelling of stem cell differentiation: the PU.1-GATA-1 interaction, J Math Biol, 64, 449, 10.1007/s00285-011-0419-3 · DOI
  12. JE Ferrell, 2012, Bistability, bifurcations, and Waddington's epigenetic landscape, Current biology, 22, R458, 10.1016/j.cub.2012.03.045 · DOI
  13. MJ Jolly, 2014, Towards elucidating the connection between epithelial-mesenchymal transitions and stemness, J R Soc Interface, 11, 20140962, 10.1098/rsif.2014.0962 · DOI
  14. C. Furusawa, 2012, A dynamical-systems view of stem cell biology, Science, 338, 215, 10.1126/science.1224311 · DOI
  15. M Pal, 2014, Non-genetic heterogeneity, criticality and cell differentiation, Phys Biol, 12, 016001, 10.1088/1478-3975/12/1/016001 · DOI
  16. L Chen, 2012, Detecting early-warning signals for sudden deterioration of complex diseases by dynamical network biomarkers, Sci Rep, 2, 18, 10.1038/srep00342 · DOI · Paper page
  17. S Tsai, 1994, Lymphohematopoietic pro-genitors immortalized by a retroviral vector harboring a dominant-negative ret-inoic acid receptor can recapitulate lymphoid, myeloid, and erythroid develop-ment, Genes Dev, 8, 2831, 10.1101/gad.8.23.2831 · DOI
  18. F Paul, 2015, Transcriptional Heterogeneity and Lineage Commitment in Myeloid Progenitors, Cell, 163, 1663, 10.1016/j.cell.2015.11.013 · DOI · Paper page
  19. D Grun, 2014, Validation of noise models for single-cell transcriptomics, Nature methods, 11, 637, 10.1038/nmeth.2930 · DOI · Paper page
  20. GK Marinov, 2014, From single-cell to cell-pool transcriptomes: stochasticity in gene expression and RNA splicing, Genome Res, 24, 496, 10.1101/gr.161034.113 · DOI · Paper page
  21. JX Zhou, 2015, Quasi-potential landscape in complex multi—stable systems, J R Soc Interface, 12, 1
  22. CH Waddington, 1956, Principles of Embryology, 10.5962/bhl.title.7217 · DOI
  23. AN Gorban, 2010, Correlations, risk and crisis: From physiology to finance, Phys A Stat Mech. its Appl, 389, 3193, 10.1016/j.physa.2010.03.035 · DOI
  24. J Wang, 2010, The potential landscape of genetic circuits imposes the arrow of time in stem cell differentiation, Biophys J, 99, 29, 10.1016/j.bpj.2010.03.058 · DOI
  25. A Giuliani, 2014, Statistical Mechanics of Gene Expression Networks: Increasing Connectivity as a Response to Stressful Condition, Adv Syst Biol, 3, 1
  26. A Giuliani, 2014, Invariant features of metabolic networks: a data analysis application on scaling properties of biochemical pathways, Physica A, 337, 157
  27. T Enver, 1998, Do stem cells play dice?, Blood, 92, 358, 10.1182/blood.V92.2.348 · DOI
  28. T Enver, 2009, Instructions writ in blood, Nature, 461, 10.1038/461183a · DOI · Paper page
  29. RL Coffman, 1999, Instruction, selection, or tampering with the odds?, Science, 284, 1283, 10.1126/science.284.5418.1283 · DOI
  30. RE Callard, 2007, Decision-making by the immune response, Immunol. Cell Biol, 85, 300, 10.1038/sj.icb.7100060 · DOI
  31. AJ Veraart, 2012, Corrigendum: Recovery rates reflect distance to a tipping point in a living system, Nature, 484, 404, 10.1038/nature11029 · DOI
  32. D Metcalf, 1999, Stem cells, pre-progenitor cells and lineage-committed cells: Are our dogmas correct?, Ann N Y Acad Sci, 872, 289, 10.1111/j.1749-6632.1999.tb08473.x · DOI
  33. M Ogawa, 1989, Hemopoietic stem cells: stochastic differentiation and humoral control of prolifera-tion, Environ Health Perspect, 80, 199, 10.1289/ehp.8980199 · DOI
  34. GS Eichler, 2003, Gene Expression Dynamics Inspector (GEDI): for integrative analysis of expression profiles, Bioinformatics, 19, 2321, 10.1093/bioinformatics/btg307 · DOI
  35. D Metcalf, 1998, Lineage commitment and maturation in hematopoietic cells: the case for extrinsic regulation, Blood, 92, 345, 10.1182/blood.V92.2.345b · DOI
  36. C Marr, 2016, Single-cell gene expression profiling and cell state dynamics: collecting data, correlating data points and connecting the dots, Curr Opin Biotechnol, 39, 207, 10.1016/j.copbio.2016.04.015 · DOI · Paper page
  37. JX Zhou, 2011, Understanding gene circuits at cell-fate branch points for rational cell reprogramming, Trends Genet, 27, 55, 10.1016/j.tig.2010.11.002 · DOI
  38. P Liu, 2007, Imperfect transcritical and pitchfork bifurcations, J Funct Anal, 251, 573, 10.1016/j.jfa.2007.06.015 · DOI

Cited by

Metadata

Authors
Mojtahedi, M., Skupin, A., Zhou, J., Castaño, I. G., Leong-Quong, R. Y. Y., Chang, H., Trachana, K., Giuliani, A., & Huang, S.
Venue
PLOS Biology
Published
2016-12-27
Volume
14
Issue
12
Pages
e2000640
ISSN
1545-7885
Cited by
350
References
10
Language
en

About these checks

Registry, retraction, Cites and statistics results are computed, not judged, with one exception: a printed reference that no index lists is read once by a model before it counts as a citation error. Claims are extracted by a model and kept only when their quote appears verbatim in the text; each cited source is then checked, with a model’s help, for whether it says what the claim relies on. Those verdicts can be wrong — read the claim’s quote and the cited source before relying on one.

Statistics and arithmetic are recomputed, not judged.

Confidence bands are shown only for fields calibrated on our annotated test set. None is calibrated yet, so no band is shown.

What counts as a claim ›How we measure confidence ›How we check this ›

Are you an author, or do you see a problem with this page? Request a review

Data sources

Bibliographic data: Semantic Scholar, OpenAlex, Crossref.

Contains information from the Semantic Scholar Academic Graph Datasets, made available under the ODC Attribution License (ODC-BY 1.0).

Semantic Scholar’s data is described in: Kinney, R., et al. (2023). The Semantic Scholar Open Data Platform. arXiv:2301.10140.