Simulated E/M Topography Displays Fractal Features.

A phenotypic landscape associated with our EMT/MET network can be reconstructed by performing a large number (M0=107) of simulations, starting from random initial conditions until the network reaches a steady state where si does not change. (No limit cycles are found; see SI Appendix for details.) In this way, we find a large number of distinct steady states that can be projected into a two-dimensional map using principal component analysis (PCA). We classify these steady states according to the expression of E-cadherin (CDH1), which we use as a reporter of the E/M phenotype (see Fig. 1B). The E/M map reconstructed from the model shows a clear separation between E and M states, with a boundary layer where E and M states coexist in very close proximity. A topographic representation of the stability of the states can be obtained by projecting H on the same two-dimensional map (Fig. 1C) showing that the boundary layer is more elevated with respect to pure E/M states, suggesting that those states are less stable. Furthermore, the map displays a very rough topography, with two main valleys separated by a large barrier populated by smaller and smaller valleys.

Given the sheer amount of distinct steady states (see Fig. 1D, Inset), we resorted to a statistical analysis and computed the probability distribution P(a) of the relative abundances of the states, where a is the fraction of times we find a given state. Fig. 1D shows that P(a) is a power law distribution, indicating that most of the states are very rarely found [when a is small, P(a) is large], but few states are found multiple times [when a is large, P(a) is small]. Alternative functional forms for P(a) are discussed in SI Appendix and shown in SI Appendix, Fig. S10. The presence of a power law is a signature of a scale-free fractal organization of the map, as is also apparent by the correlation matrix of the states. Fig 1E shows the presence of large correlated clusters subdivided into smaller and smaller clusters. In the physics of disordered systems, a hierarchical organization of the states is traditionally revealed by a broad distribution P(qαβ) of states’ overlap qαβ=∑i(siαsiβ)/N, measuring the similarity between two states {siα} and {siβ} (35). Hierarchical ground state structures have been observed in short-range Ising spin glasses (see refs. 36 and 37). When we restrict the sampling to low-H states, P(qαβ) displays a two-peak structure, indicating the presence of two classes of distinct and separate states (Fig. 1F), but when we consider all steady states, the overlap distribution becomes very broad, resembling the one observed in spin glasses, as noticed a long time ago for random Boolean networks (38⇓–40).

Simulated Phenotypic Transitions Reveal Scale-Free Stochastic Fluctuations.

Once the topography associated with the E/M landscape has been established, we investigate how the landscape changes when each one of the nodes is held fixed to si=±1, which simulates overexpression (OE) or knockdown (KD) of the corresponding gene (see SI Appendix for details). As an example, Fig. 2 A and B reports the one-dimensional projection of the topography under OE or KD of the SNAIL1 gene, a well-known inducer of the EMT. SNAIL1 OE leads to a rightward tilt of the landscape, favoring the M phenotype, while under SNAIL1 KD, the landscape tilts to the left, favoring the E state. This behavior is reminiscent of the effect of a magnetic field in a disordered magnet, where the free-energy landscape tilts in the direction of the field. If the network is initially in an E state, SNAIL1 OE can induce EMT, but the success rate and the trajectory crucially depend on the initial state (see Fig. 2C), with high-H states much more likely to undergo EMT than low-H states (see SI Appendix, Fig. S2). The variability in the outcome resulting from the OE/KD of a single gene can also be quantified by measuring the distribution of the number of nodes z affected by the process (see Fig. 2D). The distribution decays as power law P(z)∼z−τ, up to a cutoff value that increases with the H value of the initial state (see Fig. 2E), a further indication that high-H states are more susceptible to fluctuations (see also SI Appendix, Fig. S2). The avalanche exponent of the power law distribution is τ≃3/2, a value expected for mean-field avalanches in driven disordered systems (28).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

EMT/MET occurs with different probabilities through multiple paths. The model shows many forms of EMT/MET, and these occur with different probabilities. (A and B) One-dimensional PCA projection of the H landscape where (A) OE or (B) KD of SNAI1 tilts the landscape toward the M or E regions, respectively. (C) Transition map under SNAI1 OE. The model displays different forms of SNAI1-induced EMT. (D) The distribution of gene expression avalanches after individual KD/OE is a power law with exponent τ≃1.5. (E) The cutoff of the distribution depends on H, quantified here by quartiles, with high-H states producing larger avalanches. (F) EMT/MET probabilities under KD/OE conditions. The model lays out a nondeterministic picture of EMT/MET, where well-known factors such as SNAI1 (EMT) or KLF4 (MET) induce phenotypic transitions with higher probability (see Materials and Methods and SI Appendix, Fig. S2 for further details).

Using the model, it is possible to perform OE/KD on all of the nodes and estimate the probability of each node to induce EMT or MET (see Fig. 2F). Ranking the nodes as a function of their relevance for EMT, we recover well-known EMT inducers such as SNAIL1, ZEB1, or TGF-β and MET suppressors such as KLF4 and mir-200. The general pattern is that an inducer of EMT by OE also induces MET by KD and similarly for MET. We also simulate a transient version of OE/KD where a node is switched (si→−si), but it is then allowed to eventually relax back to its previous state. The results summarized in SI Appendix, Fig. S2 are similar to those obtained under stable OE/KD, for which the node variable is held fixed throughout the simulation, but the probability of EMT/MET is always smaller.

E/M Topography Inferred from Gene Expression Data Agrees with Simulations.

To confirm that the topographic representation of the E/M landscape obtained through the model provides an accurate representation of cellular phenotype, we examined the large cohort of gene expression data from human tissues provided by the GTEx project (41). To directly compare experimental data to the model, we designed a simple binarization strategy to decide whether a gene is expressed or not in a particular sample or cell. To calibrate the binarization scale, we used skin cells and fibroblasts as reference E and M states, respectively, and set a threshold based on the expression distribution of each gene in these two datasets (see Fig. 3A and SI Appendix). Genes whose expression was above the threshold were assigned to si=1 and otherwise to si=−1. The same threshold could then be used to binarize all of the 11,688 transcriptomes from different tissues present in the GTEx database.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Multitissue gene expression data display statistical features in agreement with simulations. (A) Illustration of the binarization process (see Materials and Methods for details). Gene-level expression data are casted into node-level binary data using binarization thresholds, computed using two reference samples (orange and black coloring). (B) Skin (orange) and fibroblast (black) samples from the GTEx project projected in PCA space. The E-cadherin expression probability in the model is shown with green (100%) to violet (0%) shades. Fibroblast samples tend to be in areas of very low E-cadherin expression probability. (C) Same as B but coloring the model steady states by average H. (D) Distribution of abundances, computed using all GTEx binarized samples and the 14 most relevant nodes (see SI Appendix and SI Appendix, Fig. S3 for details). (E) Clustering of 500 GTEx samples (all tissues), displaying a hierarchical structure qualitatively similar to that of the model (compare with Fig. 1E). (F) Overlap distribution over skin and fibroblast samples from the GTEx project (compare with Fig. 1F). (G) Overlap distribution over all GTEx samples (compare with Fig. 1G).

Using the topographic map of the E/M landscape constructed from simulations, we could then localize individual samples projecting their gene expression data on the map as shown in Fig. 3B. We then used the model to infer the stability of each phenotype by computing H associated to each state (Fig. 3C). When we plotted skin cells and fibroblasts on a two-dimensional map, we saw that they correctly fell into E or M regions, respectively (see Fig. 3B), but not all samples had the same value of H (see Fig. 3C). We used the same strategy to localize on the same topographic map the entire set of tissues present in the GTEx database (see SI Appendix, Figs. S4 and S5) and showed that they cover all of the available phase space. Assuming that the GTEx database contains an unbiased random sampling of all of the available states—which is a reasonable assumption given that the GTEx project provides multitissue gene expression data from healthy individuals only (41)—we analyzed the statistical properties of these states. As shown in Fig. 3D, the abundance distribution derived from GTEx data decays as a power law with an exponent that is very close to the one found numerically (compare with Fig. 1D and see SI Appendix and SI Appendix, Fig. S4 for technical details). Furthermore, clustering of the states showed a correlation matrix with hierarchical features that are in reasonable agreement with the prediction of the model (compare Fig. 3E with Fig. 1E). Finally, the overlap distribution displayed a two-peak structure when the statistics were restricted to fibroblasts and skin cells (Fig. 3F), while a single-peaked distribution was found when using all of the GTEx samples (Fig. 3G). This is in close agreement with the simulation results reported in Fig. 1 and confirms that experimental gene expression data give rise to a topographic landscape quantitatively similar to the one predicted by the model.

Tracing Bulk and Single-Cell RNA-Sequencing Trajectories Reveals the Nature of Hybrid E/M States.

The topographic representation of E/M states derived above can be used to visualize and interpret RNA-sequencing (RNA-seq) data obtained while the cells are undergoing phenotypic transformations. We first consider the classical example of TGF-β–induced EMT in a human lung adenocarcinoma cell line (42). Fig. 4A reports the trajectory of the states obtained from the bulk RNA-seq data recorded at different time points after TGF-β induction. As expected, the trajectory starts from the E region and crosses over to the M region of the map, as revealed by coloring the map according to the predicted expression of CDH1. Conversely, the trajectory obtained from RNA-seq data for DOX-induced MET during somatic cell reprogramming starts from the M valley and moves into the E valley of the landscape (30).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Single cells and bulk transcriptomic data yield trajectories through the E/M map with putative hybrid states lying on high-H regions. (A) Data from TGF-β–treated lung adenocarcinoma cell lines [GSE17708 (42)] yield a trajectory moving from the E to the M region. (B) Data from Dox-induced somatic cell reprogramming [GSE21757 (30)] display a reverse trajectory from M to E. Experimental data are shown as colored symbols with time course marked with arrows. The colored background depends on the ratio of steady states of the model that express E-cadherin at a given location in PCA coordinates, ranging from 100% (green) to 0% (violet). (C) Experimental data from single-cell embryonic-to-endoderm differentiation [GSE75748 (43)] move across the map as cells undergo EMT. The background color indicates the ratio of steady states that express E-cadherin or KLF4 (see SI Appendix, Fig. S7 for more markers). (D) Localization of single-cell gene expression data from tumor cells obtained from head and neck squamous cell carcinoma patients (45). All tumor cells correctly lie in the E region of the map, with high pEMT-scored cells located toward high-H areas. (E) The pEMT score correlates with H. Each gray dot represents a single cell. R and p denote the Pearson correlation coefficient and its associated P value (Student’s t test, two-tailed). The red line shows the average H.

Our methodology is even more revealing when applied to single-cell RNA-seq (scRNA-seq) data, as shown in Fig. 4C, which reports the time course of the states obtained from scRNA-seq data undergoing EMT during embryonic to endoderm differentiation (43) [see also SI Appendix, Fig. S6 illustrating MET during fibroblast to cardiomyocyte reprogramming in single-cell and bulk samples (44)]. As time goes on, cells originally in the E region transition to the M region between 24 and 36 h. After this time, even though EMT is apparently completed, the kinetic evolution of the cell population does not stop, and the region occupied by single-cell states shrinks. If we color the map by the predicted expression of other markers, we observe that the evolution moves cells in a low-KLF4 region (Fig. 4C; see also SI Appendix, Fig. S7 for similar maps for other markers). Hence, when applied to scRNA-seq data, our method can reveal subtle features associated with phenotypic transformations.

This last point is best illustrated by analyzing recent data (6,000 single cells) obtained from 18 head and neck squamous cell carcinoma patients (45). The original analysis revealed the presence of an aggressive cancer cell population, associated with metastasis and poor prognosis, described as partial-EMT (pEMT) (45). Classification of cells as pEMT was based on a pEMT score computed from the expression values of a set of 100 genes (45), none of which directly maps into nodes of our model. It is thus particularly remarkable to see that the projection of the scRNA-seq data on our map reveals that tumor cells are correctly located in the E region of the map and cells with high pEMT score are typically located on higher ground with respect to low pEMT cells (see Fig. 4D). This is corroborated by the strong correlation between H and the pEMT score, as reported in Fig. 4E.