Hybrid E/M Phenotype(s) and Stemness: A Mechanistic Connection Embedded in Network Topology

Abstract

Metastasis remains an unsolved clinical challenge. Two crucial features of metastasizing cancer cells are (a) their ability to dynamically move along the epithelial–hybrid–mesenchymal spectrum and (b) their tumor initiation potential or stemness. With increasing functional characterization of hybrid epithelial/mesenchymal (E/M) phenotypes along the spectrum, recent in vitro and in vivo studies have suggested an increasing association of hybrid E/M phenotypes with stemness. However, the mechanistic underpinnings enabling this association remain unclear. Here, we develop a mechanism-based mathematical modeling framework that interrogates the emergent nonlinear dynamics of the coupled network modules regulating E/M plasticity (miR-200/ZEB) and stemness (LIN28/let-7). Simulating the dynamics of this coupled network across a large ensemble of parameter sets, we observe that hybrid E/M phenotype(s) are more likely to acquire stemness relative to “pure” epithelial or mesenchymal states. We also integrate multiple “phenotypic stability factors” (PSFs) that have been shown to stabilize hybrid E/M phenotypes both in silico and in vitro—such as OVOL1/2, GRHL2, and NRF2—with this network, and demonstrate that the enrichment of hybrid E/M phenotype(s) with stemness is largely conserved in the presence of these PSFs. Thus, our results offer mechanistic insights into recent experimental observations of hybrid E/M phenotype(s) that are essential for tumor initiation and highlight how this feature is embedded in the underlying topology of interconnected EMT (Epithelial-Mesenchymal Transition) and stemness networks.

1. Introduction

Metastasis is the deadliest process in cancer, and it is involved in over 90% of all cancer-related deaths. It is a dynamic multi-stage cascade of events involving the initial detachment of cells from the primary tumor mass, dissemination of cancer cells, their exit (extravasation) at many distant organs, and finally, their ability to colonize distant organs through tumor outgrowth [1]. It is a challenging process for cells with extremely high attrition rates (>99%) such that most cells die in circulation or are unable to adapt to the foreign biochemical and/or biophysical ecology of distant organs [2,3]. Thus, phenotypic plasticity—the ability of disseminating cancer cells to adapt their phenotypes reversibly in response to their dynamic microenvironments—has been regarded as a hallmark of metastasis-initiating cells [4]. Phenotypic plasticity exists at multiple interconnected axes (a) epithelial–mesenchymal plasticity (EMP), (b) metabolic plasticity, and (c) plasticity between a cancer stem cell (CSC) and non-CSC state (i.e., stemness), among others [5]. Understanding such functional inter-dependencies from a dynamical systems level can lead to deciphering the organizing principles of underlying regulatory network modules and, eventually, therapeutic advances that can target cancer cell adaptability/plasticity [6,7].

Out of the various possible pairwise couplings among these different axes of plasticity, the one between EMP and stemness has been investigated most thoroughly through in vitro, in vivo, and in silico approaches. Initial reports, which treated EMP as a binary (all-or-none) process where cells switch between epithelial and mesenchymal phenotypes, suggested that a “full” epithelial–mesenchymal transition (EMT) was associated with increased stemness (tumor initiation potential) [8,9]. This view was challenged by later studies demonstrating that a “full” EMT was dispensable for acquiring stemness, or also perhaps an obstacle for the same [10,11]. A more nuanced recent view highlights that EMP involves multiple stable hybrid E/M states in addition to “fully” E or “fully” M ones and that these hybrid states may be maximally stem-like [12,13,14,15,16,17]. In squamous cell carcinoma cells that were categorized into six phenotypes along the EMP spectrum, stemness was shown to be acquired in earlier stages of EMP, and it did not increase as cells progressed towards a “fully” mesenchymal phenotype. Intriguingly, the metastatic potential was found to be the maximum for hybrid E/M states [16]. Furthermore, in breast cancer cells, a transition from a hybrid E/M phenotype to a “fully” mesenchymal one, as driven by constitutive ectopic expression of EMT-inducing transcription factor (EMT-TF) ZEB1, led to the loss of tumorigenicity in vitro and in vivo [13]. Three-dimensional assays to investigate collective cell invasion in breast cancer organoids showed that most leader cells co-expressed epithelial, mesenchymal, and CSC markers [18], further strengthening the association between hybrid E/M phenotype(s) and stemness. One possible reason for the higher stemness of hybrid E/M cells may be their increased propensity to give rise to epithelial and mesenchymal cells [19,20] in addition to self-renewal, the two critical traits of stem cells in development and homeostasis. While this association of hybrid E/M phenotype with enhanced stemness has been consistently reported across cancer types [21], a mechanistic underpinning facilitating the association’s emergence remains unclear.

We have previously developed mechanism-based mathematical models that investigate the coupled nonlinear dynamics of EMP and stemness modules—miR-200/ZEB and LIN28/let-7, respectively. These models have predicted that while the most likely positioning of the “stemness window” is around the mid-point of the “EMP axis” (i.e., hybrid E/M phenotypes) [15], the “window” is dynamic in nature and could move to either more epithelial or more mesenchymal ends too [22,23]. A caveat of this analysis is that it was performed on a few specific parameter sets; thus, its applicability to explain the diverse and mounting experimental evidence across cancer types, as mentioned above, remains limited.

Here, we have simulated the coupled dynamics of miR-200/ZEB and LIN28/let-7 circuit over a wide range of biologically relevant parameter sets to interrogate whether the association of a hybrid E/M phenotype with stemness is a trait embedded in the topology of this coupling itself. Our results suggest that the overlap of hybrid E/M and stemness features is an outcome of coupled network topology between the miR-200/ZEB and LIN28/let-7 feedback loops. Moreover, upon extending the network to include various phenotypic stability factors (PSFs), such as GRHL2, OVOL1/2, and NRF2, we notice that this association is maintained. Thus, our results unravel the underlying design principles of the EMP–stemness association and predict that this interconnection is likely to be seen across multiple carcinomas.

2. Results

2.1. Gene Regulatory Network Underlying E/M Plasticity and Stemness Enables the Existence of Four E/M Phenotypes

To unravel the mechanistic underpinnings of the association between stemness and epithelial–mesenchymal plasticity (EMP), we considered the regulatory interactions among key factors implicated in governing the E/M plasticity (ZEB/miR-200) and the stemness characteristics (LIN28/let-7) (Figure 1A). E/M plasticity is at least in part controlled by a mutually inhibitory feedback loop between the transcription factor family ZEB and the microRNA family miR-200 [24,25]. Moreover, ZEB can self-activate through ESRP1 and/or CD44/HA signaling [26,27,28]. This network can give rise to three phenotypes: epithelial (high miR-200, low ZEB), mesenchymal (low miR-200, high ZEB), and hybrid E/M (medium miR-200, medium ZEB) [29]. This network can be influenced by other EMT-TFs such as SNAIL which self-inhibits transcriptionally and activates the expression of ZEB, and suppresses the transcription of the microRNA miR-200 [29]. Similarly, the “stemness” of a cell is primarily controlled by a mutually inhibitory feedback loop between RNA-binding factor family LIN28 and microRNA family let-7 [30]. Both LIN28 and let-7 are known to self-activate through direct and/or indirect mechanisms [31,32,33], and the transcription factor NF-kB activates the expression of both of them [15]. LIN28 regulates the levels of OCT4 [34], whose intermediate levels are considered to be optimal for attaining stemness [35,36,37]. Interestingly, these two modules—ZEB/miR-200 and LIN28/let-7—can influence each other, such that let-7 can target ZEB translation [33], and microRNA family miR-200 can inhibit LIN28 translation [38,39]. Put together, this integrated network of EMP and stemness regulators can explain experimental observations on LIN28 driving EMP [40] and ZEB/miR-200 influencing stemness [41].

All these regulatory links in the abovementioned network may be active, with different strengths in individual cells. Further, in cancer cells, the mutational background can alter the stability and/or production rate of various species, for instance, mutant p53 gain of function, triggering ZEB1 [42]. Thus, a rigorous analysis of the emergent dynamics of the coupled EMP–stemness network would require a framework that can decode the features robust to parameter variation of a specific network topology. To achieve this goal, we simulated the coupled EMP–stemness network through random circuit perturbation RACIPE [43] (see Section 4). RACIPE takes network topology as its input, converts it into a set of coupled ordinary differential equations (ODEs) where the influence of each activation or inhibition regulatory link is represented by a shifted Hill function. Parameters corresponding to this set of equations are sampled randomly within a biologically relevant range. Therefore, it generates an ensemble of ODE models, each with a different parameter set, thus representing the intrinsic variability in kinetic parameter values in a given cell population. The set of steady-state values obtained from this ensemble can be then used to identify robust dynamical signatures emerging from this topology through appropriate statistical analysis (see Section 4).

First, we plotted a heatmap of all the steady-state values obtained via RACIPE for this coupled network and performed agglomerative hierarchical clustering (Figure 1B) on it. We observed the possible existence of four major clusters, indicated in the dendrogram plotted using Ward’s minimum variance criterion (see Section 4) (Figure S1A). We next performed K-means clustering for varied values of K and observed a peak in the average silhouette width at K = 4, endorsing the existence of four clusters in RACIPE solutions (Figure S1B). This trend was supported by all other cluster quality metrics we used (Figure S1C–E). We plotted histograms of ZEB levels for the four clusters identified in the heatmap and assigned them phenotypic labels accordingly (Figure 1C,D). The cluster with the minimum median ZEB levels is referred to as an epithelial (e) phenotype, while the cluster with the maximum ZEB median levels is referred to as a mesenchymal (m) phenotype. The other two clusters had intermediate ZEB levels, being categorized as E/M hybrids, and are assigned the labels of hybrid epithelial (he) and hybrid mesenchymal (hm). Furthermore, we performed principal component analysis (PCA) on all steady-state solutions (Figure 1E; each color represents a distinct cluster identified by K-means). Principal component 1 and 2 explain 47% and 19% of the total variance, respectively, and show roughly four clusters in the PCA plot where the two hybrid clusters (blue and pink) are sandwiched between the epithelial (purple) and the mesenchymal (yellow) clusters (Figure 1E). PC1 can be thought of as an approximate EMT axis based on the coefficients for miR-200 and ZEB. Interestingly, while miR-200 levels are very different between e and m states, its levels are similar for e and he clusters, and for hm and m clusters (Figure S2A and Figure 1D, Table S7), thus indicating that miR-200 may be a more coarse-grained readout of the EMT state of a cell. Corroborative trends, i.e., the highest levels of let-7 and lowest levels of LIN28 in the e cluster, and lowest levels of let-7 and highest levels of LIN28 in the m cluster, were observed (Figure 1D and Figure S2B). The levels of NF-kB and SNAIL levels were found to be almost similar across the clusters (Figure S2C,D), which is not surprising because neither NF-kB nor SNAIL are influenced by any of the other four nodes (miR-200, ZEB, let-7, LIN28) in the coupled EMP–stemness network.

2.2. Stemness Characteristics are Enriched for in the Hybrid E/M Phenotypes

The gene regulatory network involving the interplay between E/M plasticity and stemness can give rise to different phenotypes (steady states). Next, we investigated whether these states can co-exist for certain parameter sets, i.e., is the network capable of exhibiting multistability? Out of the ensemble of parameter sets sampled by RACIPE, we observed that a subset of parameter sets can give rise to a combination of two or more steady states per parameter set (bistability: 37.3 ± 0.95%, tristability: 24.2 ± 0.09%, etc.), while the others converged to a single steady state for the given system of ordinary differential equations (monostability: 29.6 ± 0.70%). A closer inspection of the monostable solutions revealed epithelial (e) and mesenchymal (m) states as the most predominant solutions (Figure 2A), indicating that epithelial and mesenchymal phenotypes are more likely to exist relative to the hybrid ones, thereby reminiscent of observations based on the modeling of larger EMP regulatory networks too [44,45]. This trend is also observed when considering all the solutions together irrespective of the number of states (Table S5). Consistently, among all possible bistable solutions, the most predominant one is the phase {e, m}, i.e., the co-existence of epithelial and mesenchymal states (Figure S3A). The least frequent bistable phase is {he, hm}, i.e., the one with the co-existence of both hybrid phenotypes, while the frequency of bistable phases containing one hybrid phenotype and one of the epithelial or mesenchymal ones—{e, hm}, {he, m}, {e, he}, {m, hm}—is intermediate (Figure S3A). Similarly, among all possible tristable solutions, the phases containing both epithelial and mesenchymal states—{e, he, m} and {e, hm, m}—were the most predominant (Figure S3B). Put together, these results indicate that the coupled EMP–stemness regulatory network is multistable and therefore allows for cells to switch among the epithelial, mesenchymal, and the hybrid phenotypes due to stochastic perturbations.

Next, we investigated the stemness traits of different cell states along the E/M spectrum. Intermediate levels of OCT4 have been reported to confer maximal stemness to the cells [35,36,37]. Given that OCT4 is a direct target of LIN28, we defined the “stemness window” to be a region around median LIN28 expression levels in all steady-state solutions combined and the median averaged over all the replicates of the coupled EMP–stemness circuit. Specifically, we considered the median ± one interquartile range of the LIN28 distribution to be the “biological range” of LIN28, and the middle 30% of this range centered around its median is considered to be the “stemness window”, which might allow for a population enriched in the stem-like phenotype (Figure 2B). To test the robustness of this specifically chosen stemness window, we need a semi-objective heuristic for increasingly wider “stemness windows” centered at the mid-point of the range chosen. When the size of the “stemness window” is considered to be more than one-third of the range of LIN28 levels, we see a qualitative change in the distribution of phenotypes lying in this window. This analysis underscores that hybrid E/M phenotypes are predominantly associated with stemness unless the “stemness window” starts getting too large to be of a biologically implausible size (Figure S3C).

Next, we quantified the probability that a particular solution lies within the “stemness window” given that it belongs to a particular phenotype along the E–M spectrum (p1). Interestingly, we found that hybrid (he and hm) phenotypes are more likely to be stem-like than either extreme (e and m) E/M phenotypes (Figure 2C). Alternatively, we computed the probability that a solution belongs to a particular phenotype along the E–M spectrum, given that it already lies in the “stemness window” (p2). Congruently, we observed that hybrid phenotypes (he and hm) are more enriched for in the stemness region than either epithelial (e) or mesenchymal (m) solutions, although the contribution of e and m cell states to “stemness window” cannot be ignored (Figure 2D). Put together, these observations suggest that although the hybrid phenotypes he and hm are more likely to be associated with a stem-like characteristic, stemness is not exclusive to them. Thus, a subset of “pure” epithelial and/or mesenchymal cell states can also be potentially stem-like, depending on parametric combinations. Similar trends are observed when analyzing the multistable systems only as well (i.e., excluding all solutions obtained from monostable parameter sets), endorsing the robustness in trends (Supplementary Section S4).

Having established a strong association between the hybrid E/M phenotypes and stemness, we sought to understand the mechanistic underpinnings that can facilitate this association. Through RACIPE, we simulated a hypothetical “uncoupled” circuit where the inhibitory links from miR-200 to LIN28 and those from let-7 to ZEB are absent. The ensemble of solutions obtained for this uncoupled circuit also has four clusters (see Supplementary Section S2.2). For the uncoupled circuits, when we plot the clusters in a ZEB–LIN28 plane (indicative of an EMP–stemness plane with ZEB being a proxy for the E/M axis and LIN28 as a proxy for stemness), we observed that the mid-points of these four clusters were arranged in a square—the four vertices being [high ZEB, high LIN28], [low ZEB, high LIN28], [high ZEB, low LIN28], and [low ZEB, low LIN28]—demonstrating the independence of the ZEB and LIN28 levels (Figure S4A). The [high ZEB, high LIN28] and [low ZEB, low LIN28] clusters correspond to m and e states, respectively, and the location of the former on the two-dimensional plane for the coupled circuit overlapped with that of the uncoupled circuit. However, in the case of the coupled circuit, the square geometry is disturbed; instead, more of a rhombus-like shape emerged due to the shift of the mid-points of the [low ZEB, high LIN28] and [high ZEB, low LIN28] clusters which now tend to have intermediate values of both. This analysis suggests that coupling EMP and stemness modules enables a higher likelihood of both hybrid (he, hm) clusters to be more stem-like (Figure S4A). This observation can be mechanistically interpreted by understanding the relative “activity” of the two coupling links (miR-200 inhibiting LIN28 and let-7 inhibiting ZEB). In the [high ZEB, high LIN28] state, the levels of let-7 and miR-200 are relatively low, thus, the links are effectively “inactive”. In the case of [low ZEB, high LIN28] and [high ZEB, low LIN28], one link is much more active than the other, hence we see relocations mostly along one of the two axes of the LIN28–ZEB plane. For the [low LIN28, low ZEB] state, both links are approximately equally “active”, and therefore the cluster migrates along to even lower values of both LIN28 and ZEB.

To validate this hypothesis, we quantified the strength of coupling links—the net effect was strongest for e (i.e., [low LIN28, low ZEB]) and weakest for m (i.e., [high ZEB, high LIN28]) (Figure S4B). As a control case for the link strength metric, we quantified the degree of asymmetry in the strength of the two links within EMP and stemness modules individually (see Section 4). We observed that, as expected, the effective “suppression” (the expression free topological asymmetry) of ZEB on miR-200 was the highest in the case of hm and m states relative to e and he states (Figure S4D). On the other hand, we observed that the suppression of LIN28 on let-7 was higher in the case of he and m and lower in the e and hm phenotypes (Figure S4C).

2.3. The Effect of Phenotypic Stability Factors (PSFs) on the Phenotypic Landscape of EMP

Previous work, including ours, has shown that various factors, such as GRHL2, OVOL1/2, and NRF2, can stabilize the hybrid E/M phenotype of cells [46,47,48,49]. Their knockdown in H1975—a stable hybrid E/M lung cancer in vitro model—pushed cells towards a complete EMT, thus, they were classified as “phenotypic stability factors” (PSFs) for hybrid E/M phenotypes. For a specific set of parameters, incorporating GRHL2 also promoted the association between hybrid E/M states and stemness [46]. Thus, we investigated whether the hybrid E/M–stemness interconnection was maintained for multiple PSFs across an ensemble of parameter sets.

We simulated via RACIPE the coupled EMP–stemness networks after incorporating these three PSFs, one at a time (Figure 3A). As compared to the circuit without these PSFs (i.e., the “base” coupled EMP–stemness network shown in Figure 1A), the circuits including PSFs had a higher frequency of parameter sets leading to tristability and quadrastability, at the expense of a lower frequency of parameter sets primarily leading to monostability and even bistability in the case of the GRHL2 circuit (Figure 3B). In general, all these circuits showed a trend towards increasing multistability, with the strongest effect seen for the GRHL2 circuit (Figure 3B). However, incorporating the PSFs altered neither the optimal number of clusters as identified by average silhouette widths (Figure 3C) nor the position of clusters on the ZEB–LIN28 plane that showed an almost complete overlap with those seen for the “base” circuit (Figure 3D). Therefore, the introduction of the PSFs does not seem to disrupt the core structure of the steady-state solutions, as obtained earlier.

In addition to being thought of as PSFs, GRHL2 and OVOL1/2 have been considered as Mesenchymal-Epithelial Transition (MET)-inducing transcription factors (MET-TFs), as their overexpression in carcinomas is capable of upregulating the expression levels of E-cadherin and/or revert EMT [50,51,52,53,54] as well as other EMT-associated traits, such as anoikis resistance, metabolic reprogramming [55], and immune evasion [56]. This observation is consistent with their knockdown known to drive a “full” EMT in both cancer cell lines and in developmental contexts [46,49]. Thus, with the goal to identify the effect of overexpression and downregulation of these PSFs on relative proportions of the different EMPs, we quantified the changes in the frequency of different monostable solutions—e, m, he, and hm—upon a 10-fold over-expression (oe10) and a 10-fold down-expression (de10) of these PSFs. Computationally, this provides a dynamically consistent way to characterize the effect attributable to any node on the network dynamics and RACIPE implements this perturbation by altering the rate of production of the specific node (see Section 4). We observed that increasing the expression of any one of the three PSFs enriched the epithelial (e) phenotype at the expense of the mesenchymal (m) phenotype. Similar, albeit weaker, trends were seen for an enrichment of the he phenotype at the expense of the hm phenotype (Figure 3E). This effect was more pronounced in the GRHL2 case in comparison to either OVOL or NRF2 (Figure 3E). This trend is also present overall not only for only monostable solutions, but also for all parameter sets considered together (see Supplementary Section S2.4).

Moreover, a closer inspection of steady-state solutions from bistable parameter sets showed a consistent trend across PSFs in the case of their overexpression: a most remarkable decrease in the proportion of the {hm, m} phase with the most remarkable increase in the proportion of the {he, e} phase (Figure S5A–C). Similarly, the tristable states also showed an agreement with the trend of enriching “relatively epithelial” phases, with maximum enrichment of the {e, he, m} phase and a maximum depletion of the {e, hm, m} phase upon overexpression of the PSFs (Figure S5D–F), with unremarkable changes in the other phases (see Supplementary Section S2.3 for details on the multistable distributions). Put together, these results suggest that incorporating any of these PSFs (GRHL2, OVOL, NRF2) tends to maintain the four clusters seen in the EMP–stemness coupled network earlier, and their over-expression can serve as a way to at least inhibit the progression of a complete EMT. Further analyses regarding the association of multistability and stemness can be obtained in Supplementary Section S4.

2.4. Networks Incorporating Phenotype Stability Factors (PSFs) Preserve the Association of the Hybrid E/M States with Stemness

Next, we evaluated how the association of EMP phenotypes with stemness is altered upon including PSFs. As earlier, we explored two kinds of association of the hybrid E/M phenotypes and a stem-like state: the probability of being a stem-like state given that a solution belongs to a particular E/M phenotype (p1) and the probability of belonging to a particular E/M phenotype given that the solution is a stem-like state (p2). On comparing p1 and p2 for all circuits with the two hybrid clusters (he, hm) taken together and the two non-hybrid clusters (e, m) taken together, we observe that the hybrids are much more likely to be stem-like (Figure 4A) as well as constitute a larger proportion of the stemness-enriched pool of solutions (Figure 4B), across all cases of PSFs. This observation suggests that the association of hybrid E/M phenotypes with stemness is a dynamical trait encoded in the core network topology, and this feature is maintained upon including PSFs.

A closer look at the changes in p1 and p2 values upon over-expression or downregulation of PSFs revealed small variable minor changes for p1 (Figure S6A–C), but the p2 values for the hybrid epithelial (he) phenotype increased, whereas that of hybrid mesenchymal (hm) state decreased with little change in those of e or m phenotypes (Figure 4C–E). However, at all levels of PSF expression in all circuits (from 10-fold down-expression to a 10-fold over-expression), the hybrids together (he, hm) have a more considerable value of p1 and p2 in comparison to the “fully” epithelial and mesenchymal states (e, m) (Figure 4C–E and Figure S6A–C), suggesting that hybrid E/M phenotypes are more likely to be stem-like than either extreme of the EMT axis. Further, over-expression of PSFs may pull the composition of the solutions in the “stemness window” towards the epithelial end of the EMT axis.

We additionally looked at the distribution profile of the PSFs in all the clusters (Figure S7A–I) for all the circuits with down- and over-expression. Overall, the distribution profiles are distinct for all the circuits, perhaps indicating the role of specific topological features in dictating this distribution profile. GRHL2 remains appreciably bimodal throughout the range of expression, and the e and he clusters occupy the higher end of the GRHL2 histogram, while the m and hm clusters occupy the lower end (Figure S7A–C). OVOL is bimodal at 10-fold down-expression, but the profile shifts towards unimodality upon increasing expression, possibly due to its self-inhibition. Generally, a lower OVOL level is associated with the m and hm clusters and vice versa (Figure S7D–F). On the contrary, the NRF2 distribution profile does not show any readily appreciable preference for any of the clusters at any level of expression (Figure S7G–I).

The preference of GRHL2 and OVOL for the e and he clusters may be attributed to the presence of a mutually inhibitory negative feedback loop with ZEB, which causes the PSF (GRHL2, OVOL) levels to be relatively high when the ZEB levels are low (e and he) and vice versa. The GRHL2 and OVOL coupling with the EMP circuit has two differences: (a) GRHL2 activates miR-200 but OVOL does not, and (b) GRHL2 self-activates but OVOL self-inhibits. To deconvolute the effect of these two differences, we investigated a circuit where the activation link from GRHL2 to miR-200 is knocked down (GRHL2-KD). This circuit differs from the OVOL circuit only in terms of the nature of the self-regulation. We then investigated this knockdown circuit with regard to all the investigations explored for each of the other circuits above and compared them to the trends observed for the OVOL circuit. Qualitatively, all the trends matched well, and the only remarkable difference was in the cluster-wise distribution profile of the PSF at different levels of expression (Figure 4F and Figure S8A–F), suggesting that the specific enrichments observed are primarily dictated by the interacting links between the PSF node and the other nodes, while the PSF distribution profile itself is strongly affected by the nature of the self-link.

3. Discussion

Metastasis has been indicated as a hallmark of cancer [57]. However, decades of research have highlighted that hallmarks of metastasis are quite different from those of primary tumor growth [58]. Given the ever-changing biochemical and biomechanical environments for the cell in the metastatic cascade, it is not surprising that metastasis is a highly inefficient process and that only a minuscule percentage of disseminated cells that can pass through these bottlenecks by adapting to their dynamic environments can eventually colonize distant organs and form macrometastases. Thus, the concept of Lamarckian evolution seems to be more prevalent than that of Darwinian selection (of a subset of genetic clones) in the context of metastasis. Consequently, environment-driven phenotypic plasticity has been advocated as a hallmark of metastasis [58]. Recent observations in the evolution of cancer therapy resistance endorse the role of Lamarckian induction in enabling the dynamic adaptability of cancer cells [59,60,61,62,63].

Here, we elucidate the design principles of regulatory networks connecting two crucial aspects of phenotypic plasticity during metastasis—epithelial–mesenchymal plasticity (EMP) and stemness. Our results demonstrate that a higher likelihood of hybrid epithelial/mesenchymal (E/M) phenotype(s) acquiring increased stemness is an outcome of the network topology of a coupled EMP–stemness network. These findings offer a mechanistic basis of various recent in vitro, in vivo, and clinical work suggesting an enhanced stemness and/or aggressiveness trait of hybrid E/M phenotype(s), as reported across different carcinomas. A recent quantitative analysis of about a hundred urothelial carcinomas concluded that “it is not the amount, but merely the presence of a minimum of tumor cells in hybrid E/M states that contribute to disease aggressiveness” [64]. This clinical observation is reminiscent of in vivo experiments showing that both “extremely epithelial” (generated via CRISPR/Cas9 by eliminating Zeb1 expression) and “extremely mesenchymal” cells (generated by constitutive over-expression of Zeb1), either alone or in combination, had minimal tumor initiation potential as compared to hybrid E/M cells. Put together, these results highlight the importance of individual cells co-expressing epithelial and mesenchymal markers (i.e., one or more hybrid E/M phenotype(s)) in enabling stemness. At least two factors can underlie this enhanced stemness of hybrid E/M cells: (a) their ability to self-renew, as some cells undergoing a full EMT can undergo cell cycle arrest [65,66], and (b) their ability to give rise to heterogeneous subpopulations of epithelial and mesenchymal cells [19,20] that can potentially cooperate in driving metastatic traits [67]. Thus, hybrid E/M cells can be conceptually thought of as similar to adult stem cells in a tissue, a notion supported by a hybrid E/M phenotype seen in a subset of mammary stem cells [68,69].

It should be noted that we are proposing an enrichment of hybrid E/M phenotype(s) in stem-like populations, but not an exclusive or exhaustive association between hybrid E/M cells and stemness. Not every cell undergoing EMP, even a partial EMT, is expected to be stem-like; nor is every cancer stem cell (CSC) expected to manifest the hybrid E/M phenotype(s). CSCs of varying EMP status have been seen in breast cancer [11], prostate cancer [70], cervical cancer [71], and squamous cell carcinoma [72,73], with different possible spatial localization within a tumor [74]. This possibility can be explained by our model predictions demonstrating the presence of more epithelial or more mesenchymal phenotypes within the dynamic “stemness window”, albeit at a relatively lower frequency as compared to that of hybrid E/M phenotype(s). Increased stemness may lead to a higher tumor-initiating potential.

Our results show that stemness as a function of EMP is not necessarily always monotonically increasing. Instead, stemness is likely to first increase as cells begin to go through EMT and then decrease as they cross the range of hybrid phenotype(s) to become “fully” mesenchymal. Whether the terminal value of stemness is higher or lower than the initial “fully” epithelial remains to be rigorously quantified. Similar observations have been reported for tumor aggressiveness as a function of chromosomal instability (CIN) [75,76]. Such reports caution us against falling prey to tacit assumptions about monotonic relationships between any two molecular factors and/or phenomena and reveal the fundamental insights that can be gained by mechanistic modeling of nonlinear dynamics embedded in a complex regulatory network.

To investigate the robustness of our results and quantify the influence of various PSFs on the association between hybrid E/M phenotypes and stemness, we incorporated GRHL2, OVOL2, and NRF2 in our coupled network. We observed that these PSFs maintained the enrichment of hybrid E/M phenotypes in stem-like traits, endorsing the robustness of this association. Consistently, many of these PSFs have been reported to play a functional role in maintaining more epithelial phenotypes (e, he) as well as stemness in experiments in cancer cell lines and in cellular reprogramming contexts involving the modulation of EMP [53,77,78]. GRHL2 and OVOL2 have also been proposed as MET-inducing transcription factors (MET-TFs), but our simulations suggest that their ability to induce a “complete” MET need not be universal. This prediction is validated by recent experimental observations that the ability of GRHL2 to induce a MET depended on the epigenetic background of cells [54,79]. Thus, our results propose that. similar to variations in the ability of EMT-TFs to drive an EMT [80,81], MET-TFs may induce MET to varying degrees. Moreover, the epigenetic background of cells may alter the dynamics of EMT/MET and/or CSCs/non-CSCs and the interconnection between these two axes.

The approach presented here can be expanded to answer a long-standing open question in metastasis—how do cells coordinate various axes of their plasticity during the metastatic cascade? Mechanism-based models of individual axes of plasticity are the building blocks required to answer this question from a dynamical systems perspective. There has been increasing interest in dynamical modeling of EMP [82,83,84,85,86]; similar efforts to model related aspects of EMP, such as metabolic reprogramming [87], autophagy [88,89], therapy resistance [90,91], and immune evasion [92], are being attempted too. Coupling these modules to decode the interconnections among different axes of plasticity can help unravel the survival strategies of metastasis-initiating cells and eventually contribute to designing new clinical strategies. For instance, a recent prediction we made is that breaking the positive feedback loops in a network can restrict plasticity and potentially impact metastatic potential [6]. Indeed, proof-of-principle validation for this prediction was seen when cancer cells with a compromised ZEB/miR-200 feedback loop showed significant metastatic ability in vivo [93]. Future efforts are needed to dissect the functional consequences of breaking feedback loops among different axes of plasticity (such as the ones investigated here) in eventually delaying or preventing metastasis.

4. Methods

4.1. Random Circuit Perturbation (RACIPE)

RACIPE [43] is a tool allowing a parameter-agnostic computational simulation of regulatory networks based solely on the network topology. The network is composed of nodes representing various gene products and directed links representing inhibitory and activating links between them. For each node i with its expression level at a time t as $p_i$, with incoming links from a set $\left\{A\right\}$ of nodes, we thus have the ordinary differential equation for the dynamics of $p_i$:

(1)$\frac{d{p}_i}{dt}=g{\displaystyle \prod }_{j\in \left\{A\right\}}{H}^S\left(p_j,{\mu }_{ji},{\lambda }_{ji},n_{ji}\right)-k{p}_i$

where g is the basal production rate, k is the degradation rate, and $H^S$ is the shifted Hill function with parameters n, $\mu$, $\lambda$, p representing the effect of regulation by node j on i, altering its production rate. $\lambda >1$ represents an activating link and $\lambda <1$ represents an inhibitory link, with $\lambda =1$ representing no regulation. p is the threshold parameter that controls the level of p around which the non-linearity in $H^S$ is seen. n is the Hill coefficient representing the “extent of non-linearity” in the function by modulating the steepness of the response curve. See Supplementary Figure S9 for graphs of the shifted Hill function for varied parameters and Supplementary Section S3 for the explicit form of all the model equations. The final form of the shifted Hill function is given below.

(2)$H^S\left(p,\mu ,\lambda ,n\right)=\frac{{\mu }^n}{{\mu }^n+p^n}+\lambda \frac{p^n}{{\mu }^n+p^n}$

4.2. Parameter Sampling

Given a network topology, RACIPE generates an ensemble of random models, each with a randomly sampled set of parameters for every link. The default range of random sampling has been estimated from BioNumbers [94], and here we choose the default option of uniform random sampling from the range. Instead of directly randomizing the basal production rate parameter g, RACIPE randomizes the maximal production rate G, following

(3)$G=g{\displaystyle \prod }_{j\in \left\{A^p\right\}}{\lambda }_j$

where $\left\{A^p\right\}$ is the set of all incoming activating links (production rate when all activating links are maximally active and all inhibitory links are inactive). For each parameter set, RACIPE randomizes the initial values between the smallest possible steady state (all activating regulations are inactive and all inhibitory regulations are maximally active) and the largest possible state (the opposite of the previous configuration) in log scale.

RACIPE allows perturbing the circuit by changing the expression levels of different nodes relative to the other unperturbed nodes. Over-expression and downregulation by x-fold changes the sampling range for G by x times, hence, a 10-fold over-expression of a node will sample the G of the node from a range 10 times the range of the unperturbed (normal) node, while a 10-fold down-expression shrinks the range to 1/10th of its normal range. Details of the default ranges and the ranges used in our simulation are in Supplementary Section S1.1.2.

4.3. Simulation

We generate an ensemble of 10,000 models with each run (analogous to a biological replicate, and referred to elsewhere as such), and repeat this for 5 replicates to get the standard deviation for proportion (or frequency) calculations. RACIPE, for each of the models, simulates the network starting from multiple random initial positions until the network reaches a steady state, and prints as output the levels of each node at all identified steady-state solutions. The details of the run parameters used are in Supplementary Sections S1.1.1 and S1.1.2. The output of RACIPE for levels of all nodes are in log2 scale, and we perform a gene-wise normalization of the values before further analysis. For $p_i$, representing the actual steady-state levels of node i, the final values $z_i$ are obtained after normalizing the log-scale node levels:

(4)$z_i=\frac{lo{g}_2\left(p_i\right)-\overline{lo{g}_2\left(p_{ref}\right)}}{\sigma \left(lo{g}_2\left(p_{ref}\right)\right)}$

Here, $\overline{x}$ and $\sigma \left(x\right)$ represent the mean value and the standard deviation of the respective quantities. $p_{ref}$ is taken as the corresponding node in the circuit with no over-expressed or down-expressed nodes (the “reference circuit”), allowing a uniform normalization for circuits across perturbations (see Supplementary Section S1.2.3).

4.4. Link Strength Metrics

We define a link strength metric similar to the one used earlier [95] to quantify the strength of each link independent of its activation at steady state. The link strength ($l_{ij}$) is directly proportional to an activating ${\lambda }^p$ or inversely proportional to an inhibitory ${\lambda }^m$ and it is inversely proportional to the threshold $\mu$. To normalize the threshold, we choose the normalization value G/k, the maximum possible node level at steady state.

(5)$l_{ij}=\frac{G_i}{{\lambda }_{ij}^m{\mu }_{ij}{k}_i}=\frac{G_i{\lambda }_{ij}^p}{{\mu }_{ij}{k}_i}$

For any two links $l_{ij}$ and $l_{gh}$ (usually a part of a mutually inhibitory feedback motif between two nodes), the asymmetry is defined as $lo{g}_2\left(l_{ij}/l_{gh}\right)$, representing the overall asymmetry in the effective strengths of the two links forming a mutually inhibitory feedback loop. When quantifying the coupling strength, the combined strength is defined as $lo{g}_2\left(l_{ij}\right)+lo{g}_2\left(l_{gh}\right)$.

4.5. Calculation of p1 and p2 Values

For the calculation of the conditional probability of having a feature x conditional on having a feature y, we take the proportion of all steady-state solutions having the feature y that have the feature x too. We use two such conditional probabilities in our results. p1 is the probability of lying in the stemness window conditional on the solution belonging to a particular cluster, and is calculated as the proportion of all solutions of the specific cluster lying in the stemness window. p2 is the probability of belonging to a particular cluster conditional on the solution lying in the stemness window, and is calculated as the proportion of all solutions in the stemness window belonging to a particular cluster.

4.6. Clustering

For clustering, a K-means algorithm was used on the appropriately normalized data for only the four nodes (miR-200, ZEB, LIN28, and let7), with all solutions taken together with the same weight, without distinguishing between solutions belonging to systems with different numbers of stable states. K-means finds an optimal splitting of the data into K clusters (where K is a hyperparameter that the algorithm takes as the input) by iteratively minimizing inertia, or the total within-cluster sum of squares distance from the cluster centers. K-means results have a stochastic component since they depend upon the initialization of the cluster centers which are iteratively improved by the algorithm. To offset this, at each run, we randomize the initialization multiple times and select the best outcome in terms of minimized K-means inertia (see Supplementary Section S1.2.1 for further details of the algorithm parameters). In order to select a cluster number, we run the K-means algorithm with different values for K and at each value, compute multiple cluster quality metrics, all of which attempt to quantify the optimality of the clusters obtained at each value for K, and we then plot these metrics across different values of K. A sharp change in the slope as well as an overall higher value for the average silhouette widths indicate an appropriate partitioning [96]. For our data, the peak and the sudden change in slope [97] were almost unambiguously in the favor of K = 4 in all the metrics that we used (average silhouette widths, Calinski–Harabasz index [98], Davies–Bouldin index [99], K-means inertia) (see Supplementary Section S1.2.2 and Figure S2). Other relatively subjective supportive evidence can be gained from the visual evaluation of the PCA plot and observing the dendrogram computed using agglomerative hierarchical clustering using Ward’s minimum variance criterion [100]. The clustering algorithms and the cluster quality metrics were from scikit-learn version 0.32 (see Supplementary Section S1.1.3 for more details on the software versions and the platform). Each replicate was separately clustered and normalized (see Supplementary Section S1.2.3 for details on normalization).

This phenotype assignment strategy was maintained to be consistent all throughout the study for all circuits and replicates—a separate run of the K-means algorithm was performed for each replicate and the resultant clusters were assigned phenotypes based on their median ZEB levels (also see Supplementary Section S2.1).

4.7. Significance Testing

To compute the statistical significance for all comparisons, we used two strategies. For comparisons within the same RACIPE replicate (where n is large enough), we use the two-sided non-parametric Mann–Whitney U test [101] (also known as the Wilcoxon rank-sum test). For cross-replicate comparisons (i.e., proportions), we use Welch’s t-test [102], allowing unequal variances though with the assumption of normality of the compared groups. The Mann–Whitney U test is not suitable for small samples as it discards all information about the size of the differences (and as a result, caps the highest significance that can be achieved based on n) and also the assumption of normality of the proportions across replicates is a relatively reasonable assumption since it is stochastic differences which primarily cause the inter-replicate values to differ. Furthermore, we employ the Holm–Bonferroni correction to adjust the p-values accounting for multiple comparisons to control a family-wise error rate (FWER) in the overall rejection of the null hypotheses (see Supplementary Section S2.5). Due to the large number of steady states per RACIPE replicate, for some comparisons within the same replicate, very small differences of means (or medians) give a statistically significant p-value but which is not necessarily biologically relevant. Thus, to ascribe possible biological meaning to such comparisons, for cases of an irrelevantly small “magnitude” of this difference (i.e., the effect size), we provide the ratio of the means/medians and the absolute difference of these quantities.

Author Contributions

M.K.J. conceived and supervised research; S.P. performed research; S.P. and S.S. analyzed data. All authors contributed to writing and editing of the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by a Ramanujan Fellowship awarded to M.K.J. by the Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India (SB/S2/RJN-049/2018).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

All relevant codes are available at https://github.com/Stochastic13/emt-stemness.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Four phenotypes enabled by the coupled Epithelial Mesenchymal Plasticity (EMP)–stemness circuit. (A) The core EMP–stemness circuit (base circuit). Arrows represent activating links and hammerheads represent inhibitory links. (B) Heatmap showing the steady-state expression of all solutions of a single replicate including the core EMP (ZEB, miR-200) and stemness (LIN28, let-7) nodes. (C) Expression histograms (all solutions from a single Random Circuit Perturbation (RACIPE) replicate. For n = 5, see Supplementary Table S2 in Section S2.1) showing ZEB levels in the four clusters—epithelial (e), hybrid epithelial (he), hybrid mesenchymal (hm), and mesenchymal (m). The colors of the labels shown in the heatmap are for e, he, hm, and m clusters displayed here. (D) Expression levels of EMP and the stemness genes in the four clusters. In each node, the horizontal line represents the median value, the darkest region represents the middle 30% (35th to 65th percentile) of the data, and the region next in lightness represents the next 40% (15th to 35th and 65th to 85th percentiles), while the lightest region represents the next 20% of the data (5th to 15th and 85th to 95th percentiles) (E) Principal Component Analysis (PCA) plot colored according to the K-means cluster assignment. (Principal Component 1 (PC1) = 0.48 × miR200 − 0.54 × ZEB − 0.21 × SNAIL − 0.49 × LIN28 + 0.44 × let7 + 0.057 × NF-kB; Principal Component 2 (PC2) = 0.27 × miR200 − 0.08 × ZEB − 0.66 × SNAIL + 0.14 × LIN28 − 0.49 × let7 − 0.47 × NF-kB). For (C,D): Significance bars added to show statistically significant differences (Mann–Whitney U test). Legend for the p-values: *** p < 0.0001; for effect sizes, refer to Table S7 in Supplementary Section S2.5.

Figure 2. Hybrid E/M phenotypes are enriched for stemness. (A) Proportion of all monostable parameter sets belonging to different phenotypes—e, he, hm, m. (B) Expression histogram for LIN28 (all solutions of a single RACIPE replicate) for the four phenotypes with the stemness window (black vertical lines) as defined by the middle 30% of the “LIN28 biological range” (median ± one interquartile range averaged over n = 5 RACIPE replicates). Significance bars added to show statistically significant differences (Mann–Whitney U test) (C) Value of p1—the probability of a solution lying in the stemness window, conditional on the solution belonging to a particular phenotype for all the phenotypes (n = 5 RACIPE replicates). (D) Same as (C) except for p2—the probability of a solution belonging to a particular phenotype conditional on it lying in the stemness window. For (A,C,D): Mean values show average across five RACIPE replicates and the error bars represent the across-replicate standard deviations. Two-tailed Student’s t-test with unequal variance (Welch’s t-test) was performed. Legend for the p-values: + p > 0.01, * 0.01 > p > 0.001, ** 0.001 > p > 0.0001, *** p < 0.0001; for effect sizes, refer to Table S7 in Supplementary Section S2.5.

Figure 3. Including phenotypic stability factors (PSFs) preserves the cluster structure seen for the EMP–stemness circuit. (A) Different PSFs (GRHL2, OVOL2, NRF2) added to the base (control) EMP–stemness circuit. All new nodes are shown in pink; the PSF node added is marked with a star. (B) Multistability profiles for all the PSF circuits compared to the base circuit. (C) The average silhouette widths at different numbers of clusters for all PSF circuits with the cluster membership assignment via a K-means algorithm. Note the distinct peak at four clusters across all circuits. (D) Cluster centers in the ZEB–LIN28 plane (EMP–stemness plane) for a single RACIPE replicate of each of the circuits show that the cluster centers are well preserved across the circuits. Gene-wise normalization is performed independently for all circuits (see Supplementary Section S1.2.3). (E) Phenotypic distributions for all monostable solutions upon over-expressing and down-expressing the different PSF nodes. (de10: 10-fold down-expression, ref: Reference unperturbed circuit, oe10: 10-fold over-expression) ((B,C,E): The values are averaged across n = 5 RACIPE replicates and the error bars represent the across-replicate standard deviations).

Figure 4. PSF motifs maintain the association between hybrid E/M states and stemness. (A) Sum of p1 values for the two hybrid (he, hm) and two non-hybrid (e,m) phenotypes, for different circuits including PSFs. (B) Same as (A) but for p2. (C–E) p2 values for four phenotypes for cases of over-expression and downregulation of each of the PSFs individually: NRF2 (C), GRHL2 (D), and OVOL (E). (F) Stacked expression histograms (colored in accordance with the phenotypes) showing the change in the PSF node distribution in the GRHL2-KD circuit (GRHL2-KD) (left panel) vs. the OVOL circuit (right panel) at 10-fold over-expression. GRHL2-KD involves knockdown of the activating link from GRHL2 to miR-200, causing the only difference between the two presented circuits to be the nature of the self-link of the PSF node. All other links remain the same. (de10: 10-fold down-expression, ref: Reference unperturbed circuit, oe10: 10-fold overexpression) (A–E): The values are averaged across replicates and the error bars represent the across-replicate standard deviations; refer to Table S8 in Supplementary Section S2.5 for the statistical testing of these quantities.

Supplementary Materials

The following are available online at https://www.mdpi.com/2077-0383/10/1/60/s1, Figure S1: Optimal cluster number by varied cluster quality, Figure S2: Four E/M phenotypes enabled by the base circuit, Figure S3: Stemness window and the bistable/tristable phase distributions of the coupled EMP–stemness circuit, Figure S4: Link strength analysis, Figure S5: Distribution of bistable and tristable systems upon incorporating PSFs with the coupled EMP–stemness circuit, Figure S6: Stemness probability (p1) changes with PSF expression, Figure S7: PSF expression profiles at different expression levels, Figure S8: GRHL2-KD circuit results for the investigated quantities, Figure S9: Shifted Hill function across different parameters, Table S1: Parameters ranges and the simulation options, Table S2: Median (iqr) for ZEB and LIN28 in different replicates of the base circuit rounded off to 3 significant figures, Table S3: Average Silhouette widths for the Uncoupled E/M and Stemness circuits, Table S4: Proportion (averaged across all replicates) of filtered x-stable parameter sets based on > 1 solutions of the set being assigned the same phenotype (standard deviations in the parentheses), Table S5: Proportion (averaged across all replicates) of all solutions irrespective of the number of stable states of the corresponding parameter set belonging to different phenotypes (with standard deviations in parenthesis), Table S6: Statistical Testing results for differences of various calculated proportions in the base circuit considering all the replicates (i.e., the proportion calculated for the 5 replicates forms one comparison group), Table S7: Statistical Testing results for differences of the node distributions across the clusters for a single replicate of the base circuit, Table S8: Statistical Testing results for the PSF circuits with the proportions calculated for all the replicates of the de10 circuit forming one comparison group while the proportions of all the replicates of the oe10 circuit forming the other comparison group, Table S9: p1 and p2 for solutions of multi-stable systems (after removing monostable systems), Table S10: Different configurations of the association of stemness with the different phenotypes in bistable phases for base circuit, Table S11: Different configurations of the association of stemness with the different phenotypes in tristable phases for base circuit, Table S12: Different configurations of the association of stemness with the different phenotypes in bistable phases for GRHL2 circuit (reference un-perturbed), Table S13: Different configurations of the association of stemness with the different phenotypes in tristable phases for GRHL2 circuit(reference un-perturbed), Table S14: Different configurations of the association of stemness with the different phenotypes in bistable phases for OVOL circuit(reference un-perturbed), Table S15: Different configurations of the association of stemness with the different phenotypes in tristable phases for OVOL circuit(reference un-perturbed), Table S16: Different configurations of the association of stemness with the different phenotypes in bistable phases for NRF2 circuit(reference un-perturbed), Table S17: Different configurations of the association of stemness with the different phenotypes in tristable phases for NRF2 circuit(reference un-perturbed).

References

1. Lambert, A.W.; Pattabiraman, D.R.; Weinberg, R.A. Emerging Biological Principles of Metastasis. Cell; 2017; 168, pp. 670-691. [DOI: https://dx.doi.org/10.1016/j.cell.2016.11.037]

2. Fidler, I.J. Metastasis: Quantitative Analysis of Distribution and Fate of Tumor Emboli Labeled. J. Natl. Cancer Inst.; 1970; 45, pp. 773-782. [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/5513503]

3. Celià-Terrassa, T.; Jolly, M.K. Cancer Stem Cells and Epithelial-to-Mesenchymal Transition in Cancer Metastasis. Cold Spring Harb. Perspect. Med.; 2020; 10, a036905. [DOI: https://dx.doi.org/10.1101/cshperspect.a036905] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31570380]

4. Celià-Terrassa, T.; Kang, Y. Distinctive properties of metastasis-initiating cells. Genes Dev.; 2016; 30, pp. 892-908. [DOI: https://dx.doi.org/10.1101/gad.277681.116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27083997]

5. Jia, D.; Li, X.; Bocci, F.; Tripathi, S.; Deng, Y.; Jolly, M.K.; Onuchic, J.N.; Levine, H. Quantifying Cancer Epithelial-Mesenchymal Plasticity and its Association with Stemness and Immune Response. J. Clin. Med.; 2019; 8, 725. [DOI: https://dx.doi.org/10.3390/jcm8050725]

6. Hari, K.; Sabuwala, B.; Subramani, B.V.; La Porta, C.; Zapperi, S.; Font-Clos, F.; Jolly, M.K. Identifying inhibitors of epithelial-mesenchymal plasticity using a network topology based approach. NPJ Syst. Biol. Appl.; 2020; 6, 15. [DOI: https://dx.doi.org/10.1038/s41540-020-0132-1]

7. Williams, E.D.; Gao, D.; Redfern, A.; Thompson, E.W. Controversies around epithelial–mesenchymal plasticity in cancer metastasis. Nat. Rev. Cancer; 2019; 19, pp. 716-732. [DOI: https://dx.doi.org/10.1038/s41568-019-0213-x]

8. Mani, S.A.; Guo, W.; Liao, M.-J.; Eaton, E.N.; Ayyanan, A.; Zhou, A.Y.; Brooks, M.; Reinhard, F.; Zhang, C.C.; Shipitsin, M. et al. The epithelial-mesenchymal transition generates cells with properties of stem cells. Cell; 2008; 133, pp. 704-715. [DOI: https://dx.doi.org/10.1016/j.cell.2008.03.027]

9. Morel, A.-P.; Lièvre, M.; Thomas, C.; Hinkal, G.; Ansieau, S.; Puisieux, A. Generation of breast cancer stem cells through epithelial-mesenchymal transition. PLoS ONE; 2008; 3, e2888. [DOI: https://dx.doi.org/10.1371/journal.pone.0002888]

10. Celià-Terrassa, T.; Meca-Cortés, Ó.; Mateo, F.; De Paz, A.M.; Rubio, N.; Arnal-Estapé, A.; Ell, B.J.; Bermudo, R.; Díaz, A.; Guerra-Rebollo, M. et al. Epithelial-mesenchymal transition can suppress major attributes of human epithelial tumor-initiating cells. J. Clin. Investig.; 2012; 122, pp. 1849-1868. [DOI: https://dx.doi.org/10.1172/JCI59218]

11. Liu, S.; Cong, Y.; Wang, D.; Sun, Y.; Deng, L.; Liu, Y.; Martin-Trevino, R.; Shang, L.; McDermott, S.P.; Landis, M.D. et al. Breast cancer stem cells transition between epithelial and mesenchymal states reflective of their normal counterparts. Stem Cell Rep.; 2014; 2, pp. 78-91. [DOI: https://dx.doi.org/10.1016/j.stemcr.2013.11.009] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24511467]

12. Grosse-Wilde, A.; Fouquier d’ Herouei, A.; McIntosh, E.; Ertaylan, G.; Skupin, A.; Kuestner, R.E.; del Sol, A.; Walters, K.-A.; Huang, S. Stemness of the hybrid epithelial/mesenchymal state in breast cancer and its association with poor survival. PLoS ONE; 2015; 10, e0126522. [DOI: https://dx.doi.org/10.1371/journal.pone.0126522]

13. Kröger, C.; Afeyan, A.; Mraz, J.; Ng, E.; Reinhardt, F.; Khodor, Y.L. Acquisition of a hybrid E / M state is essential for tumorigenicity of basal breast cancer cells. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 7353-7362. [DOI: https://dx.doi.org/10.1073/pnas.1812876116] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30910979]

14. Bierie, B.; Pierce, S.E.; Kroeger, C.; Stover, D.G.; Pattabiraman, D.R.; Thiru, P.; Liu Donaher, J.; Reinhardt, F.; Chaffer, C.L.; Keckesova, Z. et al. Integrin-β4 identifies cancer stem cell-enriched populations of partially mesenchymal carcinoma cells. Proc. Natl. Acad. Sci. USA; 2017; 114, pp. E2337-E2346. [DOI: https://dx.doi.org/10.1073/pnas.1618298114] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28270621]

15. Jolly, M.K.; Huang, B.; Lu, M.; Mani, S.A.; Levine, H.; Ben-Jacob, E. Towards elucidating the connection between epithelial-mesenchymal transitions and stemness. J. R. Soc. Interface; 2014; 11, 20140962. [DOI: https://dx.doi.org/10.1098/rsif.2014.0962] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25339690]

16. Pastushenko, I.; Brisebarre, A.; Sifrim, A.; Fioramonti, M.; Revenco, T.; Boumahdi, S.; Van Keymeulen, A.; Brown, D.; Moers, V.; Lemaire, S. et al. Identification of the tumour transition states occurring during EMT. Nature; 2018; 556, pp. 463-468. [DOI: https://dx.doi.org/10.1038/s41586-018-0040-3] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29670281]

17. Ombrato, L.; Malanchi, I. The EMT Universe: Space between Cancer Cell Dissemination and Metastasis Initiation. Crit. Rev. Oncog.; 2014; 19, pp. 349-361. [DOI: https://dx.doi.org/10.1615/CritRevOncog.2014011802]

18. Quan, Q.; Wang, X.; Lu, C.; Ma, W.; Wang, Y.; Xia, G.; Wang, C.; Yang, G. Cancer stem-like cells with hybrid epithelial/mesenchymal phenotype leading the collective invasion. Cancer Sci.; 2020; 111, pp. 467-476. [DOI: https://dx.doi.org/10.1111/cas.14285]

19. Ruscetti, M.; Dadashian, E.L.; Guo, W.; Quach, B.; Mulholland, D.J.; Park, J.W.; Tran, L.M.; Kobayashi, N.; Bianchi-Frias, D.; Xing, Y. et al. HDAC inhibition impedes epithelial-mesenchymal plasticity and suppresses metastatic, castration-resistant prostate cancer. Oncogene; 2016; 35, pp. 3781-3795. [DOI: https://dx.doi.org/10.1038/onc.2015.444]

20. Tripathi, S.; Chakraborty, P.; Levine, H.; Jolly, M.K. A mechanism for epithelial-mesenchymal heterogeneity in a population of cancer cells. PLoS Comput. Biol.; 2020; 16, e1007619. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1007619]

21. Jolly, M.K.; Somarelli, J.A.; Sheth, M.; Biddle, A.; Tripathi, S.C.C.; Armstrong, A.J.J.; Hanash, S.M.M.; Bapat, S.A.A.; Rangarajan, A.; Levine, H. Hybrid epithelial/mesenchymal phenotypes promote metastasis and therapy resistance across carcinomas. Pharmacol. Ther.; 2019; 194, pp. 161-184. [DOI: https://dx.doi.org/10.1016/j.pharmthera.2018.09.007] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30268772]

22. Jolly, M.K.; Jia, D.; Boareto, M.; Mani, S.A.; Pienta, K.J.; Ben-Jacob, E.; Levine, H. Coupling the modules of EMT and stemness: A tunable “stemness window” model. Oncotarget; 2015; 6, pp. 25161-25174. [DOI: https://dx.doi.org/10.18632/oncotarget.4629] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26317796]

23. Bocci, F.; Jolly, M.K.; George, J.T.; Levine, H.; Onuchic, J.N. A mechanism-based computational model to capture the interconnections among epithelial-mesenchymal transition, cancer stem cells and Notch-Jagged signaling. Oncotarget; 2018; 9, pp. 29906-29920. [DOI: https://dx.doi.org/10.18632/oncotarget.25692] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30042822]

24. Bracken, C.P.; Gregory, P.A.; Kolesnikoff, N.; Bert, A.G.; Wang, J.; Shannon, M.F.; Goodall, G.J. A double-negative feedback loop between ZEB1-SIP1 and the microRNA-200 family regulates epithelial-mesenchymal transition. Cancer Res.; 2008; 68, pp. 7846-7854. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-08-1942]

25. Park, S.-M.M.; Gaur, A.B.; Lengyel, E.; Peter, M.E. The miR-200 family determines the epithelial phenotype of cancer cells by targeting the E-cadherin repressors ZEB1 and ZEB2. Genes Dev.; 2008; 22, pp. 894-907. [DOI: https://dx.doi.org/10.1101/gad.1640608]

26. Preca, B.-T.; Bajdak, K.; Mock, K.; Lehmann, W.; Sundararajan, V.; Bronsert, P.; Matzge-Ogi, A.; Orian-Rousseau, V.; Brabletz, S.; Brabletz, T. et al. A novel ZEB1/HAS2 positive feedback loop promotes EMT in breast cancer. Oncotarget; 2017; 8, pp. 11530-11543. [DOI: https://dx.doi.org/10.18632/oncotarget.14563]

27. Jolly, M.K.; Preca, B.-T.; Tripathi, S.C.; Jia, D.; George, J.T.; Hanash, S.M.; Brabletz, T.; Stemmler, M.P.; Maurer, J.; Levine, H. Interconnected feedback loops among ESRP1, HAS2, and CD44 regulate epithelial-mesenchymal plasticity in cancer. APL Bioeng.; 2018; 2, 031908. [DOI: https://dx.doi.org/10.1063/1.5024874]

28. Preca, B.-T.; Bajdak, K.; Mock, K.; Sundararajan, V.; Pfannstiel, J.; Maurer, J.; Wellner, U.; Hopt, U.T.; Brummer, T.; Brabletz, S. et al. A self-reinforcing CD44s/ZEB1 feedback loop maintains EMT and stemness properties in cancer cells. Int. J. Cancer; 2015; 137, pp. 2566-2577. [DOI: https://dx.doi.org/10.1002/ijc.29642]

29. Lu, M.; Jolly, M.K.; Levine, H.; Onuchic, J.N.; Ben-Jacob, E. MicroRNA-based regulation of epithelial-hybrid-mesenchymal fate determination. Proc. Natl. Acad. Sci. USA; 2013; 110, pp. 18144-18149. [DOI: https://dx.doi.org/10.1073/pnas.1318192110]

30. Yang, X.; Lin, X.; Zhong, X.; Kaur, S.; Li, N.; Liang, S.; Lassus, H.; Wang, L.; Katsaros, D.; Montone, K. et al. Double-negative feedback loop between reprogramming factor LIN28 and microRNA let-7 regulates aldehyde dehydrogenase 1-positive cancer stem cells. Cancer Res.; 2010; 70, pp. 9463-9472. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-10-2388]

31. Hafner, M.; Max, K.E.A.; Bandaru, P.; Morozov, P.; Gerstberger, S.; Brown, M.; Molina, H.; Tuschl, T. Identification of mRNAs bound and regulated by human LIN28 proteins and molecular requirements for RNA recognition. RNA; 2013; 19, pp. 613-626. [DOI: https://dx.doi.org/10.1261/rna.036491.112] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23481595]

32. Wilbert, M.L.; Huelga, S.C.; Kapeli, K.; Stark, T.J.; Liang, T.Y.; Chen, S.X.; Yan, B.Y.; Nathanson, J.L.; Hutt, K.R.; Lovci, M.T. et al. LIN28 binds messenger RNAs at GGAGA motifs and regulates splicing factor abundance. Mol. Cell; 2012; 48, pp. 195-206. [DOI: https://dx.doi.org/10.1016/j.molcel.2012.08.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22959275]

33. Zisoulis, D.G.; Kai, Z.S.; Chang, R.K.; Pasquinelli, A.E. Autoregulation of microRNA biogenesis by let-7 and Argonaute. Nature; 2012; 486, pp. 541-544. [DOI: https://dx.doi.org/10.1038/nature11134] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22722835]

34. Qiu, C.; Ma, Y.; Wang, J.; Peng, S.; Huang, Y. Lin28-mediated post-transcriptional regulation of Oct4 expression in human embryonic stem cells. Nucleic Acids Res.; 2010; 38, pp. 1240-1248. [DOI: https://dx.doi.org/10.1093/nar/gkp1071]

35. Niwa, H.; Miyazaki, J.; Smith, A.G. Quantitative expression of Oct-3/4 defines differentiation, dedifferentiation or self-renewal of ES cells. Nat. Genet.; 2000; 24, pp. 372-376. [DOI: https://dx.doi.org/10.1038/74199]

36. Karwacki-Neisius, V.; Göke, J.; Osorno, R.; Halbritter, F.; Ng, J.H.; Weiße, A.Y.; Wong, F.C.K.; Gagliardi, A.; Mullin, N.P.; Festuccia, N. et al. Reduced Oct4 expression directs a robust pluripotent state with distinct signaling activity and increased enhancer occupancy by Oct4 and Nanog. Cell Stem Cell; 2013; 12, pp. 531-545. [DOI: https://dx.doi.org/10.1016/j.stem.2013.04.023]

37. Theunissen, T.W.; Van Oosten, A.L.; Castelo-Branco, G.; Hall, J.; Smith, A.; Silva, J.C.R. Nanog overcomes reprogramming barriers and induces pluripotency in minimal conditions. Curr. Biol.; 2011; 21, pp. 65-71. [DOI: https://dx.doi.org/10.1016/j.cub.2010.11.074]

38. Sato, H.; Shien, K.; Tomida, S.; Okayasu, K.; Suzawa, K.; Hashida, S.; Torigoe, H.; Watanabe, M.; Yamamoto, H.; Soh, J. et al. Targeting the miR-200c/LIN28B axis in acquired EGFR-TKI resistance non-small cell lung cancer cells harboring EMT features. Sci. Rep.; 2017; 7, 40847. [DOI: https://dx.doi.org/10.1038/srep40847]

39. Kong, D.; Banerjee, S.; Ahmad, A.; Li, Y.; Wang, Z.; Sethi, S.; Sarkar, F.H. Epithelial to mesenchymal transition is mechanistically linked with stem cell signatures in prostate cancer cells. PLoS ONE; 2010; 5, e12445. [DOI: https://dx.doi.org/10.1371/journal.pone.0012445]

40. Liu, Y.; Li, H.; Feng, J.; Cui, X.; Huang, W.; Li, Y.; Su, F.; Liu, Q.; Zhu, J.; Lv, X. et al. Lin28 induces epithelial-to-mesenchymal transition and stemness via downregulation of let-7a in breast cancer cells. PLoS ONE; 2013; 8, e83083. [DOI: https://dx.doi.org/10.1371/journal.pone.0083083]

41. Brabletz, S.; Brabletz, T. The ZEB/miR-200 feedback loop—A motor of cellular plasticity in development and cancer?. EMBO Rep.; 2010; 11, pp. 670-677. [DOI: https://dx.doi.org/10.1038/embor.2010.117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20706219]

42. Dong, P.; Karaayvaz, M.; Jia, N.; Kaneuchi, M.; Hamada, J.; Watari, H.; Sudo, S.; Ju, J.; Sakuragi, N. Mutant p53 gain-of-function induces epithelial-mesenchymal transition through modulation of the miR-130b-ZEB1 axis. Oncogene; 2013; 32, pp. 3286-3295. [DOI: https://dx.doi.org/10.1038/onc.2012.334] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22847613]

43. Huang, B.; Lu, M.; Jia, D.; Ben-Jacob, E.; Levine, H.; Onuchic, J.N. Interrogating the topological robustness of gene regulatory circuits by randomization. PLoS Comput. Biol.; 2017; 13, e1005456. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1005456] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28362798]

44. Font-Clos, F.; Zapperi, S.; Porta, C.A.M. La Topography of epithelial–mesenchymal plasticity. Proc. Natl. Acad. Sci. USA; 2018; 115, pp. 5902-5907. [DOI: https://dx.doi.org/10.1073/pnas.1722609115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29784817]

45. Steinway, S.N.; Zañudo, J.G.T.; Michel, P.J.; Feith, D.J.; Loughran, T.P.; Albert, R. Combinatorial interventions inhibit TGFβ-driven epithelial-to-mesenchymal transition and support hybrid cellular phenotypes. NPJ Syst. Biol. Appl.; 2015; 1, 15014. [DOI: https://dx.doi.org/10.1038/npjsba.2015.14]

46. Jolly, M.K.; Tripathi, S.C.; Jia, D.; Mooney, S.M.; Celiktas, M.; Hanash, S.M.; Mani, S.A.; Pienta, K.J.; Ben-Jacob, E.; Levine, H. Stability of the hybrid epithelial/mesenchymal phenotype. Oncotarget; 2016; 7, pp. 27067-27084. [DOI: https://dx.doi.org/10.18632/oncotarget.8166]

47. Hong, T.; Watanabe, K.; Ta, C.H.; Villarreal-Ponce, A.; Nie, Q.; Dai, X. An Ovol2-Zeb1 Mutual Inhibitory Circuit Governs Bidirectional and Multi-step Transition between Epithelial and Mesenchymal States. PLoS Comput. Biol.; 2015; 11, e1004569. [DOI: https://dx.doi.org/10.1371/journal.pcbi.1004569]

48. Bocci, F.; Tripathi, S.C.; Vilchez, M.S.A.; George, J.T.; Casabar, J.; Wong, P.; Hanash, S.; Levine, H.; Onuchic, J.; Jolly, M. NRF2 activates a partial Epithelial-Mesenchymal Transition and is maximally present in a hybrid Epithelial/Mesenchymal phenotype. Integr. Biol.; 2019; 11, pp. 251-263. [DOI: https://dx.doi.org/10.1093/intbio/zyz021]

49. Watanabe, K.; Villarreal-Ponce, A.; Sun, P.; Salmans, M.L.; Fallahi, M.; Andersen, B.; Dai, X. Mammary morphogenesis and regeneration require the inhibition of EMT at terminal end buds by Ovol2 transcriptional repressor. Dev. Cell; 2014; 29, pp. 59-74. [DOI: https://dx.doi.org/10.1016/j.devcel.2014.03.006]

50. Somarelli, J.A.; Shetler, S.; Jolly, M.K.; Wang, X.; Bartholf Dewitt, S.; Hish, A.J.; Gilja, S.; Eward, W.C.; Ware, K.E.; Levine, H. et al. Mesenchymal-epithelial transition in sarcomas is controlled by the combinatorial expression of microRNA 200s and GRHL2. Mol. Cell. Biol.; 2016; 36, pp. 2503-2513. [DOI: https://dx.doi.org/10.1128/MCB.00373-16]

51. Cieply, B.; Riley, P., IV; Pifer, P.M.; Widmeyer, J.; Addison, J.B.; Ivanov, A.V.; Denvir, J.; Frisch, S.M. Suppression of the epithelial-mesenchymal transition by grainyhead-like-2. Cancer Res.; 2012; 72, pp. 2440-2453. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-11-4038] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/22379025]

52. Chung, V.Y.; Tan, T.Z.; Tan, M.; Wong, M.K.; Kuay, K.T.; Yang, Z.; Ye, J.; Muller, J.; Koh, C.M.; Guccione, E. et al. GRHL2-miR-200-ZEB1 maintains the epithelial status of ovarian cancer through transcriptional regulation and histone modification. Sci. Rep.; 2016; 6, 19943. [DOI: https://dx.doi.org/10.1038/srep19943] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26887977]

53. Saxena, K.; Srikrishnan, S.; Celia-Terrassa, T.; Jolly, M.K. OVOL1/2: Drivers of epithelial differentiation in development, disease and reprogramming. Cells Tissues Organs; 2020; in press

54. Chung, V.Y.; Tan, T.Z.; Ye, J.; Huang, R.-L.; Lai, H.-C.; Kappei, D.; Wollmann, H.; Guccione, E.; Huang, R.Y.-J. The role of GRHL2 and epigenetic remodeling in epithelial–mesenchymal plasticity in ovarian cancer cells. Commun. Biol.; 2019; 2, 272. [DOI: https://dx.doi.org/10.1038/s42003-019-0506-3]

55. Farris, J.C.; Pifer, P.M.; Zheng, L.; Gottlieb, E.; Denvir, J.; Frisch, S.M. Grainyhead-like 2 reverses the metabolic changes induced by the oncogenic epithelial-mesenchymal transition: Effects on anoikis. Mol. Cancer Res.; 2016; 14, pp. 528-538. [DOI: https://dx.doi.org/10.1158/1541-7786.MCR-16-0050]

56. MacFawn, I.; Wilson, H.; Selth, L.A.; Leighton, I.; Serebriiskii, I.; Bleackley, R.C.; Elzamzamy, O.; Farris, J.; Pifer, P.M.; Richer, J. et al. Grainyhead-like-2 confers NK-sensitivity through interactions with epigenetic modifiers. Mol. Immunol.; 2019; 105, pp. 137-149. [DOI: https://dx.doi.org/10.1016/j.molimm.2018.11.006]

57. Hanahan, D.; Weinberg, R.A. The Hallmarks of Cancer. Cell; 2000; 100, pp. 57-70. [DOI: https://dx.doi.org/10.1016/S0092-8674(00)81683-9]

58. Welch, D.R.; Hurst, D.R. Defining the Hallmarks of Metastasis. Cancer Res.; 2019; 79, pp. 3011-3027. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-19-0458]

59. Pisco, A.O.; Brock, A.; Zhou, J.; Moor, A.; Mojtahedi, M.; Jackson, D.; Huang, S. Non-Darwinian dynamics in therapy-induced cancer drug resistance. Nat. Commun.; 2013; 4, 2467. [DOI: https://dx.doi.org/10.1038/ncomms3467]

60. Salgia, R.; Kulkarni, P. The genetic/non-genetic duality of drug “resistance.”. Trends Cancer; 2018; [DOI: https://dx.doi.org/10.1016/j.trecan.2018.01.001]

61. Jolly, M.K.; Kulkarni, P.; Weninger, K.; Orban, J.; Levine, H. Phenotypic plasticity, bet-hedging, and androgen independence in prostate cancer: Role of non-genetic heterogeneity. Front. Oncol.; 2018; 8, 50. [DOI: https://dx.doi.org/10.3389/fonc.2018.00050] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29560343]

62. Robles-Espinoza, C.D. Beating tumour drug resistance: “Lamarckian” induction in the spotlight. Pigment Cell Melanoma Res.; 2019; 32, pp. 6-8. [DOI: https://dx.doi.org/10.1111/pcmr.12744] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30339326]

63. Marine, J.-C.; Dawson, S.-J.; Dawson, M.A. Non-genetic mechanisms of therapeutic resistance in cancer. Nat. Rev. Cancer; 2020; in press [DOI: https://dx.doi.org/10.1038/s41568-020-00302-4] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33033407]

64. Godin, L.; Balsat, C.; Van Eycke, Y.; Allard, J.; Royer, C.; Remmelink, M.; Pastushenko, I.; Haene, N.D.; Blanpain, C.; Salmon, I. et al. A Novel Approach for Quantifying Cancer Cells Showing Hybrid Epithelial / Mesenchymal States in Large Series of Tissue Samples: Towards a New Prognostic Marker. Cancers; 2020; 12, 906. [DOI: https://dx.doi.org/10.3390/cancers12040906] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32276404]

65. Vega, S.; Morales, A.V.; Ocaña, O.H.; Valdés, F.; Fabregat, I.; Nieto, M.A. Snail blocks the cell cycle and confers resistance to cell death. Genes Dev.; 2004; 18, pp. 1131-1143. [DOI: https://dx.doi.org/10.1101/gad.294104] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15155580]

66. Comaills, V.; Kabeche, L.; Morris, R.; Zou, L.; Daniel, A.; Yu, M.; Madden, M.W.; Licausi, J.A.; Boukhali, M.; Tajima, K. et al. Genomic Instability Is Induced by Persistent Proliferation of Cells Undergoing Epithelial-to- Mesenchymal Transition. Cell Rep.; 2016; 17, pp. 2632-2647. [DOI: https://dx.doi.org/10.1016/j.celrep.2016.11.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27926867]

67. Tsuji, T.; Ibaragi, S.; Shima, K.; Hu, M.G.; Katsurano, M.; Sasaki, A.; Hu, G.F. Epithelial-mesenchymal transition induced by growth suppressor p12 CDK2-AP1 promotes tumor cell local invasion but suppresses distant colony growth. Cancer Res.; 2008; 68, pp. 10377-10386. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-08-1444]

68. Colacino, J.A.; Azizi, E.; Brooks, M.D.; Harouaka, R.; Fouladdel, S.; McDermott, S.P.; Lee, M.; Hill, D.; Madde, J.; Boerner, J. et al. Heterogeneity of human breast stem and progenitor cells as revelaed by transcriptional profiling. Stem Cell Rep.; 2018; 10, pp. 1596-1609. [DOI: https://dx.doi.org/10.1016/j.stemcr.2018.03.001]

69. Tong, T.; Wang, Y.; Brooks, M.D.; Lee, C.; Scott, C.; Balzano, L.; Wicha, M.S.; Colacino, J.A. Hybrid stem cell states: Insights into the relationship between mammary development and breast cancer using single-cell transcriptomics. Front. Cell Dev. Biol.; 2020; in press [DOI: https://dx.doi.org/10.3389/fcell.2020.00288]

70. Galoczova, M.; Nenutil, R.; Pokorna, Z.; Vojtesek, B.; Coates, P. TAp63 and ΔNp63 (p40) in prostate adenocarcinomas: ΔNp63 associates with a basal-like cancer stem cell population but not with metastasis. Virchows Arch.; 2020; in press [DOI: https://dx.doi.org/10.1007/s00428-020-02944-z]

71. He, Y.; Xiao, M.; Fu, H.; Chen, L.; Qi, L.; Liu, D.; Guo, P.; Chen, L.; Luo, Y.; Xiao, H. et al. cPLA2α reversibly regulates different subsets of cancer stem cells transformation in cervical cancer. Stem Cells; 2020; 38, pp. 487-503. [DOI: https://dx.doi.org/10.1002/stem.3157] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32100928]

72. Biddle, A.; Liang, X.; Gammon, L.; Fazil, B.; Harper, L.J.; Emich, H.; Costea, D.E.; Mackenzie, I.C. Cancer stem cells in squamous cell carcinoma switch between two distinct phenotypes that are preferentially migratory or proliferative. Cancer Res.; 2011; 71, pp. 5317-5326. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-11-1059] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21685475]

73. Amôr, N.G.; Buzo, R.F.; Ortiz, R.C.; Lopes, N.M.; Saito, L.M.; Mackenzie, I.C.; Rodini, C.O. In vitro and in vivo characterization of cancer stem cell subpopulations in oral squamous cell carcinoma. J. Oral Pathol. Med.; 2020; in press

74. Bocci, F.; Gearhart-Serna, L.; Boareto, M.; Ribeiro, M.; Ben-Jacob, E.; Devi, G.R.; Levine, H.; Onuchic, J.N.; Jolly, M.K. Toward understanding cancer stem cell heterogeneity in the tumor microenvironment. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 148-157. [DOI: https://dx.doi.org/10.1073/pnas.1815345116]

75. Roylance, R.; Endesfelder, D.; Gorman, P.; Burrell, R.A.; Sander, J.; Tomlinson, I.; Hanby, A.M.; Speirs, V.; Richardson, A.L.; Birkbak, N.J. et al. Relationship of extreme chromosomal instability with long-term survival in a retrospective analysis of primary breast cancer. Cancer Epidemiol. Biomark. Prev.; 2011; 20, pp. 2183-2194. [DOI: https://dx.doi.org/10.1158/1055-9965.EPI-11-0343]

76. Birkbak, N.J.; Eklund, A.C.; Li, Q.; McClelland, S.E.; Endesfelder, D.; Tan, P.; Tan, I.B.; Richardson, A.L.; Szallasi, Z.; Swanton, C. Paradoxical relationship between chromosomal instability and survival outcome in cancer. Cancer Res.; 2011; 71, pp. 3447-3452. [DOI: https://dx.doi.org/10.1158/0008-5472.CAN-10-3667]

77. Chen, W.; Yi, J.K.; Shimane, T.; Mehrazarin, S.; Lin, Y.L.; Shin, K.H.; Kim, R.H.; Park, N.H.; Kang, M.K. Grainyhead-like 2 regulates epithelial plasticity and stemness in oral cancer cells. Carcinogenesis; 2016; 37, pp. 500-510. [DOI: https://dx.doi.org/10.1093/carcin/bgw027]

78. Kim, D.; Choi, B.; Ryoo, I.; Kwak, M.-K. High NRF2 level mediates cancer stem cell-like properties of aldehyde dehydrogenase (ALDH)-high ovarian cancer cells: Inhibitory role of all-trans retinoic acid in ALDH/NRF2 signaling. Cell Death Dis.; 2018; 9, 896. [DOI: https://dx.doi.org/10.1038/s41419-018-0903-4]

79. Jia, W.; Tripathi, S.; Chakraborty, P.; Chedere, A.; Rangarajan, A.; Levine, H.; Jolly, M.K. Epigenetic feedback and stochastic partitioning during cell division can drive resistance to EMT. Oncotarget; 2020; 11, pp. 2611-2624. [DOI: https://dx.doi.org/10.18632/oncotarget.27651]

80. Subbalakshmi, A.R.; Sahoo, S.; Biswas, K.; Jolly, M.K. A computational systems biology approach identifies SLUG as a mediator of partial Epithelial-Mesenchymal Transition (EMT). bioRxiv; 2020; [DOI: https://dx.doi.org/10.1101/2020.09.03.278085]

81. Jolly, M.K.; Ware, K.E.; Gilja, S.; Somarelli, J.A.; Levine, H. EMT and MET: Necessary or permissive for metastasis?. Mol. Oncol.; 2017; 11, pp. 755-769. [DOI: https://dx.doi.org/10.1002/1878-0261.12083] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28548345]

82. Jolly, M.K.; Levine, H. Computational systems biology of epithelial-hybrid-mesenchymal transitions. Curr. Opin. Syst. Biol.; 2017; 3, pp. 1-6. [DOI: https://dx.doi.org/10.1016/j.coisb.2017.02.004]

83. Devaraj, V.; Bose, B. Morphological State Transition Dynamics in EGF-Induced Epithelial to Mesenchymal Transition. J. Clin. Med.; 2019; 8, 911. [DOI: https://dx.doi.org/10.3390/jcm8070911] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31247884]

84. Xin, Y.; Cummins, B.; Gedeon, T. Multistability in the epithelial-mesenchymal transition network. BMC Bioinform.; 2020; 21, 71. [DOI: https://dx.doi.org/10.1186/s12859-020-3413-1]

85. Shomar, A.; Barak, O.; Brenner, N. Local and global features of genetic networks supporting a phenotypic switch. PLoS ONE; 2020; 15, e0238433. [DOI: https://dx.doi.org/10.1371/journal.pone.0238433]

86. Sha, Y.; Wang, S.; Zhou, P.; Nie, Q. Inference and multiscale model of epithelial-to-mesenchymal transition via single-cell transcriptomic data. Nucleic Acids Res.; 2020; 48, pp. 9505-9520. [DOI: https://dx.doi.org/10.1093/nar/gkaa725]

87. Jia, D.; Lu, M.; Jung, K.H.; Park, J.H.; Yu, L.; Onuchic, J.N.; Kaipparettu, B.A.; Levine, H. Elucidating cancer metabolic plasticity by coupling gene regulation with metabolic pathways. Proc. Natl. Acad. Sci. USA; 2019; 116, pp. 3909-3918. [DOI: https://dx.doi.org/10.1073/pnas.1816391116]

88. Xu, S.; Ware, K.; Ding, Y.; Kim, S.; Sheth, M.; Rao, S.; Chan, W.; Armstrong, A.; Eward, W.; Jolly, M. et al. An Integrative Systems Biology and Experimental Approach Identifies Convergence of Epithelial Plasticity, Metabolism, and Autophagy to Promote Chemoresistance. J. Clin. Med.; 2019; 8, 205. [DOI: https://dx.doi.org/10.3390/jcm8020205]

89. Tavassoly, I.; Parmar, J.; Shajahan-Haq, A.N.; Clarke, R.; Baumann, W.T.; Tyson, J.J. Dynamic modeling of the interaction between autophagy and apoptosis in mammalian cells. CPT Pharmacomet. Syst. Pharmacol.; 2015; 4, pp. 263-272. [DOI: https://dx.doi.org/10.1002/psp4.29]

90. Gómez Tejeda Zañudo, J.; Scaltriti, M.; Albert, R. A network modeling approach to elucidate drug resistance mechanisms and predict combinatorial drug treatments in breast cancer. Cancer Converg.; 2017; 1, 5. [DOI: https://dx.doi.org/10.1186/s41236-017-0007-6]

91. Goldman, A.; Majumder, B.; Dhawan, A.; Ravi, S.; Goldman, D.; Kohandel, M.; Majumder, P.K.; Sengupta, S. Temporally sequenced anticancer drugs overcome adaptive resistance by targeting a vulnerable chemotherapy-induced phenotypic transition. Nat. Commun.; 2015; 6, 6139. [DOI: https://dx.doi.org/10.1038/ncomms7139] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25669750]

92. Tripathi, S.C.; Peters, H.L.; Taguchi, A.; Katayama, H.; Wang, H.; Momin, A.; Jolly, M.K.; Celiktas, M.; Rodriguez-Canales, J.; Liu, H. et al. Immunoproteasome deficiency is a feature of non-small cell lung cancer with a mesenchymal phenotype and is associated with a poor outcome. Proc. Natl. Acad. Sci. USA; 2016; 113, pp. E1555-E1564. [DOI: https://dx.doi.org/10.1073/pnas.1521812113] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26929325]

93. Celià-Terrassa, T.; Bastian, C.; Liu, D.D.; Ell, B.; Aiello, N.M.; Wei, Y.; Zamalloa, J.; Blanco, A.M.; Hang, X.; Kunisky, D. et al. Hysteresis control of epithelial-mesenchymal transition dynamics conveys a distinct program with enhanced metastatic ability. Nat. Commun.; 2018; 9, 5005. [DOI: https://dx.doi.org/10.1038/s41467-018-07538-7] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30479345]

94. Milo, R.; Jorgensen, P.; Moran, U.; Weber, G.; Springer, M. BioNumbers—The database of key numbers in molecular and cell biology. Nucleic Acids Res.; 2010; 38, pp. D750-D753. [DOI: https://dx.doi.org/10.1093/nar/gkp889]

95. Duddu, A.S.; Sahoo, S.; Hati, S.; Jhunjhunwala, S.; Jolly, M.K. Multistability in cellular differentiation enabled by a network of three mutually repressing master regulators. J. R. Soc. Interface; 2020; 17, 20200631. [DOI: https://dx.doi.org/10.1098/rsif.2020.0631]

96. Rousseeuw, P.J. Silhouettes: A graphical aid to the interpretation and validation of cluster analysis. J. Comput. Appl. Math.; 1987; 20, pp. 53-65. [DOI: https://dx.doi.org/10.1016/0377-0427(87)90125-7]

97. Hardy, A. On the number of clusters. Comput. Stat. Data Anal.; 1996; 23, pp. 83-96. [DOI: https://dx.doi.org/10.1016/S0167-9473(96)00022-9]

98. Caliñski, T.; Harabasz, J. A Dendrite Method for Cluster Analysis. Commun. Stat.-Theory Methods; 1974; 3, pp. 1-27. [DOI: https://dx.doi.org/10.1080/03610927408827101]

99. Davies, D.L.; Bouldin, D.W. A Cluster Separation Measure. IEEE Trans. Pattern Anal. Mach. Intell.; 1979; 1, pp. 224-227. [DOI: https://dx.doi.org/10.1109/TPAMI.1979.4766909]

100. Ward, J.H. Hierarchical Grouping to Optimize an Objective Function. J. Am. Stat. Assoc.; 1963; 58, pp. 236-244. [DOI: https://dx.doi.org/10.1080/01621459.1963.10500845]

101. Mann, H.B.; Whitney, D.R. On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other. Ann. Math. Stat.; 1947; 18, pp. 50-60. [DOI: https://dx.doi.org/10.1214/aoms/1177730491]

102. Welch, B.L. The Generalization of Student’s Problem when Several Different Population Variances are Involved. Biometrika; 1947; 34, pp. 28-35. [DOI: https://dx.doi.org/10.2307/2332510] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20287819]