Research ArticleCell Biology

TGF-β–induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops

See allHide authors and affiliations

Sci. Signal.  30 Sep 2014:
Vol. 7, Issue 345, pp. ra91
DOI: 10.1126/scisignal.2005304

Stepping Through the Transitions in Cell State

Cells can change their functional and morphological characteristics in response to changes in the environment. Epithelial cells can exhibit a type of plasticity called epithelial-to-mesenchymal transition (EMT), in which they transform from cells with tight contacts to neighboring cells and function as a barrier into motile cells that are not connected to their neighbors. This process is important during physiological processes, such as wound healing, and can also contribute to metastatic progression of cancer. Zhang et al. performed population and single-cell analysis of changes in the abundance of core regulators of EMT in a mammary epithelial cell line culture to test theoretical models by which EMT could occur. Their analysis confirmed that two sequential feedback loops, involving transcription factors and microRNAs, function as two cascading switches to control the discrete steps in EMT and identified the first step in this two-step process as the only readily reversible one for these cells.


The process of epithelial-to-mesenchymal transition (EMT) is an essential type of cellular plasticity associated with a change from epithelial cells that function as a barrier consisting of a sheet of tightly connected cells to cells with properties of mesenchyme that are not attached to their neighbors and are highly motile. This phenotypic change occurs during development and also contributes to pathological processes, such as cancer progression. The molecular mechanisms controlling the switch between the fully epithelial and fully mesenchymal phenotypes and cells that have characteristics of both (partial EMT) are controversial, and multiple theoretical models have been proposed. To test these theoretical models, we systematically measured the changes in the abundance of proteins, mRNAs, and microRNAs (miRNAs) that represent the core regulators of EMT induced by transforming growth factor–β1 (TGF-β1) in the human breast epithelial cell line MCF10A at the population and single-cell levels. We provide experimental confirmation for a model of cascading switches in phenotypes associated with TGF-β1–induced EMT of MCF10A cells that involves two double-negative feedback loops: one between the transcription factor SNAIL1 and the miR-34 family and another between the transcription factor ZEB1 and the miR-200 family. Furthermore, our data showed that whereas the transition from epithelial to partial EMT was reversible for MCF10A cells, the transition from partial EMT to mesenchymal was mostly irreversible at high concentrations of TGF-β1.


Epithelial-to-mesenchymal transition (EMT) plays important roles in normal physiological processes, such as development and wound healing, and also in pathological processes, such as fibrosis and cancer metastasis (1). During EMT, cells gain the ability to invade and migrate through a loss of epithelial characteristics, such as tight cellular adhesions, and acquisition of mesenchymal attributes, such as loss of cellular adhesion, production of secreted proteases, and increased motility (13). In addition to the epithelial (E) and mesenchymal (M) states (phenotypes), there exists an intermediate phenotype in the process of EMT, known as the partial EMT (P) state (410). The P state retains some characteristics of epithelial cells but also shows features of mesenchymal cells (11). Thus, EMT is a three-state transition process (Fig. 1A). Notably, experimental studies suggest that the E-to-P (E→P) transition is reversible, and the P state is a metastable phenotype, whereas the P-to-M (P→M) transition has been reported to be either reversible or irreversible, depending on the cell types (1214). Like EMT, the E→P transition and the reverse process P→E may be relevant to pathological processes, such as cancer metastasis (15, 16).

Fig. 1 Schematic illustration of TGF-β–induced EMT.

(A) Three states of EMT. (B) The core regulatory network of TGF-β–induced EMT. The input of the system is exogenous TGF-β, which induces SNAIL1 expression. SNAIL1 and miR-34 participate in a double-negative feedback loop (17, 18). SNAIL1 also inhibits its own expression (36). SNAIL1 stimulates expression of ZEB1 and inhibits expression of miR-200. Another similar double-negative feedback loop is formed between ZEB1 and miR-200 (1921). Furthermore, miR-200 inhibits the autocrine production of TGF-β (13), forming another feedback loop. E-cadherin, an epithelial cell marker, is inhibited by SNAIL1 and ZEB1, whereas N-cadherin and vimentin, markers of mesenchymal cells, are promoted by SNAIL1 and ZEB1. (C) Different functional roles of SNAIL1 and ZEB1 modules are proposed in the CBS model (blue lines) and TCS model (orange lines). SNAIL1 and miR-34 form a bistable switch (blue line) in the CBS model but a monostable noise filter (orange line) in the TCS model. ZEB1/miR-200 is bistable (blue line) in the CBS model but tristable (orange line) in the TCS model. (D) The predicted bifurcation diagram. The thick blue, green, and red lines correspond to E, P, and M states, respectively. The blue and purple arrow lines show how cell fates flip or maintain when the exogenous TGF-β dose increases or decreases.

Transforming growth factor–β (TGF-β) triggers EMT. The core regulatory network of TGF-β–induced EMT is composed of two transcription factors (SNAIL1 and ZEB1), two families of microRNAs (miRNAs; miR-34 and miR-200), and TGF-β (Fig. 1B). Several feedback loops exist: a SNAIL1/miR-34 double-negative feedback loop (17, 18), a ZEB1/miR-200 double-negative feedback loop (1921), and an autocrine TGF-β/miR-200 feedback loop (13). How these interlinked feedback loops coordinately modulate EMT is unknown. Despite the relatively small number of participants in the core regulatory network, a quantitative and systematic study is needed to understand this complex process, because the involved molecular species are interconnected and the functional roles of each one must be placed in the context of the overall regulatory network.

Two disparate mathematical models were proposed to elucidate the underlying mechanism of TGF-β–induced EMT. Through deterministic analyses and stochastic simulations, Tian et al. proposed the CBS model consisting of “Cascading Bistable Switches” (22). In this model, EMT proceeds after a sequential two-step program, in which an epithelial cell first transits to P then to M state, depending on the strength and duration of TGF-β stimulation. Mechanistically, two bistable switches are coupled to govern the system. The SNAIL1/miR-34 feedback loop is responsible for the first reversible switch and regulates the initiation of EMT (Fig. 1C, left panel, blue lines). The ZEB1/miR-200 feedback loop forms the second switch, controlling the establishment of the M state (Fig. 1C, right panel, blue lines). Furthermore, the autocrine TGF-β/miR-200 feedback loop may contribute to the irreversibility of the second switch by maintaining the mesenchymal state.

Lu et al. proposed the TCS model for “Ternary Chimera Switch” (23). In this model, instead of acting as a binary switch, the SNAIL1/miR-34 loop forms a monostable module and functions as a noise-buffering integrator of internal and external signals (Fig. 1C, left panel, orange line). These authors indicated that the negative autoregulation of SNAIL1 is the biological evidence for this function (Fig. 1B). The ZEB1/miR-200 module, which includes the established double-negative feedback loop and a putative positive autoregulatory loop from ZEB1 to its own gene, forms a ternary switch with three steady states corresponding to the E, P, and M states (Fig. 1C, right panel, orange lines).

Both models predict that the three cell phenotypes will occur through a similar process in response to exogenous stimulation with TGF-β (Fig. 1D), although the molecular mechanisms of generating the cell phenotypes are different. Additionally, the two models lead to qualitatively different predictions regarding the characteristics of the cells in the population exposed to TGF-β (see table S1). The CBS model predicts that upon changing the TGF-β concentration, the abundance of SNAIL1 defines two populations of cells. The SNAIL1low population consists of the E cells, and the SNAIL1high population consists of the P or M cells. The abundance of ZEB1 also defines two populations of cells. The ZEB1low population consists of the E or P cells, and the ZEB1high population consists of the M cells. At certain TGF-β concentrations, cells with low and cells with high SNAIL1 (or ZEB1) coexist in the population. The TCS model predicts that SNAIL1 smoothly increases with TGF-β; thus, at any TGF-β concentration, the abundance of SNAIL1 only defines one population of cells, whereas the abundance of ZEB1 defines three populations of cells: ZEB1low corresponding to the E cells, ZEB1medium corresponding to the P cells, and ZEB1high corresponding to the M cells.

Here, we performed quantitative measurement of the dynamics of EMT induced by TGF-β1 in MCF10A cells to test the predictions from these two models. We identified and characterized the three phenotypes at the single-cell level. Both the temporal and steady-state results revealed two-step dynamics of EMT, which is consistent with both models. However, both the dynamics and response curves of SNAIL1 and ZEB1 showed characteristics of a bistable switch, as predicted by the CBS model. In summary, our findings provide a new mechanistic understanding of the TGF-β–induced EMT.


The MCF10A population response to increasing concentrations of TGF-β1 supports the CBS model of EMT

TGF-β induces EMT in several cell lines, including the human mammary epithelial cell line MCF10A and the dog kidney epithelial cell line MDCK (Madin-Darby canine kidney) (24). We cultured MCF10A cells in medium supplemented with various concentrations of human recombinant TGF-β1 (hTGF-β1) for 7 days, changing the medium every other day to maintain hTGF-β1 concentrations, and analyzed the morphology of the cells (Fig. 2A) and the abundance of molecular markers of epithelial (E-cadherin) or mesenchymal (vimentin) phenotypes (Fig. 2B). Control cells formed a contiguous monolayer and adhered to each other tightly, showing characteristics of epithelial cells (Fig. 2A). In response to increased concentrations of hTGF-β1, an increasing proportion of the cells progressed to spindle-like morphology and lost their tight connections. At the highest concentration of hTGF-β1 tested (4 ng/ml), few cells morphologically resembled epithelia (Fig. 2A). Immunofluorescence analysis showed that without hTGF-β1, abundant E-cadherin (labeled green) localized to the plasma membrane, and little vimentin (labeled red) was detected. At the low concentration of hTGF-β1 (0.5 ng/ml), we observed a mixture of cells with varying abundance of E-cadherin, and some cells also showed faint vimentin staining. With increasing hTGF-β1 concentration, vimentin abundance increased and E-cadherin abundance decreased. These observations confirmed that hTGF-β1 induces EMT in MCF10A cells.

Fig. 2 The TGF-β1 dose response of EMT in MCF10A cells.

(A) Morphology of cells treated with various concentrations of hTGF-β1 as indicated for 7 days, replacing the medium every other day (magnification, ×100). (B) Double immunofluorescence staining of E-cadherin (epithelial marker, green) and vimentin (mesenchymal marker, red) with same treatment as in (A). Nuclei were stained with 4′,6-diamidino-2-phenylindole (DAPI) (blue). (C) RNA abundance of the core regulatory molecules measured by quantitative polymerase chain reaction (PCR). (D) Protein abundance of core transcription factors measured by Western blot. For (C) and (D), data represent means ± SD (n = 3); star denotes P < 0.05 (Mann-Whitney test) and at least threefold changes from the control group.

To explore functions of the key players in EMT, we analyzed changes in the abundance of mRNAs encoding SNAIL1 and ZEB1 or the miRNAs miR-34 and miR-200 in response to the different concentrations of hTGF-β1 (Fig. 2C). In the absence of added hTGF-β1, the cells were high in miR-34a and miR-200c and low in SNAIL1 and ZEB1. In response to either 0.5 or 1 ng/ml hTGF-β1, miR-34 abundance decreased, SNAIL1 abundance increased, and miR-200c and ZEB1 were largely unchanged. With 2 or 4 ng/ml hTGF-β1, the abundance of both miRNAs was low, and the transcripts for the two transcription factors, SNAIL1 and ZEB1, were abundant (Fig. 2C). The overall transcript profile indicated two steps, with each step accompanied by an abrupt increase in the mRNA of the transcription factor and a reduction in the abundance of the miRNA. Measuring protein abundance in these populations of cells showed a trend of change consistent with the changes in the mRNA abundance of the two transcription factors (Fig. 2D), with SNAIL1 changing at the first step and ZEB1 changing at the second step. Consistent with the immunofluorescence data (Fig. 2B), the abundance of E-cadherin decreased at the same concentrations at which vimentin increased. As predicted by both the CBS and TCS models, these cell population data indicated that TGF-β–induced EMT is a two-step process that depends on the concentration of TGF-β. The changes in the abundance of the transcription factors, their transcripts, and the miRNAs with varying TGF-β concentrations and the lack of change in the abundance of miR-200 and ZEB1 at the first step are consistent with the predictions of the CBS model (22).

Three steady states corresponding to epithelial and partial EMT and mesenchymal phenotypes are distinguishable during TGF-β–induced EMT

The experiments performed so far represent population averages that could mask coexisting phenotypes, that is, the presence of a small proportion of cells that are different from the majority of cells in the population. To confirm the existence of three cell fates in hTGF-β1–induced EMT, we exposed cells to different concentrations of hTGF-β1 for specified periods, fixed and stained the cells with fluorescent antibodies recognizing E-cadherin and vimentin, and then performed flow cytometric analysis. These data revealed four distinct groups of cells when the abundance of E-cadherin was plotted against the abundance of vimentin (Fig. 3A). Most of the cells were present in the following three regions: E-cadherinhigh/vimentinlow, E-cadherinmedium/vimentinmedium, and E-cadherinlow/vimentinhigh, corresponding to the E, P, and M states, respectively. In the absence of exogenous hTGF-β1, most cells were in the E state, E-cadherinhigh/vimentinlow. In the presence of 0.5 ng/ml hTGF-β1, two distinct groups were detected: Some of the cells were E-cadherinhigh/vimentinlow, consistent with epithelial characteristics, and a second group of cells had less E-cadherin and exhibited an increase in vimentin abundance, consistent with cells in the P state. As the cells were exposed to increasing concentrations of hTGF-β1, an increasing proportion of the cells exhibited the mesenchymal characteristic represented by E-cadherinlow/vimentinhigh with a concordant decrease in the proportion of cells in the epithelial cell E-cadherinhigh/vimentinlow group. Exposure of cells to 4.0 ng/ml hTGF-β1 resulted in a near-complete population of cells in the M state. We repeated the analysis with cells cultured for 7, 10, or 15 days and obtained similar fractions of cells for each population (Fig. 3B), indicating that the cells had reached a steady-state distribution in each experimental condition.

Fig. 3 Three steady states corresponding to epithelial and partial EMT and mesenchymal phenotypes are distinguishable during hTGF-β1–induced EMT.

(A) The abundance of E-cadherin and vimentin was analyzed by flow cytometry. Cells were treated with the indicated concentrations of hTGF-β1 for 7 days, replacing the medium every other day, and double-stained with antibodies recognizing E-cadherin and vimentin. Pictures are representative of three biological replicates. (B) Dependence of the relative percentages of three subpopulations of cells on the dose of hTGF-β1 cultured for 7, 10, and 15 days (means ± SD, n = 3). (C) Protein abundance of subpopulations of cells sorted by fluorescence-activated cell sorting and then analyzed by Western blot. Data are representative of three experiments.

To confirm the protein abundances observed by flow cytometry for each group, we isolated a fraction of cells from each group and performed Western blot analysis (Fig. 3C). E-cadherin abundance decreased and the vimentin abundance increased as the cells were classified as E, P, and M states. The SNAIL1 abundance was low in cells classified as E state and high both in cells classified as P or M state by flow cytometry. In comparison, the ZEB1 abundance was high only for cells classified as M state. These results further confirmed the two-step process of EMT and indicated that SNAIL1 is important for the transition from E to P, whereas ZEB1 is important for the transition from P to M, consistent with the CBS model.

EMT involves a reversible switch followed by an irreversible switch

The results above indicated that MCF10A cells, starting from the epithelial phenotype, converted to partial EMT and then to mesenchymal phenotypes as hTGF-β1 concentration was increased, consistent with the inducing branches of the bifurcation diagram predicted by the CBS model (blue arrow lines in Fig. 1D). To investigate whether the predicted reversing branches (purple arrow lines in Fig. 1D) were also correct, we performed “reculture” experiments. MCF10A cells were incubated with a specific concentration of hTGF-β1 for 7 days (replacing the hTGF-β1–containing medium every other day), then the medium was removed, and the cells were washed and recultured with lower concentrations of hTGF-β1 for another 10 days. The abundance of E-cadherin and vimentin served as indicators of the E, P, and M states and were measured on the 7th day after the exposure to first concentration of hTGF-β1 and on the 10th day after the exposure to the subsequent reduced concentration of hTGF-β1 (Fig. 4A). Pretreatment with 0.1 ng/ml hTGF-β1 did not induce any noticeable change in the abundance of E-cadherin or vimentin, indicating that this concentration was insufficient to induce EMT or the transition from E to P. Cells pretreated with 0.5 and 1 ng/ml hTGF-β1 maintained the E-cadherin and vimentin coexpression pattern, indicative of the P state, after 17 days in the presence of either of these two concentrations of hTGF-β1, but reverted back to the epithelial phenotype—high E-cadherin and undetectable vimentin—when recultured in the absence of hTGF-β1. Thus, the transition between E and P phenotypes was reversible. In comparison, the high-vimentin and low–E-cadherin abundance pattern of cells pretreated with 4 ng/ml hTGF-β1 was well maintained even when recultured in the absence of hTGF-β1, indicating that the P-to-M transition was largely irreversible for MCF10A cells. Plotting the Western blot data of vimentin abundance onto the bifurcation diagram from the theoretical analysis showed that the experimental results are well explained by the CBS model (Fig. 4B).

Fig. 4 hTGF-β1–induced EMT in MCF10A is composed of a reversible switch between the E and P states and a largely irreversible switch between the P and M states.

(A) Abundance of E-cadherin and vimentin in MCF10A cells treated first with various concentrations of hTGF-β1 for 7 days and then with the indicated concentrations of hTGF-β1 for another 10 days. (B) Comparison between measured and predicted dose-response curves of vimentin over hTGF-β1. Measured data are means ± SD obtained from three replicates. (C) Flow cytometric analysis of the sorted cells in the P state with 10 additional days of culture with various concentrations of hTGF-β1. (Upper panel) Cells were initially exposed to hTGF-β1 (0.5 ng/ml) for 7 days and then sorted using the live-cell staining procedure, and cells from the E-cadherinmedium/vimentinmedium group were exposed to the indicated concentrations of hTGF-β1, then restained for vimentin and E-cadherin and analyzed by flow cytometry. (Lower panel) Cells were sorted from the E-cadherin low/vimentin high group that were treated with 4 ng/ml TGF-β for 7 days, and then cultured in the absence of added TGF-β for 10 days. Figures are representative of three independent experiments. (D) Subpopulation changes in cell cultures that have been treated with TGF-β inhibitors or control. Cells were pretreated with 4 ng/ml hTGF-β1 for 7 days, then TGF-β1 was removed, and SB-431542, SB-525334, or dimethyl sulfoxide (DMSO) was added in the culture. Nothing was added into the No-DMSO group but culture medium. Cells were harvested and analyzed after 20 days. Data represent means ± SD obtained from three replicates.

To further confirm the reversibility at the single-cell level, we sorted living cells in the P or M state by flow cytometry and recultured them for additional 10 days with various concentrations of hTGF-β1 (Fig. 4C). About half of the sorted P cells were in the E-cadherinhigh/vimentinlow group, representing cells in the E state, when recultured with 0.5 ng/ml hTGF-β1, and most cells transited back to the epithelial region when recultured in the absence of hTGF-β1 (Fig. 4C, upper graphs). In contrast, most of the sorted M cells remained E-cadherinlow/vimentinhigh even after 3 weeks of reculturing in the absence of hTGF-β1, and none of the cells were E-cadherinhigh/vimentinlow (Fig. 4C, lower graph and Fig. 4D), indicating that unlike the P state, the M state was not readily reversible.

The CBS model suggests that the autocrine TGF-β feedback loop contributes to maintenance of the M phenotype. We performed reculture experiments with cells initially in the M state and then recultured in the presence of TGF-β receptor 1 (TGFBRI) inhibitors and in the absence of hTGF-β1. The addition of the TGFBRI inhibitors SB-431542 or SB-525334 resulted in an increase in the proportion of cells that reverted to the P or E state compared to the control cells that were not exposed to the inhibitors but were recultured in the absence of hTGF-β1 (Fig. 4D). The result is consistent with a study of EMT in MDCK cells, which also exhibited an M→E conversion in the presence of TGFBRI inhibitors (13). Therefore, the data and the computational analysis indicated that TGF-β1–induced EMT of MCF10A cells is composed of a readily reversible step between the E and P states and a much less reversible step between the P and M states.

Flow cytometry data and theoretical analysis confirmed that each of the SNAIL1/miR-34 and ZEB1/miR-200 modules forms a bistable switch

Both the CBS and TCS models predicted a two-step EMT process. The CBS model predicted that each of the modules, SNAIL1/miR-34 and ZEB1/miR-200, was responsible for one step, respectively. The TCS model predicted that only the ZEB1/miR-200 module was responsible for the two steps; thus, one should observe populations of cells with at least three different abundance states of ZEB1 and miR-200 during EMT. Therefore, we measured the distribution of SNAIL1 and ZEB1 with various hTGF-β1 concentrations using flow cytometry.

When cultured for 7 days in the absence of hTGF-β1, flow cytometric analysis showed that the cells formed a single cluster when stained for SNAIL1 and that the addition of 0.5 ng/ml hTGF-β1 resulted in two clusters of cells with different amounts of SNAIL1 (Fig. 5A). This distribution of two subpopulations was also evident when the cells were exposed to 1 ng/ml hTGF-β1, but a single population of cells with high SNAIL1 abundance emerged with increasing concentrations of hTGF-β1. Comparison with results in Fig. 3A, which show the abundance of the E, P, and M state markers, revealed that SNAIL1low and SNAIL1high cells coexisted only when cells in both the E and P states coexisted (0, 0.5, and 1.0 ng/ml hTGF-β1).

Fig. 5 Flow cytometry data and theoretical analysis confirm that the SNAIL1/miR-34 module and the ZEB1/miR-200 module each form a bistable switch.

(A and B) Flow cytometric analysis of SNAIL1 abundance (A) and ZEB1 abundance (B) in MCF10A cells exposed to the indicated concentration of hTGF-β1 for 7 days. The y axis reflects the front scatter characteristics (FSC). Data are plotted on a linear scale (FCS, y axes) and a log scale (protein abundance, x axes). (C) Left: Nullclines of SNAIL1 and miR-34. Right: Bifurcation diagram of SNAIL1 as a function of the concentration of inducing TGF-β. (D) Left: Nullclines of ZEB1 and miR-200. Right: Bifurcation diagram of ZEB1 as a function of the abundance of SNAIL1.

Flow cytometric analyses of ZEB1 abundance showed that the change in ZEB1 abundance did not follow the same pattern as that of SNAIL1 (Fig. 5B). A single population with a low abundance of ZEB1 existed when the cells were cultured in the absence of hTGF-β1 or in the presence of 0.5 ng/ml hTGF-β1. Similarly, a single population with a high abundance of ZEB1high existed when the cells were cultured in 4 ng/ml hTGF-β1. Unlike SNAIL1, which showed two coexisting subpopulations at 0.5 and 1 ng/ml hTGF-β1, coexisting subpopulations of cells with high and low abundance of ZEB1 were present in cells cultured with 1 or 2 ng/ml hTGF-β1. These results—especially the observation that cells exposed to 0.5 ng/ml hTGF-β1, a condition that corresponds to those of cells in the E and P states, coexist as SNAIL1low and SNAIL1high but only ZEB1low—agree with the CBS model that SNAIL1/miR-34 and ZEB1/miR-200 form two separate binary switches.

We also tested for roles in reversibility by performing reculture experiments. Cells were cultured with various concentrations of hTGF-β1 for 7 days, then recultured in the absence of hTGF-β1 for another 10 days and analyzed SNAIL1 abundance by flow cytometry (fig. S1A). Cells cultured in the absence of hTGF-β1 after culture in conditions that produced a coexisting population of SNAIL1high and SNAIL1low cells (0.5 ng/ml hTGF-β1) (Fig. 5A) resulted in complete loss of the SNAIL1high population and single population of SNAIL1low (fig. S1A). However, cells pretreated with 1 or 2 ng/ml hTGF-β1 maintained both SNAIL1low and SNAIL1high populations when recultured in the absence of hTGF-β1, which is consistent with a population of the cells having transitioned into the M state. Most cells pretreated with 4 ng/ml hTGF-β1 remained as SNAIL1high, with only a small fraction of cells existing as SNAIL1low. Because the results in Fig. 3A indicated that some cells coexisted in the P state and the M state with 1, 2, and even 4 ng/ml hTGF-β1, we hypothesized that the SNAIL1high cells were mixtures of P and M cells under these conditions and that the P but not the M cells converted to form the SNAIL1low population. To confirm this hypothesis, we separated cells in the P or M states (on the basis of E-cadherin and vimentin abundance) after treating the MCF10A cells with 2 ng/ml hTGF-β1 for 7 days and then recultured the cells in the absence of hTGF-β1 for 10 days. Most of the P cells converted back to SNAIL1low, whereas M cells remained as SNAIL1high (fig. S1B).

The agreement between the experimental results and the CBS model prompted us to reexamine several theoretical concerns of the CBS model raised by Lu et al. (23). The CBS model from Tian et al. (22) used a phenomenological Hill function without much molecular detail to represent the miR-34–mediated inhibition of SNAIL1 and did not include the self-inhibition of SNAIL1. Lu et al. explicitly included miRNA-mRNA interactions in the TCS model and concluded that the SNAIL1/miR-34 module might not have sufficient nonlinearity for the bistability predicted by the CBS model. Therefore, we revised the CBS model to include the self-inhibition of SNAIL1 and detailed binding and unbinding reactions involved in the miRNA-mediated inhibition of mRNA (see texts S1 and S2 for details of the model). For the latter, we considered the number of the binding sites explicitly. Indeed, we analyzed the nullclines of the isolated SNAIL1/miR-34 module and identified three intersections corresponding to three steady states of the corresponding ordinary differential equations (Fig. 5C, left). Further mathematical analysis showed that two of them are stable (solid circles), representing the E state and the P or M state, and the other one is unstable (open circle). The bifurcation diagram of SNAIL1 also revealed bistable regions when plotted against exogenous TGF-β concentration, corresponding to the coexistence region between epithelial and partial EMT cells. SNAIL1 transmits the signal to the ZEB1/miR-200 module by promoting ZEB1 transcription (6, 7) and repressing miR-200 transcription. Similar bifurcation analysis of the ZEB1/miR-200 module regulated by SNAIL1 (Fig. 5D) confirmed the possible bistable states of ZEB1 upon changing the SNAIL1 abundance.

Mechanistically, bistability is ensured by nonlinearity and the presence of positive feedback loops. For the SNAIL1/miR-34 module, several molecular mechanisms possibly contribute to the nonlinearity. Two conserved SNAIL1-binding sites are found in the miR-34 promoter (17). Endogenous miR-34a produces a nonlinear threshold response from its target gene in cancer cells (25). Thus, the population and single-cell analyses, along with the theoretical analysis, confirmed that the SNAIL1/miR-34 module and the ZEB1/miR-200 module each form a bistable switch.

Time course measurements revealed two sequential jumps among the three cell types

The steady-state results revealed that MCF10A cells exist in three possible phenotypes when treated with different concentrations of hTGF-β1. In response to a given hTGF-β1 concentration, how did cells transform from one phenotype to another? Did they change gradually, or did they transit sharply among the different phenotypes? How rapid is the transition among the phenotypes? To address these questions, we monitored the temporal dynamics of hTGF-β1–induced EMT.

We exposed MCF10A cells to 4 ng/ml hTGF-β1 for various lengths of time. Morphological analysis showed that the cells changed from being tightly packed and spherical on day 0 to partially losing cell-cell contacts on day 3 and exhibiting a detached spindle shape by day 7 (Fig. 6A). Analysis of transcript and miRNA abundance (fig. S2A) and of protein abundance (fig. S2B) in the population at the different days showed that expression of miR-34a decreased by day 1 and the expression of SNAIL1 started to increase by day 1 and peaked by day 4 (fig. S2A, left). The inhibition of expression of miR-200c and increase in ZEB1 expression occurred on day 5 (fig. S2A, right). The abundance of E-cadherin decreased starting on day 1 and by day 6 was undetectable, whereas the change in the abundance of vimentin displayed a reversed trend (fig. S2B). The trends in the changes in E-cadherin and vimentin abundance were also observed with immunofluorescence staining (fig. S3). These results suggested that temporally TGF-β–induced EMT is also a two-step process. In the first step, cells transit from the E state to the P state, and the SNAIL1/miR-34 bistable switch is activated, resulting in a reduction in E-cadherin and an increase in vimentin. In the second step, cells transit from the P state to the M state, and the ZEB1/miR-200 switch is activated, resulting in the loss of E-cadherin and a robust increase in vimentin.

Fig. 6 Time course measurements revealed two sequential jumps among three cell phenotypes.

(A) Morphological changes of MCF10A cells exposed to 4 ng/ml hTGF-β1 for the time indicated (magnification, ×100). (B) Time-course flow cytometry data of cells with 4 ng/ml hTGF-β1 treatment. Cells were double-stained with antibodies recognizing E-cadherin and vimentin. (C) Temporal changes in the subpopulation percentages during EMT induced with the indicated hTGF-β1 concentrations. Data shown represent the means ± SD, with n = 3 biological replicates.

The two-step dynamics of hTGF-β1–induced EMT was also evident in the single-cell analysis (Fig. 6B). Throughout the time, we observed three clusters of cells located on the E-cadherin–vimentin plane (Fig. 6B) as the clusters seen in the steady-state studies (Fig. 3A). Over time, the percentage of cells in each cluster changed with a sequence from E to P and then from P to M states. This suggests that the dynamics could be well represented by a three-state model, E→P→M.

To examine the effect of increasing concentration of hTGF-β1 on the dynamics of the transitions, we performed time-course flow cytometry measurements with different hTGF-β1 concentrations (fig. S4 and Fig. 6C). Without hTGF-β1, cells remained in the E state. With 0.5 ng/ml hTGF-β1, a fraction of E cells converted to P cells, but none of the cells reached the M state. With higher hTGF-β1 concentrations, all the three states were at least transiently populated. The rate at which the transition from E to P and P to M occurred increased as the concentration of hTGF-β1 increased, but the transition from E to P always included a period when the population included cells in the M state, reflecting a two-step sequential hopping process among the three states.

In summary, our results revealed that hTGF-β1–induced EMT is temporally also a two-step program. The process can be well represented by a three-state model, first transiting from E to P and then to M state, with the fraction of cells in each state and the timing of the transition dependent on the concentration of hTGF-β1.


Three cell phenotypes have been proposed to exist during the process of EMT, with the partial EMT state having both epithelial and mesenchymal features (2629). Here, we performed quantitative measurements of TGF-β1–induced EMT of the MCF10A cell line. The E-cadherin–vimentin flow cytometry data showed that the cells formed three self-organized clusters. Cells of the three clusters had (high, low), (medium, medium), and (low, high) abundance patterns of E-cadherin and vimentin, respectively, consistent with the E, P, and M phenotypes. Furthermore, at the highest concentrations of TGF-β tested, epithelial cells first converted to partial EMT and then to the mesenchymal phenotype, showing sequential stepwise dynamics. Additionally, as the concentration of TGF-β increased and more of the cells were of the M state, the reversibility of the process decreased (Fig. 7). At low TGF-β concentrations, the E and M states are the most stable states but are separated by high barriers and by a metastable P state. Increasing TGF-β concentration destabilizes the E state relative to the P state. At sufficiently high TGF-β concentrations, the M state becomes the most stable one.

Fig. 7 Metaphorical landscape illustration of TGF-β–induced EMT.

Each well represents a stable or metastable phenotype. The lower the well, the more stable (and the more likely to be populated) the phenotype is.

Although existence of the P state is well known, its molecular nature has not been systematically studied and compared to the E or M state. Some early studies called these cells with a “reversible scatter” phenotype that was preferentially induced by certain conditions (30, 31). Our study showed that the P state could be maintained for a long period (at least about 2 weeks) in the presence of midrange concentrations of hTGF-β1 (Fig. 3C), and thus, the P state can be a stable state under some conditions. The two mathematical models of Tian et al. and Lu et al. differ in the molecular mechanism for generating the P state (22, 23), more specifically, whether the SNAIL1/miR-34 module can generate bistability. Our analysis confirmed that a binary SNAIL1/miR-34 switch generates the P state, consistent with the CBS model. We reached this conclusion by evaluating the qualitatively different predictions from the two models, which is experimentally more reliable than comparing quantitative different predictions (see text S3).

The CBS model is also consistent with other reported experimental observations. SNAIL1 is produced at the initiation of EMT, whereas ZEB1 is produced at a later phase (32, 33). Depletion of SNAIL1 prevents the production of ZEB1 (34), impairing EMT (35, 36). Knockdown of ZEB1 in mesenchymal cells partially restores epithelial features (3740). These observations agree with the CBS model that SNAIL1 triggers initiation of EMT by down-regulating E-cadherin transcription, followed by further reduction in E-cadherin transcription by ZEB1 (11, 41) to induce complete EMT and maintain the mesenchymal phenotype (42). Mechanistically, we showed that the SNAIL1/miR-34 module functions as a bistable switch for the transition between E and P states and that the ZEB1/miR-200 module serves as second bistable switch for the transition from the P to M state.

The TCS model emphasizes a possible noise-filtering function of SNAIL1 self-inhibition. Our time course measurements showed that the abundance of SNAIL1 peaked and exhibited a slight decrease in abundance at later times, which could be the result of self-inhibition. This observation is consistent with that of Tran et al. who also used the MCF10A cell line (43). The possible functional roles of this inhibition need further investigation. Self-inhibition can accelerate system response time (44). Alternatively, the transient overshoot of SNAIL1 abundance may contribute to the differential recruitment of other proteins, such as histone modification enzymes, to its target gene’s promoters (45).

Further experiments are needed to examine whether the CBS model accurately represents the mechanism of EMT for other cell types. Whether the sequential bistable modules contribute to partial EMT induced by stimulation with ligands other than TGF-β and the relation between EMT, partial EMT, and tumor invasion and cancer metastasis remain open areas for investigation. Huang et al. examined the protein abundance patterns of 43 ovarian cancer cell lines under normal culturing conditions, without the addition of EMT-inducing stimuli (42). They identified a spectrum of four EMT phenotypes: the epithelial (E), intermediate epithelial (IE), intermediate mesenchymal (IM), and mesenchymal (M) states. The E and IE phenotypes of Huang et al. are similar to the E and P states in this work. Their IM and M states may correspond to a finer classification of the M state as described in this work. The patterns of marker abundance associated with the four states are consistent with the two-stage dynamics that we described here. However, there are some differences between the Huang et al. data and our results. For example, in the cancer cell lines, the SNAIL1 abundance peaked in the IE state but was low in the IM and M states, whereas here, we observed that SNAIL1 abundance remained high in both the P and M states. Besides the potential genetic and epigenetic heterogeneities among different cancer cell lines, another possible explanation is that SNAIL1 is essential for initializing EMT but drops after the mesenchymal phenotype has been fixed. Unlike the established cell lines cultured in the absence of exogenous EMT-inducing stimulus, we studied mesenchymal cells that had been induced by exogenous TGF-β1 for 1 week. Additional events, such as DNA methylation, may occur after the time that we studied (13, 29). One may test this hypothesis by monitoring SNAIL1 abundance after an extended period of culturing the induced mesenchymal cells.

The irreversibility of the M state is conditional. Our study involved the application of exogenous TGF-β1 concentration as the EMT-regulating event. A mesenchymal-to-epithelial transition can take place if other factors, such as the ZEB1/miR-200 ratio (13), are manipulated while removing the exogenous TGF-β1 from the culture medium. In vivo, the situation is even more complex and involves multiple pathways. Additionally, the irreversibility of the M state is likely to be different for different cell types because cells have different functional requirements. Indeed, like the MCF10A cells that we used here, Gregory et al. reported that MDCK cells maintained a mesenchymal morphology for several months after exogenous TGF-β1 removal (13), whereas Gal et al. observed that NMuMG cells reversed to the epithelial phenotype (14). Gregory et al. also reported more efficient reversion of the induced mesenchymal cells to epithelial cells in response to TGF-β1 receptor inhibition of MDCK cells than the reversion that we observed with MCF10A cells. For the MCF10A cells, additional pathways, such as the interleukin 6 (IL-6) signaling pathway (46), may contribute to the maintenance of the mesenchymal phenotype. Differences in the histone and DNA epigenetic states may also contribute to this cell type–dependent stability of the mesenchymal state (29). One may also expect differences between normal and tumor cells. Genetic mutations in tumor cells may lead to different effective model parameters and, thus, different reversibility behaviors.

Here, we focused on the “core” regulatory network of TGF-β–induced EMT. However, biological networks are highly connected and have redundant pathways. For example, ZEB1 and ZEB2 share overlapping functions. ZEB1-knockout mice produced the EMT programs required for entry into the postnatal period (47). Other proteins and miRNAs participate in EMT regulation, such as the miR-205 and miR-30 families (48, 49). Therefore, the SNAIL1/miR-34a bistable switch is probably one of the many regulatory modules that are involved in controlling EMT. Analyzing an expanding network may help reconcile various experimental observations. Additionally, more intermediate phenotypes may be identified if additional markers are used to characterize the process. Indeed, Lee et al. showed bistable states formed by the double-negative feedback loop of the metastasis suppressor RKIP (RAF kinase inhibitory protein) and the transcription factor BACH1 (50). Furthermore, other studies suggest that the EMT dynamics are regulated by a number of pathways, such as the IL-6 and TGF-β pathways (51, 52), and by the three-dimensional environment experienced by the cells (53). Controlled quantitative studies will shed light on how these factors work together to affect EMT.

During EMT, a cell globally changes its genetic program at epigenetic, transcriptional, and translational levels, with relevant processes occurring on the subsecond scale of molecular binding and unbinding events to changes in chromatin structure and posttranslational modification changes occurring in minutes to hours, and changes in gene expression and morphology occurring on a scale of days to weeks (15, 29, 5458). The complexity of the process, in terms of both the number of components involved and the time scales spanned, requires systems biology approaches. In this work, we adopted an integrated mathematical modeling and experimental approach to study EMT dynamics quantitatively. We suggest that this systematic approach can be effective in elucidating how a large number of molecular species orchestrate temporally and spatially to control this important cellular process.


Cell line

The MCF10A cell line was originally obtained from the American Type Culture Collection and was a gift from I. Lazar (Virginia Tech). Cells were maintained in Dulbecco’s modified Eagle’s medium/F12 (1:1) medium (Invitrogen) supplemented with 5% horse serum (Invitrogen), human epidermal growth factor (20 ng/ml) (PeproTech), hydrocortisone (0.5 μg/ml) (Sigma), cholera toxin (0.1 μg/ml) from Vibrio cholerae (Sigma), insulin (10 μg/ml) from bovine pancreas (Sigma), and 1× antibiotic-antimycotic mixture (Gibco). All cells were cultured in 5% CO2–humidified atmosphere at 37°C. MCF10A cells were seeded 1 day before induction and reached to 70 to 80% confluence. Recombinant human TGF-β1 (Cell Signaling) was added to the abovementioned medium at designated concentrations to induce EMT. For reculture experiments, the cells were incubated in medium containing the first concentration of hTGF-β1 for the designated time, then the original medium was removed, the cells were washed with Dulbecco’s phosphate-buffered saline (Sigma) three times, and fresh medium containing the same or a different concentration of TGF-β1 was added. For TGF-β inhibition experiment, the TGF-β receptor inhibitor SB-431542 (5 μM) or SB-525334 (5 μM) (Selleck) was added to the medium. DMSO (0.5 μl/ml) was used as TGF-β receptor inhibitor control. For all experiments, the medium together with the added chemicals or hTGF-β1 was replaced every other day.

Quantitative real-time PCR

Total RNA was isolated using miRNeasy Mini Kit (Qiagen) according to the manufacturer’s protocol. Reverse transcription of genes was performed by using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems). miRNA reverse transcription was performed using the protocol of Kramer et al. (59). Quantitative real-time PCR (RT-PCR) was performed as previously described (60) with the primers shown in table S2. Products were detected with SYBR Green Master Mix (Applied Biosystems) on a Bio-Rad C1000 thermal cycler (Bio-Rad). RNA abundance was normalized against 18S ribosomal RNA. Relative RNA quantification was obtained using the 2−ΔΔCt method (61).

Western blot

MCF10A cells were harvested and lysed by RIPA (radioimmunoprecipitation assay) lysis buffer (Sigma) containing Halt protease inhibitor cocktail (Thermol). Protein solution (40 μg) of each sample was used for Western blot. The targeted proteins were detected with antibodies recognizing the following proteins: actin (Santa Cruz Biotechnology, sc-130656), E-cadherin (Santa Cruz Biotechnology, sc-7870), vimentin (Santa Cruz Biotechnology, sc-5565), SNAIL1 (Santa Cruz Biotechnology, sc-28199), or ZEB1 (Santa Cruz Biotechnology, sc-25388). After incubation with horseradish peroxidase–conjugated bovine anti-rabbit immunoglobulin G (Santa Cruz Biotechnology, sc-2370), the immunolabeled proteins were detected by Western Blotting Luminol Reagent (Santa Cruz Biotechnology) and developed on premium x-ray film (Phenix). Protein concentrations were quantified using Image Studio Lite Software. β-Actin was used to normalize target protein abundance.

Flow cytometry

Flow cytometry studies with MCF10A cells were performed as described (6265) with modifications for cell sorting and culturing. Cells for analysis without the need for culturing after sorting were fixed with 4% paraformaldehyde, incubated with primary antibodies (see Western blot section above) at room temperature for 1 hour, and then incubated with secondary antibody conjugated to fluorescent dye, Alexa Fluor 488 (Santa Cruz Biotechnology, sc-362262) or Alexa Fluor 647 (Santa Cruz Biotechnology, sc-362292), for an additional hour. For cell sorting and culturing after sorting, cells were incubated directly without fixation with E-cadherin antibody conjugated to Alexa Fluor 488 (Santa Cruz Biotechnology, sc-21791; Cell Signaling, 3199) and vimentin antibody conjugated to Alexa Fluor 647 (Santa Cruz Biotechnology, sc-6260). Finally, cells were analyzed, sorted, or both with a BD FACSAria. Flow cytometry data analysis was performed by using FlowJo v10 software.

Immunofluorescence microscopy

MCF10A cells were cultured in eight-chamber slides for immunofluorescence microscopy. Targeted proteins were labeled by using indirect immunofluorescence method (66) and observed using an AxioVision Fluorescence microscope (Carl Zeiss). Images were collected, digitized, and analyzed with AxioVision software.

Statistical methods

Data were first analyzed with the Shapiro-Wilk normality test. For samples that fit the normal distribution, we performed Student’s t test for significant analysis. Otherwise, Mann-Whitney test was performed for significant analysis.

Mathematical modeling

Details of the mathematical modeling approach can be found in the supporting text. The revised CBS model is similar to that of Tian et al. (22) with some key modifications. Inspired by the TCS model, we included the negative autoregulation of SNAIL1. Second, instead of using a Hill function, we modeled the regulation of mRNA by miRNA with mass-action dynamics by considering the miRNA-mRNA complexes and the number of miRNA binding sites on its targets. For simplicity, we used a quasi-equilibrium approximation of miRNA-mRNA binding and unbinding with the conservation conditions of both the total abundance of mRNA and miRNA (see text S1 and figs. S5 and S6 for details of the model for the miRNA regulation of miRNA). The positive autoregulation of ZEB1 is not included, because there is no direct evidence for it. The revised CBS model is described in text S2, and the model parameters and initial conditions of the variables are given in tables S3 and S4, respectively. A computer XPP code is also provided as model S1 (EMT2.ode). Furthermore, text S3 and table S4 provide discussions on qualitatively different predictions that can be used to distinguish between the CBS and TBS models for EMT.


Text S1. General model for the miRNA-mediated regulation of mRNA.

Text S2. The revised CBS model for EMT system.

Text S3. Qualitatively, but not quantitatively, different predictions can distinguish the CBS and TBS models for EMT.

Fig. S1. Flow cytometric analysis of TGF-β1–induced cells recultured in the absence of TGF-β1 confirms that the SNAIL1/miR-34 module functions as a bistable switch.

Fig. S2. Temporal dynamics of TGF-β1–induced EMT in MCF10A cells.

Fig. S3. Immunofluorescence of E-cadherin and vimentin at different time points in response to 4 ng/ml TGF-β1.

Fig. S4. Analyses of the three subgroups at the indicated time points by flow cytometry in response to different concentrations of TGF-β1.

Fig. S5. Generic model for the regulation of mRNA by miRNA.

Fig. S6. Full reaction diagram of miRNA regulation of mRNA with four binding sites when considering all the possible miRNA-mRNA complexes.

Table S1. Comparison between the original CBS and TCS models.

Table S2. Primers used for quantitative RT-PCR.

Table S3. Parameters used in the present model studies.

Table S4. Initial conditions of the variables in the model.

Model S1. Computer code of the theoretical model.

References (6771)


Acknowledgments: We thank I. Lazar for providing the MCF10A cell line, M. Makris for flow cytometry, and S. Kumar, J. J. Tyson, G. Yao, X. Shen, and S. Weiss for many helpful discussions. Funding: The research was supported by the U.S. National Science Foundation (DMS-0969417) and the Fralin Life Sciences Institute at Virginia Tech. Author contributions: J.X. conceived the study. J.X., X.-J.T., J.Z., H.Z., Y.T., F.B., and S.E. designed the study. J.Z. and R.L. conducted the experiments, and X.-J.T. performed the model studies. X.-J.T., J.Z., H.Z., Y.T., J.X., and S.E. analyzed the data. X.-J.T., J.Z., and J.X. wrote the paper with inputs from all authors. Competing interests: The authors declare that they have no competing interests.

Correction: Five parameters in Table S3 and in Model S1 are updated to reflect the values used to produce Figure 4B.

View Abstract

Navigate This Article