The epithelial-to-mesenchymal transition (EMT) of primary cancer contributes to the acquisition of lethal properties, including metastasis and drug resistance. Blocking or reversing EMT could be an effective strategy to improve cancer treatment. However, it is still unclear how to achieve complete EMT reversal (rEMT), as cancer cells often transition to hybrid EMT states with high metastatic potential. To tackle this problem, we employed a systems biology approach and identified a core-regulatory circuit that plays the primary role in driving rEMT without hybrid properties. Perturbation of any single node was not sufficient to completely revert EMT. Inhibition of both SMAD4 and ERK signaling along with p53 activation could induce rEMT in cancer cells even with TGFβ stimulation, a primary inducer of EMT. Induction of rEMT in lung cancer cells with the triple combination approach restored chemosensitivity. This cell-fate reprogramming strategy based on attractor landscapes revealed potential therapeutic targets that can eradicate metastatic potential by subverting EMT while avoiding hybrid states.

Significance:

Network modeling unravels the highly complex and plastic process regulating epithelial and mesenchymal states in cancer cells and discovers therapeutic interventions for reversing epithelial-to-mesenchymal transition and enhancing chemosensitivity.

### Graphical Abstract

Cancer is the second leading cause of death worldwide (1) and metastasis is closely associated with such cancer-related deaths in patients because high plasticity of cancer cells contributes to local invasion, dissemination, and distant metastasis (2). A critical process leading to such plasticity is the epithelial-to-mesenchymal transition (EMT; ref. 3).

The EMT process starts with epithelial cells losing polarity and their ability to adhere to each other and the extracellular matrix (4, 5). These changes result in the cells acquiring properties of mesenchymal cells, which includes invasive and migratory properties. Rather than binary state transitions between epithelial and mesenchymal states, it is now appreciated that the process of EMT results in multiple cellular states, including hybrid states in which cells have properties of both epithelial and mesenchymal cells (6, 7). Indeed, partial EMT generates a spectrum of hybrid cellular states with varying proportions of epithelial and mesenchymal features (8–10). When partial EMT occurs in cancer cells, the cells exhibit heterogeneity at the gene expression level as well as the phenotype level. Cancer cells displaying hybrid phenotypes have a wide range of stemness levels and collectively migrate in small clusters, which are associated with drug resistance and metastasis, respectively (11). Due to the plasticity of EMT, improper drug perturbations have resulted in unexpected hybrid states, which still retain high stemness and metastatic potential (12).

In addition, extracellular signals provided by the tumor microenvironment (TME) can enhance the occurrence of EMT (13) and contribute to the induction of partial EMT (14). Among many extracellular factors, TGFβ is often secreted from cancers as well as other adjacent cells of TME and is known as a primary inducer of EMT (15). EMT is controlled by this extracellular signal (16), and its behavior and interactions become more complicated after drug perturbations, therefore controlling EMT cells can be led to further unexpected results after treatment (17).

In this study, we constructed a logical regulatory network model to tackle the complexity of the EMT process. By implementing the EMT network model and conducting complex network analysis, we systematically analyzed relative stability of attractor states and unraveled feedback loops (FBL), which are critical in the transition to hybrid states. With this approach, we found key regulators that govern such core-regulatory circuit and induce transition to a reversed EMT (rEMT) state by avoiding multiple hybrid states. We predicted and experimentally validated that inhibition of SMAD4 and the kinases ERK1 and ERK2 (collectively ERKs) in the presence of p53 activation can induce rEMT effectively and overcome chemotherapeutic drug resistance in lung cancer cells regardless of TGFβ input signal. Our systems study for cell-fate reprogramming provides insights into the identification of therapeutic targets that can optimize the effectiveness of chemotherapy by avoiding undesired or unexpected hybrid states in cancer.

### Constructing a logical network model for EMT progression

The EMT regulatory network model was constructed on the basis of key proteins and their interactions from signaling databases and published literature (18–20). Distinct molecules in our network model were selected and categorized into eight signaling pathways according to their functions in biological systems. We also described the desired epithelial and mesenchymal phenotypes by their corresponding node activities. Detailed methods used to select major signaling pathways and their representative nodes within our network model are described in Supplementary Text S1. Not all edges were direct, thus interaction modes were categorized according to the relevant literature. Detailed information for each node and edge are provided in Supplementary Table S1 and Supplementary Data S1. We developed a mathematical model using a Boolean network modeling framework. In the Boolean model, the value of each node represents its activity by ‘ON’ or ‘OFF’ for active or inactive state, respectively. Our Boolean network model was established on the basis of our experimental data and other information from relevant literature. Boolean logic reconstruction and validation are described in detail in Supplementary Text S2. Significantly, our model can explain the regulation of epithelial cell plasticity through p53 and p38/AP-1 pathways (detailed comparison with other studies is described in Supplementary Text S3). All of the Boolean functions were written in accordance with the BoolNet format in Supplementary Table S2.

### Robustness of random sampling

Due to the high computational complexity of considering all initial states (231), we randomly sampled the initial states. Detailed methods for the robustness of random sampling are in Supplementary Method S1.

### Attractor landscape analysis

The network state is determined by the set of node states. Initial states repeatedly transition within the state space over time and eventually converge to an equilibrium state called an attractor. Attractor landscape consists of all attractors and their basin of attraction represented by a set of states with their transition trajectories converging to each attractor. With this, we computed the activities of each node by averaging their values in the attractors, weighted by the corresponding basin sizes. Such averaged activity of each node can represent the gene expression level or protein activity. Detailed methods for attractor landscape analysis are in Supplementary Method S2. Because attractors are corresponded to distinct cellular phenotypes (21, 22), we describe the epithelial, hybrid, and mesenchymal phenotypes using E-cadherin (E-cad) and ZEB1 nodes as markers as follows: epithelial (E-cad+/ZEB1), hybrid (E-cad+/ZEB1+, E-cad/ZEB1), and mesenchymal (E-cad/ZEB1+). During Boolean network simulations, each node state was updated synchronously according to its logical rule (more discussed in Supplementary Text S5). We used a Python package ‘PyBoolNet’ with Python 3.7 for such Boolean simulation.

### Qualitative individual input–output relationship

To elucidate the biological phenomena represented within our network model, we qualitatively investigated individual input–output relationships. Detailed methods for individual input–output relationship are in Supplementary Method S3.

### Perturbation analysis

Using our Boolean network model, we simulated knockout (KO) or overexpression (OE) perturbations by fixing the values of nodes to ‘OFF’ or ‘ON’. To mimic our experimental model, all of our main simulation results were obtained while fixing Ras to ON, reflective of the genetic mutation profile of A549 cell line. Detailed methods for perturbation analysis are provided in Supplementary Method S4. After simulating the perturbation effects in our network, we computed an ‘rEMT score’ based on cosine-similarity between each vector and a predefined vector for the desired rEMT state, which consists of average node activities on nine nodes within our network model that are functionally associated with EMT state; E-cad, miR-200, miR-34, ZEB1, Snail, Twist1, EpCAM, c-Myc, and Thy-1, as follows:
where each ZEB1 and E-cad is the average node activity of its node, respectively. x is the average node activity profiles of nine nodes, and y is the desired node activity profiles of the nine nodes in Supplementary Table S1.

### Analysis of network dynamics

After converging into an attractor, each node state is either zero or one. By integrating its regulatory relationship between the nodes, each edge function can be characterized as frustrated or unfrustrated (23). From that, we computed the frustration of a major attractor x, following an equation of |$\sum\nolimits_{i,j} {{F_e}} ( {{x_{ij}}} ) = \sum\nolimits_{i,j} { - {J_{ij}}{x_i}{x_j}}$| (23), where J is regulatory relationship matrix, and xi is the state value of node i. We further demonstrated this process using an example network in Supplementary Fig. S1A and S1B. Then we defined ‘molecular state ambiguity’, as the sum of normalized frustration influence of incoming edges, rin(xj), and outgoing edges, rout(xj), for node j (Supplementary Fig. S1C and S1D). Detailed methods for this simulation are in Supplementary Method S5.

### Analysis of regulatory dynamics with edge function representation

To unravel a subnetwork or interlinked FBLs (24) that is mainly involved in phenotype transition, we analyzed regulatory dynamics under perturbation. Detailed methods for the representation of regulatory dynamics after perturbations are in Supplementary Method S6.

### Survival analysis

Detailed methods for survival analysis of lung patient samples are in Supplementary Method S7.

### Analysis of cellular response to chemotherapy

Detailed methods for analysis of cellular response to chemotherapy are in Supplementary Method S8.

### Reagents and antibodies

Detailed methods for reagents and antibodies are in Supplementary Method S9.

### Cell culture

A549 (KCLB No. 10185) and A427 (KCLB No. 30053) cells were obtained from Korean Cell Line Bank (KCLB). Cells were cultured in DMEM (WelGENE Inc.) with 10% FBS (WelGENE Inc.) and antibiotics (100 units/mL of penicillin, 100 μg/mL streptomycin, and 0.25 μg/mL of Fungizone; Life Technologies Corp.) at 37°C in a humidified atmosphere containing 5% CO2. Cells were periodically tested for Mycoplasma contamination (e-Myco plus Mycoplasma PCR Detection Kit, iNtRON Biotechnology).

### Virus production and establishment of stable cell lines

Detailed methods for virus production and establishment of stable cell lines are in Supplementary Method S10.

### Total RNA extraction and qRT-PCR

Detailed methods for total RNA extraction and qRT-PCR are in Supplementary Method S11 and Supplementary Table S5.

### Western blot analysis

Cells were lysed in lysis buffer [20 mmol/L HEPES, pH 7.2, 150 mmol/L NaCl, 0.5% Triton X-100, 10% glycerol, protease/phosphatase inhibitor cocktail (Thermo Fisher)]. The lysates were centrifuged at 13,000 rpm for 20 minutes at 4°C and the supernatants were separated by SDS-PAGE, followed by Western blot analysis.

### Cell growth and viability analysis

To analyze the effect of chemotherapeutic drugs on cell viability, transfected cells were seeded in 96-well plate (4 × 103 cells/well). After 24 hours, drugs [doxorubicin (DOX), 5-fluoroucil (5-FU), Oxaliplatin (OX), or irinotecan (CPT-11)] were treated with a combination of TGFβ (5 ng/mL) without or with Nutlin-3a (1 μmol/L, Sigma) for 5 days. After seeding, cells were imaged using IncuCyte ZOOM (Sartorius). Cell growth or viability was analyzed from confluence measurement of adherent living cells using IncuCyte software.

### Crystal violet assay

Detailed methods for crystal violet assay are in Supplementary Method S12.

### Data availability

The Genomics of Drug Sensitivity in Cancer dataset (Release 8.1, Oct 2019) was downloaded from https://www.cancerrxgene.org/downloads/bulk_download. The Cancer Genome Atlas (TCGA) datasets used in this study are Lung Adenocarcinoma (TCGA, PanCancer Atlas) dataset: https://www.cbioportal.org/study/summary?id=luad_tcga_pan_can_atlas_2018 and Lung Squamous Cell Carcinoma (TCGA, PanCancer Atlas) dataset: https://www.cbioportal.org/study/summary?id=lusc_tcga_pan_can_atlas_2018. All datasets of the Cell Model Passports, including ‘Model Annotation’, ‘Gene Annotation’, ‘Mutation Data’, and ‘Expression Data’ (version: ‘_20191101.zip’) are available at https://cellmodelpassports.sanger.ac.uk/downloads. All codes for reproducing our analyses shown in the main figures are available at https://github.com/namheee/rEMT.

### Logical regulatory network modeling of EMT

The molecular mechanisms that drive the EMT process are involved with multiple EMT-TFs, and their interactions become more complicated after TGFβ exposure (16), a primary inducer of EMT. Thus, the mechanism of EMT or rEMT is challenging to decipher intuitively. For this, we applied a logical Boolean network modeling approach to systematically explore the complex signaling dynamics of a wide spectrum of EMT, including hybrid states, and identify optimal molecular targets that effectively enable rEMT. This approach has been used to investigate biological processes with complex signaling networks (25). Through performing comparison analyses of nodes and pathways across epithelial and mesenchymal cell lines, we could select the major signaling pathways for EMT (Materials and Methods). With this, we identified key regulators of EMT by investigating the publicly available databases STRING (18), SIGNOR (26), and Kyoto Encyclopedia of Genes and Genomes (19) and also integrated all relevant relationships between the molecules from previous studies (8–10, 20, 27). The resulting EMT model consists of 31 nodes connected by 92 edges and cross-talk between the TGFβ (16), Wnt (28), Notch (29), MAPK (30), PI3K-Akt (16), NF-κB (16), and p53 (31) signaling pathways (Fig. 1A). The input node of the network is TGFβ and the output nodes are E-cadherin (E-cad) and ZEB1. More detailed information for constructing logical regulatory network is described in Materials and Methods.

Figure 1.

Construction of the EMT regulatory network model. A, Regulatory network for TGFβ-driven EMT. The network consists of 31 nodes and 92 edges. Each node is colored according to its related signaling pathway (Supplementary Table S1). Sixty-three edges are activation (red arrows) and 29 edges are inhibition (blue blunted arrows). B, Qualitative input–output relationships in the EMT network model corresponding experimental validation.

Figure 1.

Construction of the EMT regulatory network model. A, Regulatory network for TGFβ-driven EMT. The network consists of 31 nodes and 92 edges. Each node is colored according to its related signaling pathway (Supplementary Table S1). Sixty-three edges are activation (red arrows) and 29 edges are inhibition (blue blunted arrows). B, Qualitative input–output relationships in the EMT network model corresponding experimental validation.

Close modal

To elucidate biological dynamics underlying the EMT Boolean model, we qualitatively investigated individual input–output relationship and compared the results with quantitative transcript analysis from in vitro results. We showed that a lung cancer cell line (A549) displays hybrid state with increased levels of E-cad and ZEB1. After exogenous TGFβ exposure, A549 cells showed complete EMT state (A549-T, mesenchymal state) with decreased E-cad and increased ZEB1 (Supplementary Fig. S2A). We also assessed changes in transcript abundance associated with the epithelial and mesenchymal phenotypes in cells, and found consistent results with the protein abundance results (Supplementary Fig. S2B). We represented the mutation profile of A549 cells by setting the value of Ras to the ON state for all simulations (Materials and Methods) and confirmed that the qualitative effect of TGFβ within the network model is consistent with the published literature and our experiments (Fig. 1B). TGFβ exhibited a positive relationship with the EMT-TFs Snail (32), Twist1 (33), and ZEB1 (33); TGFβ promoted ERK phosphorylation (16) and c-Myc activity (34) by enhancing their activities through the MAPK pathway. Moreover, TGFβ decreased the activities of E-cad and EpCAM (35). Simulated results of steady states from other molecules were also in accordance with published data (Supplementary Fig. S2C; Supplementary Table S3). Thus, our EMT model reproduced experimentally observed behaviors of cancer cells and their EMT processes that depend on TGFβ.

### Limitations of rEMT by p53 activation

Advanced cancers include heterogeneous EMT states with hybrid and mesenchymal states after exposure to a TME signal, TGFβ (14). We thereby aim to find effective targets (indicative of a potential drug target) that can achieve epithelial states regardless of any cancer population (epithelial, hybrid, mesenchymal, or mixed states). To identify potential targets that can effectively revert EMT in the TGFβ-OFF and -ON conditions, we performed simulations with either the activity of each node set to OFF or ON and analyzed attractor landscapes (Materials and Methods). In addition, using the E-cad and ZEB1 nodes as phenotypic markers, we succinctly defined three states of the model as cellular phenotypes: epithelial (E-cad+/ZEB1), hybrid (E-cad+/ZEB1+, E-cad/ZEB1), and mesenchymal (E-cad/ZEB1+). With this perturbation analysis, we systematically screened potential targets using a quantitative phenotype score, the rEMT score, based on cosine-similarity between simulation results and desired states of the selected nine nodes that are functionally associated with the epithelial feature (E-cad, miR-200, miR-34, EpCAM) and the mesenchymal and stemness features (ZEB1, Snail, Twist1, c-Myc, Thy-1; Materials and Methods). In the absence of TGFβ, OE of p53 or KO of Mdm2 led to the epithelial phenotype with a high rEMT score of 0.86. However, in the presence of TGFβ, no single-node perturbation reached the epithelial phenotype, and p53 OE did not show effective rEMT with a score of 0.13 (Fig. 2A). Furthermore, we found that, in the absence of TGFβ, p53 OE in A549 cells can induce the epithelial phenotype (E-cad+/ZEB1) with some hybrid properties remaining, which we called the epithelial-like hybrid state (Epi/H). However, in the presence of TGFβ, p53 OE in A549-T cells did not restore E-cad activity but rather arrived in another state with the hybrid phenotype (E-cad/ZEB1; Fig. 2B; Supplementary Fig. S3A and Supplementary Data S2).

Figure 2.

Limitations of rEMT by p53 activation. A, Single-node perturbation analysis for each TGFβ condition. All simulations were plotted in bar plots representing rEMT score per each perturbation. Perturbation types are denoted by (-) and (+), representing KO and OE, respectively. B, Schematic description of rEMT after single-node perturbations under TGFβ-OFF (blue arrow) or -ON (red arrow) condition. The circle is colored according to its phenotype depending on E-cad and ZEB1 activities: epithelial (E, purple), hybrid (H, orange), mesenchymal (M, red). C and D, A549 cells were treated without or with Nutlin-3a (5 μmol/L) for 48 hours in the absence or presence of TGFβ (5 ng/mL). Western blots for E-cad, ZEB1, and p53 are shown with GAPDH as a loading control (C). Transcript analysis by qRT-PCR (n = 3, mean ± SD; D). E, Frustration of network state for each perturbation under TGFβ-OFF or -ON condition. All simulations were plotted in PCA plots in the x-y axis. The z-axis corresponds to frustration (Materials and Methods). Each data point is colored by its phenotype according to E-cad and ZEB1 activities: epithelial (E, purple), hybrid (H, orange), mesenchymal (M, red). Each plane represents a best-fit plane for three-dimensional simulation data using linear regression. F, Relative molecular state ambiguity in the presence of TGFβ. The top 20 nodes are identified. Each bar displays ΔΔrin or ΔΔrout, which is the difference of Δrin(xj) or Δrout(xj), for each node j, at the p53-perturbed states between TGFβ-OFF and -ON conditions (Materials and Methods).

Figure 2.

Limitations of rEMT by p53 activation. A, Single-node perturbation analysis for each TGFβ condition. All simulations were plotted in bar plots representing rEMT score per each perturbation. Perturbation types are denoted by (-) and (+), representing KO and OE, respectively. B, Schematic description of rEMT after single-node perturbations under TGFβ-OFF (blue arrow) or -ON (red arrow) condition. The circle is colored according to its phenotype depending on E-cad and ZEB1 activities: epithelial (E, purple), hybrid (H, orange), mesenchymal (M, red). C and D, A549 cells were treated without or with Nutlin-3a (5 μmol/L) for 48 hours in the absence or presence of TGFβ (5 ng/mL). Western blots for E-cad, ZEB1, and p53 are shown with GAPDH as a loading control (C). Transcript analysis by qRT-PCR (n = 3, mean ± SD; D). E, Frustration of network state for each perturbation under TGFβ-OFF or -ON condition. All simulations were plotted in PCA plots in the x-y axis. The z-axis corresponds to frustration (Materials and Methods). Each data point is colored by its phenotype according to E-cad and ZEB1 activities: epithelial (E, purple), hybrid (H, orange), mesenchymal (M, red). Each plane represents a best-fit plane for three-dimensional simulation data using linear regression. F, Relative molecular state ambiguity in the presence of TGFβ. The top 20 nodes are identified. Each bar displays ΔΔrin or ΔΔrout, which is the difference of Δrin(xj) or Δrout(xj), for each node j, at the p53-perturbed states between TGFβ-OFF and -ON conditions (Materials and Methods).

Close modal

To validate these results, we increased p53 activity through a MDM2 inhibitor, Nutlin-3a, in A549 cells and compared the abundance of proteins associated with the epithelial or mesenchymal phenotypes in cells without or with exogenous TGFβ exposure (Fig. 2C). Consistent with the previous studies (31, 36), we found that, in the absence of TGFβ, upregulation of p53 activity increased E-cad and decreased ZEB1. However, in the presence of TGFβ, upregulation of p53 activity only decreased ZEB1 abundance without any change in E-cad abundance. We also assessed changes in transcript abundance for CDH1 (encoding E-cad) and three mesenchymal markers: ZEB1, SNAI1 (encoding Snail), and TWIST1 (encoding Twist1). Consistent with the protein abundance results, we found that activating p53 resulted in an Epi/H state in the absence of TGFβ, but in a hybrid state in the presence of TGFβ (Fig. 2D). Similarly, SMAD4 KO, which is most effective for rEMT in the presence of TGFβ, still represented the hybrid phenotype (E-cad+/ZEB1+) from the in vitro experiments (Supplementary Fig. S3B). Finally, these results showed that any single-node perturbation is not sufficient to completely revert EMT, and despite the same perturbations, different responses are observed in the TGFβ-OFF and -ON conditions.

To examine what the difference is between the network dynamics depending on the input condition, we explored which regulatory relationships are mainly involved in determining hybrid states in the presence of TGFβ. For this, we first computed the relative stability of attractors based on network frustration (Materials and Methods; refs. 23, 37). In addition, we performed principal component analysis (PCA) using the average node activity profiles and showed that ambiguous cellular states, hybrid states, are detected in the region with high frustration, which is consistent with the previous studies (Fig. 2E; refs. 23, 37). In the absence of TGFβ, p53 OE induced the Epi/H state with lower frustration, whereas p53 OE induced the hybrid state with higher frustration in the presence of TGFβ. Furthermore, to identify key genes that give rise to more ambiguous relationships, resulting in higher frustration, we defined the molecular state ambiguity (Δr) and analyzed relative molecular state ambiguity (ΔΔr) by computing the difference in Δr between the TGFβ-OFF and -ON conditions (Fig. 2F, Materials and Methods). From this analysis, we hypothesized that top-ranked genes with higher molecular state ambiguities in the presence of TGFβ would be associated with determining hybrid states even after p53 activation.

### Identification of SMAD4 inhibition to restore the epithelial characteristics

On the basis of the representation of regulatory connections and average node activity profiles, we extracted a subnetwork that is involved in phenotype transition (detailed explanation for this process in Supplementary Fig. S4A–S4C and Materials and Methods). By adding top-ranked genes in Fig. 2F and comparing attractors between TGFβ-OFF and -ON conditions, we identified two positive FBLs that are involved in determining hybrid states within the network model (Fig. 3A). In the absence of TGFβ, we showed that the average activity of some nodes resembled the epithelial state, whereas other nodes resembled a hybrid state (Fig. 3A, left). Because PI3K, NF-κB, and Snail were inactive, p53 OE induced high E-cad activity through a double-negative FBL between Snail and miR-34. Moreover, p53 OE reduced ZEB1 activity along with a path involved with miR-200 by a functional double-negative FBL between these two nodes. In the presence of TGFβ, however, activities of most nodes were associated with the mesenchymal phenotype (Fig. 3A, right). In this TGFβ-ON state, despite high miR-34 activity due to p53 OE, the double-negative FBL between Snail and miR-34 did not reduce the activity of Snail, because the upstream nodes of Snail such as NOTCH, NF-κB, SMAD4, and β-catenin signaling remained active due to TGFβ. However, the activity of ZEB1 was reduced by a double-negative FBL between miR-200 and ZEB1 after p53 OE. Consequently, the effect of p53 OE in the presence of TGFβ signal changed the state of ZEB1, but not of E-cad and other mesenchymal markers, which led to a hybrid state of E-cad/ZEB1. We further discussed the functions of such FBLs by comparing with other studies in Supplementary Text S4 and Supplementary Table S4. In addition, when comparing regulatory connections between TGFβ-OFF and -ON, this TGFβ-ON state showed more frustrated outgoing edges from these FBLs. Thus, we hypothesized that reducing ambiguous outgoing signals from two positive FBLs is critical for avoiding hybrid states. On the basis of this analysis, we attempted to find optimal targets for rEMT.

Figure 3.

Identifying a combination target for rEMT through network analysis. A, Different regulatory dynamics after p53 OE. Networks (top) represent the average node activity profiles (high, above 0.5; low, below 0.5). Edges that are unfrustrated are red (activation) or blue (inhibition); frustrated edges are gray. Each subnetwork (bottom) highlights key components involved in EMT transition. The brown box encapsulates a circuit consisting of top-ranked genes from Fig. 2F. The red up-arrow and blue down-arrow indicate active and inactive, respectively. All edges except the ‘FBL-forming edges’ are indicated by dashed lines. B, Relative stability of the circuit after double-node perturbation analysis in the presence of TGFβ. The z-axis corresponds to the sum of Δrout of the nodes in the circuit from Fig. 3A. C, Top-ranked combination targets for rEMT. Each bar displays rEMT score (purple) and Δrout of the circuit (orange). D, The effect of p53 activation and SMAD4 knockdown. In silico: subnetwork after p53 OE and SMAD4 KO. In vitro: A549 cells with the indicated knockdown constructs were treated without or with Nutlin-3a in the absence or presence of TGFβ for 48 hours. Protein abundance was monitored by Western blotting with antibodies recognizing the indicated proteins. α-Actinin was used as a loading control. Scramble shRNA expressing cells, shScr; SMAD4 shRNA-expressing cells, shS4.

Figure 3.

Identifying a combination target for rEMT through network analysis. A, Different regulatory dynamics after p53 OE. Networks (top) represent the average node activity profiles (high, above 0.5; low, below 0.5). Edges that are unfrustrated are red (activation) or blue (inhibition); frustrated edges are gray. Each subnetwork (bottom) highlights key components involved in EMT transition. The brown box encapsulates a circuit consisting of top-ranked genes from Fig. 2F. The red up-arrow and blue down-arrow indicate active and inactive, respectively. All edges except the ‘FBL-forming edges’ are indicated by dashed lines. B, Relative stability of the circuit after double-node perturbation analysis in the presence of TGFβ. The z-axis corresponds to the sum of Δrout of the nodes in the circuit from Fig. 3A. C, Top-ranked combination targets for rEMT. Each bar displays rEMT score (purple) and Δrout of the circuit (orange). D, The effect of p53 activation and SMAD4 knockdown. In silico: subnetwork after p53 OE and SMAD4 KO. In vitro: A549 cells with the indicated knockdown constructs were treated without or with Nutlin-3a in the absence or presence of TGFβ for 48 hours. Protein abundance was monitored by Western blotting with antibodies recognizing the indicated proteins. α-Actinin was used as a loading control. Scramble shRNA expressing cells, shScr; SMAD4 shRNA-expressing cells, shS4.

Close modal

We performed double-node perturbation analysis and also computed the sum of Δrout of the nodes in the two positive FBLs to confirm our hypothesis (Materials and Methods). From this, we first showed that most perturbations representing E-cad+/ZEB1 have lower Δrout of the FBLs (Fig. 3B). Then top-ranked results were ordered by rEMT score combined with Δrout of the FBLs (Fig. 3C; Supplementary Data S2). As expected, direct TGFβ KO is the most effective. Next, without inhibiting this exogenous signal, we predicted that SMAD4 KO along with p53 OE is most effective at restoring the epithelial phenotype in the presence of TGFβ (rEMT score is 0.76). The combined perturbation directly downregulated Snail and upregulated miR-34 activity, enabling increased E-cad activity and decreased ZEB1 activity through the double-negative FBL between Snail and miR-34 (Fig. 3D, top). Then we validated the effect of this combination through in vitro experiments. We made a stable SMAD4 knockdown cell line and increased the activity of p53 with Nutlin-3a. These cells had a higher abundance of E-cad and a lower abundance of ZEB1 than other cells in the absence or presence of TGFβ (Fig. 3D, bottom). We showed that SMAD4 inhibition with p53 activation can revert mesenchymal cells into epithelial cells represented by E-cad+/ZEB1 regardless of TGFβ input signal through in silico and in vitro experiments.

### Network analysis to achieve rEMT while avoiding hybrid characteristics

We hypothesized that if we had successfully induced epithelial phenotypes, the reversed cells would have recovered their chemosensitivity because EMT is known to be deeply related to drug resistance (38, 39). To examine whether the effect of this combination induces chemosensitive epithelial states as well, we activated p53 with Nutlin-3a and reduced the expression of SMAD4 in A549 cells. We then exposed the cells to various chemotherapeutic drugs in the absence or presence of TGFβ. Unexpectedly, we found that p53 activation and SMAD4 knockdown, individually or together, did not significantly affect cell viability in the absence or presence of TGFβ with no sensitivity to 5-FU or DOX (Fig. 4A) or to OX or CPT-11 (Supplementary Fig. S5A). On the basis of this result, we concluded that cells with the increased E-cad activity and decreased ZEB1 activity are not sufficient to characterize a chemosensitive rEMT state (Fig. 4B).

Figure 4.

Molecular markers defining the epithelial phenotype with chemosensitivity. A, The effect of p53 activation and SMAD4 knockdown on the sensitivity of cells to chemotherapy. Cell viability was analyzed by IncuCyte over time after indicated drug treatments (n = 3). B, Schematic description of rEMT after double-node perturbations for TGFβ-OFF (blue arrow) or -ON (red arrow) condition. C, Average node activity profiles after indicated perturbations for TGFβ-OFF (top) or -ON (bottom) condition. D, Average node activity profiles of all epithelial states in the TGFβ-OFF (white) or -ON (gray) condition. E, The effect of p53 activation and SMAD4 knockdown on gene expression in the presence of TGFβ. Data are presented as mean ± SD. ns, nonsignificant, P > 0.05; ***, P ≤ 0.001; ****, P ≤ 0.0001 (Student t test). Scr or shScr, Scramble shRNA; shS4, SMAD4 shRNA; Veh, Vehicle.

Figure 4.

Molecular markers defining the epithelial phenotype with chemosensitivity. A, The effect of p53 activation and SMAD4 knockdown on the sensitivity of cells to chemotherapy. Cell viability was analyzed by IncuCyte over time after indicated drug treatments (n = 3). B, Schematic description of rEMT after double-node perturbations for TGFβ-OFF (blue arrow) or -ON (red arrow) condition. C, Average node activity profiles after indicated perturbations for TGFβ-OFF (top) or -ON (bottom) condition. D, Average node activity profiles of all epithelial states in the TGFβ-OFF (white) or -ON (gray) condition. E, The effect of p53 activation and SMAD4 knockdown on gene expression in the presence of TGFβ. Data are presented as mean ± SD. ns, nonsignificant, P > 0.05; ***, P ≤ 0.001; ****, P ≤ 0.0001 (Student t test). Scr or shScr, Scramble shRNA; shS4, SMAD4 shRNA; Veh, Vehicle.

Close modal

To identify molecular markers that can define rEMT with chemosensitivity, we first examined the activities of nine nodes representing rEMT score through p53 OE and/or SMAD4 KO simulations. In the SMAD4 KO with p53 OE, high Twist1 activity and reinforced EpCAM activity were observed in the simulation results for all TGFβ conditions (Fig. 4C). Furthermore, we investigated if high activities of Twist1 and EpCAM are also observed in all other epithelial states. Consistent with Fig. 4C, Twist1 and EpCAM had high activities regardless of TGFβ condition (Fig. 4D). Twist1 is one of the EMT-TFs that is a molecular marker of drug resistance in several types of cancers (40). Although EpCAM is an epithelial cell adhesion molecule, EpCAM-high cells have shown stem-cell properties (41), and downregulation of EpCAM have promoted chemosensitizing effect in many cancers (42, 43). On the basis of these, our results indicate that Twist1 and EpCAM would be potential molecular markers for drug resistance that results from the emergence of cancer stem cell–like properties in the epithelial group.

We further assessed whether Twist1 and/or EpCAM are related to chemotherapy responsiveness of p53 wild-type cells using public data, including TCGA lung patient samples (44) and Cell Model Passports (45). By analyzing various mesenchymal and stemness markers (EPCAM, MYC, THY1, TWIST1, and SNAI1), we confirmed that, except for EPCAM and TWIST1, there is no statistically significant difference in the survival rate of patients with wild-type p53 among the epithelial group (Supplementary Fig. S5B). Moreover, patients with wild-type p53 in the epithelial group and low TWIST1 and EPCAM had a significantly higher survival rate than those with high TWIST1 and EPCAM (Supplementary Fig. S5C). Similarly, in p53 wild-type cell lines, the epithelial group with high expression of both TWIST1 and EPCAM showed a significantly higher IC50 z-score to chemotherapy (Supplementary Fig. S5D, detailed results can be found in the supplementary information). Finally, we confirmed these molecular signatures with A549 cells. We measured mRNA expression levels of EMT-associated genes and stemness markers. Consistent with the simulation data in Fig. 4C, despite the reduced expression levels of ZEB1, VIM, and CDH2, TWIST1 expression remained high and EPCAM expression was increased after knockdown of SMAD4 in the presence of TGFβ (Fig. 4E; Supplementary Fig. S5E). Together, our in silico, bioinformatics, and in vitro results suggest that reducing Twist1 and EpCAM would be essential to eliminate stem cell–like properties of cancer cells in the epithelial group to achieve a complete rEMT state with chemosensitivity.

### Identification of ERK inhibition to restore the epithelial phenotype and overcome drug resistance

By performing triple-node perturbation analysis, we defined two epithelial groups using additional markers: Epi/H state represents Epi/H with high Twist1 and EpCAM activities, and epithelial state represents complete rEMT state with low Twist1 and EpCAM activities. As expected, regardless of TGFβ signal, the epithelial state showed the desired rEMT state in all markers including E-cad, ZEB1, Twist1, and EpCAM (Fig. 5A). To gain a mechanistic understanding of how cells stabilize in the epithelial state without any hybrid properties, we analyzed the difference in Δrout or Δrin of all nodes between Epi/H and epithelial states (Materials and Methods). We found that some molecules including Snail, miR-200, Twist1, and SMAD4 have more ambiguous relationships in the Epi/H state (Fig. 5B). We ranked the relative molecular state ambiguity of each FBL and identified a core-regulatory circuit (core-circuit) composed of top-ranked FBLs, which exhibit high molecular state ambiguities in the Epi/H state (Fig. 5C). With two FBLs in Fig. 3A, we found a three-node positive FBL between Snail, Twist, and AP-1 as well as another two-node positive FBL between SMAD4 and ERK, which are related to p38 signaling pathway (46). We discussed the known functions of these FBLs by comparing with other studies in Supplementary Text S4. Finally, we hypothesized that ambiguous outgoing signals from this core-circuit are the primary source of lingering hybrid characteristics in the epithelial state.

Figure 5.

Identifying ERK inhibition in combination with p53 activation and SMAD4 inhibition to achieve a chemosensitive rEMT. A, Comparison with Epi/H (green) and epithelial (light green) states in the epithelial group for TGFβ-OFF (top) or -ON (bottom) condition. Each bar plot illustrates average activities of the indicated marker genes. Each error bar represents the SD within Epi/H and epithelial state. B, Relative molecular state ambiguity in the Epi/H state for TGFβ-OFF (left) or -ON (right) condition. Each bar displays ΔΔr representing the sum of ΔΔrin(xj) (light gray) and ΔΔrout(xj) (dark gray) of each node, which is computed by the difference of Δrin or Δrout between Epi/H and epithelial states. C, FBLs with higher molecular state ambiguities in the Epi/H state for TGFβ-OFF (blue) or -ON (red) condition. Top-ranked FBLs are identified as potential core-circuit that could have the most significant role in reprogramming a mesenchymal state to a chemosensitive rEMT state. D, Relative stability of the core-circuit for each perturbation under TGFβ-OFF (top) or -ON (bottom) condition. The z-axis corresponds to the sum of Δrout of the nodes in the core-circuit from the results in C. A dashed circle encapsulates each simulation result with p53 OE (black), SMAD4 KO along with p53 OE (black), and ERK KO along with p53 OE and SMAD4 KO (red). E, Perturbation analysis for rEMT. The results were ranked by chemosensitive rEMT score combined with Δrout of the core-circuit (orange). Each bar representing chemosensitive rEMT score is colored according to Epi/H (green) and epithelial (light green) states. F, The core-circuit for chemosensitive rEMT. The subnetwork highlights the key regulatory components responsible for reducing hybrid properties in the epithelial group. Colored arrows indicate the direction of each node activity. Frustrated edges are gray. Each shaded region encapsulates each FBL that forms the core-circuit. All edges except the ‘FBLs-forming edge’ are indicated by dashed lines. G, Attractor landscape analysis. Each landscape represents attractors based on calculated basin sizes, and the region is colored according to its phenotype: Epi, light green; Epi/H, green; H, orange; Mes, red. Epi, epithelial; H, hybrid; Mes, mesenchymal.

Figure 5.

Identifying ERK inhibition in combination with p53 activation and SMAD4 inhibition to achieve a chemosensitive rEMT. A, Comparison with Epi/H (green) and epithelial (light green) states in the epithelial group for TGFβ-OFF (top) or -ON (bottom) condition. Each bar plot illustrates average activities of the indicated marker genes. Each error bar represents the SD within Epi/H and epithelial state. B, Relative molecular state ambiguity in the Epi/H state for TGFβ-OFF (left) or -ON (right) condition. Each bar displays ΔΔr representing the sum of ΔΔrin(xj) (light gray) and ΔΔrout(xj) (dark gray) of each node, which is computed by the difference of Δrin or Δrout between Epi/H and epithelial states. C, FBLs with higher molecular state ambiguities in the Epi/H state for TGFβ-OFF (blue) or -ON (red) condition. Top-ranked FBLs are identified as potential core-circuit that could have the most significant role in reprogramming a mesenchymal state to a chemosensitive rEMT state. D, Relative stability of the core-circuit for each perturbation under TGFβ-OFF (top) or -ON (bottom) condition. The z-axis corresponds to the sum of Δrout of the nodes in the core-circuit from the results in C. A dashed circle encapsulates each simulation result with p53 OE (black), SMAD4 KO along with p53 OE (black), and ERK KO along with p53 OE and SMAD4 KO (red). E, Perturbation analysis for rEMT. The results were ranked by chemosensitive rEMT score combined with Δrout of the core-circuit (orange). Each bar representing chemosensitive rEMT score is colored according to Epi/H (green) and epithelial (light green) states. F, The core-circuit for chemosensitive rEMT. The subnetwork highlights the key regulatory components responsible for reducing hybrid properties in the epithelial group. Colored arrows indicate the direction of each node activity. Frustrated edges are gray. Each shaded region encapsulates each FBL that forms the core-circuit. All edges except the ‘FBLs-forming edge’ are indicated by dashed lines. G, Attractor landscape analysis. Each landscape represents attractors based on calculated basin sizes, and the region is colored according to its phenotype: Epi, light green; Epi/H, green; H, orange; Mes, red. Epi, epithelial; H, hybrid; Mes, mesenchymal.

Close modal

To systematically examine the most effective target, we redefined the chemosensitive rEMT state as cells with high E-cad, miR-200, miR-34, and low ZEB1, Snail, Twist1, EpCAM, c-Myc, and Thy-1. Then top-ranked results were ordered using chemosensitive rEMT scores combined with the sum of the Δrout of nodes in the core-circuit (Materials and Methods). We confirmed that the relative stability of the core-circuit in the epithelial state is lower than in the Epi/H state regardless of the TGFβ condition (Fig. 5D). Finally, we found that adding ERK inhibition, through KO of the ERK node or the upstream MEK node, to p53 OE and SMAD4 KO is the most effective combination to achieve chemosensitive rEMT state with low relative stability of the core-circuit (Fig. 5E).

### Further investigation for transition to a chemosensitive rEMT state

We further analyzed this epithelial state resulting from ERK KO and indicated that the core-circuit, comprised of these four positive FBLs, contributes to transition to the epithelial state. Adding ERK KO reduced Snail, Twist1, and AP-1 activities, contributing to decreased EpCAM and c-Myc activities and reduced the number of frustrated outgoing edges in the network (Fig. 5F). We also revealed that perturbations modify the distribution of attractors, including epithelial, Epi/H, hybrid, and mesenchymal phenotypes. In the presence of TGFβ, after p53 OE, only the basin of hybrid phenotype was increased to 63%. Using SMAD4 KO along with p53 OE, the basin of the Epi/H phenotype was increased to 85%. Finally, after ERK KO along with p53 OE and SMAD4 KO, the basin of the epithelial phenotype was increased to 93%, and all five attractors, including point and cyclic attractors, showed high epithelial feature (the average epithelial node activity is 0.96) and low mesenchymal and stemness features (the average mesenchymal and stemness node activity is 0.04; Fig. 5G; Supplementary Data S2). On the basis of these network analyses, we hypothesized that combined inhibition of SMAD4 and ERK along with p53 activation are capable of driving mesenchymal cells to a chemosensitive rEMT state (epithelial phenotype).

### Validation of combined inhibition of SMAD4 and ERKs along with p53 activation

To validate our simulation results, we activated p53 using Nutlin-3a, blocked SMAD4 using short hairpin RNA (shRNA), and inhibited ERKs using a pharmacologic inhibitor of MEK or shRNA targeting ERKs. Activation of p53 and simultaneous inhibition of SMAD4 and ERKs increased the abundance of E-cad and decreased the abundance of ZEB1 in the absence or presence of TGFβ (Fig. 6A). We also examined the effect of p53 activation on cells with knockdown of SMAD4 and ERK1/2. In these cells, p53 activation also increased the abundance of E-cad and decreased the abundance of ZEB1 even in the presence of TGFβ (Fig. 6B). To demonstrate that the combination of p53 activation with SMAD4 and ERK1/2 knockdown induced epithelial state–like expression profile, we examined the expression of stemness and mesenchymal genes in the presence of TGFβ. The expression of EPCAM and TWIST1 persisted in cells with p53 activation and SMAD4 knockdown, but p53 activation with SMAD4 and ERK1/2 knockdown impaired the expressions of the mesenchymal markers SNAI1, TWIST1, ZEB1, VIM, and CDH2 and the stemness markers EPCAM, MYC, and THY1 (Fig. 6C; Supplementary Fig. S6A).

Figure 6.

Combination effect of SMAD4 and ERK inhibition with p53 activation in achieving chemosensitive rEMT. A and B, Western blots for the abundance of the indicated proteins in A549 cells. GAPDH was used as a loading control. p53 activation and SMAD4 knockdown cells were treated with MEK inhibitor (U0126) in the absence or presence of TGFβ for 48 hours (A). p53 was activated with Nutlin-3a in cells with SMAD4 and ERK knockdown in the absence or presence of TGFβ for 48 hours (B). C, The effect of p53 activation and SMAD4 and ERK knockdown on gene expression in the presence of TGFβ. Transcript analysis by qRT-PCR is presented relative to that in the scramble shRNA condition. D, Sensitivity to chemotherapy of cells with p53 activation and SMAD4 and ERK knockdown. A549 cells with the indicated knockdown constructs were treated with the indicated drugs for 5 days in the absence or presence of TGFβ. Cell viability was analyzed by IncuCyte. Representative images of crystal violet staining are shown below. Data are presented as mean ± SD of three biological replicates. ns, nonsignificant, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001 (Student t test). Scr or shScr, scramble shRNA; shS4, SMAD4 shRNA; shERK1/2, shERKs; MEK_i, MEK inhibitor; N3a, Nutlin-3a; Veh, Vehicle.

Figure 6.

Combination effect of SMAD4 and ERK inhibition with p53 activation in achieving chemosensitive rEMT. A and B, Western blots for the abundance of the indicated proteins in A549 cells. GAPDH was used as a loading control. p53 activation and SMAD4 knockdown cells were treated with MEK inhibitor (U0126) in the absence or presence of TGFβ for 48 hours (A). p53 was activated with Nutlin-3a in cells with SMAD4 and ERK knockdown in the absence or presence of TGFβ for 48 hours (B). C, The effect of p53 activation and SMAD4 and ERK knockdown on gene expression in the presence of TGFβ. Transcript analysis by qRT-PCR is presented relative to that in the scramble shRNA condition. D, Sensitivity to chemotherapy of cells with p53 activation and SMAD4 and ERK knockdown. A549 cells with the indicated knockdown constructs were treated with the indicated drugs for 5 days in the absence or presence of TGFβ. Cell viability was analyzed by IncuCyte. Representative images of crystal violet staining are shown below. Data are presented as mean ± SD of three biological replicates. ns, nonsignificant, P > 0.05; *, P ≤ 0.05; **, P ≤ 0.01; ***, P ≤ 0.001 (Student t test). Scr or shScr, scramble shRNA; shS4, SMAD4 shRNA; shERK1/2, shERKs; MEK_i, MEK inhibitor; N3a, Nutlin-3a; Veh, Vehicle.

Close modal

To validate that such reversed cells can overcome drug resistance as expected, we measured cell viability after treating the cells with 5-FU or DOX in the absence or presence of TGFβ. We confirmed that p53 activation and SMAD4 and ERK1/2 knockdown cells have reduced cell viability compared to other cells in the absence or presence of TGFβ when treated with 5-FU or DOX (Fig. 6D) or with OX or CPT-11 (Supplementary Fig. S6B). We further conducted perturbation analysis with other lung cancer cell line, A427, by reflecting its mutation profile and found that our combination target is also top-ranked in A427 as in A549 (Supplementary Fig. S7A and S7B). Consistent with the results from A549 cells, our combination target increased the abundance of E-cad and decreased the abundance of ZEB1 in A427 cells, and such reversed cells could overcome drug resistance as well (Supplementary Fig. S7C and S7D). These results showed that blocking SMAD4 and ERKs along with p53 activation can revert EMT fully without any hybrid characteristic and finally overcome chemoresistance in the presence of an EMT-promoting TME signal.

rEMT is considered as a promising approach for cancer therapy, thus multiple attempts have been made to revert EMT (31, 47). However, most studies have been unsuccessful at least in part because they failed to systematically analyze the complexity of the process. EMT and rEMT are not binary state transitions between fully epithelial and mesenchymal states, but rather involve frequent transitions between multiple hybrid states during the processes (12). This plasticity leads cells to acquire stemness characteristics, which can cause chemotherapy resistance (48, 49). Here, we applied a systems biological approach of logical regulatory network analysis along with cell experiments to revert EMT within cancer cells while avoiding multiple hybrid states (Fig. 7A). By employing a Boolean network modeling framework called coarse-grained regulatory network modeling, which can cover larger molecular networks without detailed estimation of kinetic parameters (25), we unraveled highly complex and plastic EMT and rEMT processes and discovered possible therapeutic interventions for achieving rEMT within the network model.

Figure 7.

A schematic representation of rEMT while avoiding hybrid states. A, A schematic view of achieving rEMT. This conceptual model summarizes our study; cells representing mesenchymal phenotypes can revert to cells with epithelial phenotypes while avoiding multiple hybrid states and finally respond to chemotherapeutic agents (top). Gray dashed arrow indicates the EMT induction through TGFβ signal secreted from the cells in the TME. Each colored arrow corresponds to each perturbation, p53 activation (red), SMAD4 inhibition (dark yellow), and ERK inhibition (purple), respectively. The bottom panel represents the epithelial (purple), the mesenchymal (red), and stemness (dark yellow) features of each state during the rEMT process. B, The core-circuit that plays the primary role in driving rEMT without undesired hybrid states. p53 activation did not induce E-cad activation and resulted in the hybrid phenotype (left) due to a malfunction of FBL between Snail and miR-34 (light brown). The SMAD4 inhibition with p53 activation induced E-cad activation, but still induced Twist1 and EpCAM activities (middle) through two other malfunctioning FBLs (light brown). Combined inhibition of SMAD4 and ERKs along with p53 activation finally induced a chemosensitive rEMT while overcoming cellular heterogeneity and also restoring drug sensitivity (right). All edges except the ‘FBL-forming edges’ are indicated by dashed lines.

Figure 7.

A schematic representation of rEMT while avoiding hybrid states. A, A schematic view of achieving rEMT. This conceptual model summarizes our study; cells representing mesenchymal phenotypes can revert to cells with epithelial phenotypes while avoiding multiple hybrid states and finally respond to chemotherapeutic agents (top). Gray dashed arrow indicates the EMT induction through TGFβ signal secreted from the cells in the TME. Each colored arrow corresponds to each perturbation, p53 activation (red), SMAD4 inhibition (dark yellow), and ERK inhibition (purple), respectively. The bottom panel represents the epithelial (purple), the mesenchymal (red), and stemness (dark yellow) features of each state during the rEMT process. B, The core-circuit that plays the primary role in driving rEMT without undesired hybrid states. p53 activation did not induce E-cad activation and resulted in the hybrid phenotype (left) due to a malfunction of FBL between Snail and miR-34 (light brown). The SMAD4 inhibition with p53 activation induced E-cad activation, but still induced Twist1 and EpCAM activities (middle) through two other malfunctioning FBLs (light brown). Combined inhibition of SMAD4 and ERKs along with p53 activation finally induced a chemosensitive rEMT while overcoming cellular heterogeneity and also restoring drug sensitivity (right). All edges except the ‘FBL-forming edges’ are indicated by dashed lines.

Close modal

We found that p53 activation induces a hybrid state through a loss of function of the double-negative FBL between Snail and miR-34 within the cells that acquired the mesenchymal phenotype from TGFβ stimulation (Fig. 7B, left). After inhibiting SMAD4 along with p53 activation, cells were able to escape from hybrid states by restoring the above FBL with reduced Snail activity and increased miR-34 activity, and it ultimately induced some node activities to those of the epithelial state. However, the rEMT state achieved by targeting SMAD4 and p53 still had hybrid properties that are deeply associated with drug resistance due to high Twist1 and EpCAM activities (Fig. 7B, middle). We discovered that a functional positive FBL between Snail, Twist, and AP-1, as well as another functional positive FBL between SMAD4 and ERK are necessary to drive cancer cells to a chemosensitive rEMT state with fully suppressed mesenchymal and stem-cell properties (Fig. 7B, right). As a result, combined inhibition of SMAD4 and ERKs with p53 activation was required to inhibit activities of all nodes involved in these FBLs. From a systems perspective using computer simulation combined with cell and molecular experiment, we finally unraveled that this functional core-circuit is crucial in inducing cell-fate transition to the chemosensitive rEMT state.

rEMT can suppress early stages of metastasis and force cells to maintain a controllable tumorigenic state without any invasive property at the primary site. Although EMT can be triggered by various conditions, including growth factors, inflammatory cytokines, and hypoxia, the scope of this study was limited to investigating only the primary condition of the EMT-promoting TME signal, TGFβ. Hence, it is necessary to expand our network model to further investigate the role of various other TME signals. Moreover, the clinical application of rEMT should be carefully evaluated because inducing rEMT after metastasis or during metastasis may contribute to the establishment of tumors at secondary sites. Another limitation of this study is that the presented approach is only applicable to p53 wild-type tumors and would not be able to induce rEMT of p53 KO tumors.

A key contribution of this study is the presentation of a quantitative framework for cell-fate reprogramming that avoids undesired or unexpected states even with an external physiological signal. We found that regulation of p53, SMAD4, and ERKs in combination may be used as a new treatment option for patients with primary cancer with metastatic and stemness phenotypes or chemoresistant cancer. Our attractor landscape-based reprogramming study can be further applied to identifying hidden regulatory mechanisms or therapeutic targets to overcome cellular heterogeneity and induce reversal of cancer cells to normal-like cells without undesired states for clinical application.

K.H. Cho reports patents for 10-2020-0187110 and PCT/KR2021/020161 pending, licensed, and with royalties paid from biorevert, Inc. No disclosures were reported by the other authors.

N. Kim: Resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft. C.Y. Hwang: Resources, validation, investigation, visualization, methodology, writing–original draft. T. Kim: Investigation, methodology. H. Kim: Resources, data curation, investigation. K.H. Cho: Conceptualization, formal analysis, supervision, funding acquisition, validation, investigation, project administration, writing–review and editing.

This work was supported by the National Research Foundation of Korea (NRF) funded by the Korea Government, the Ministry of Science and ICT (2020R1A2B5B03094920), and by Electronics and Telecommunications Research Institute (ETRI) grants (22ZS1100, Core Technology Research for Self-Improving Integrated Artificial Intelligence System). It was also supported by the Bio & Medical Technology Development Program of NRF funded by the Ministry of Science & ICT (2021M3A9I4024447) and the Korea Health Technology R&D Project through the Korea Health Industry Development Institute (KHIDI) and Korea Dementia Research Center (KDRC), funded by the Ministry of Health & Welfare and Ministry of Science and ICT, Republic of Korea (grant number: HU21C0060). The authors thank Nancy R. Gough for editorial service and Jonghoon Lee, Jae Il Joo, Sea Rom Choi, and Corbin Hopper for their critical reading and comments.

The publication costs of this article were defrayed in part by the payment of publication fees. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Note: Supplementary data for this article are available at Cancer Research Online (http://cancerres.aacrjournals.org/).

1.
Sung
H
,
Ferlay
J
,
Siegel
RL
,
Laversanne
M
,
Soerjomataram
I
,
Jemal
A
, et al
.
Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries
.
CA Cancer J Clin
2021
;
71
:
209
49
.
2.
van Zijl
F
,
Krupitza
G
,
Mikulits
W
.
Initial steps of metastasis: cell invasion and endothelial transmigration
.
Mutat Res
2011
;
728
:
23
34
.
3.
Nieto
MA
,
Huang
RY
,
Jackson
RA
,
Thiery
JP
.
EMT: 2016
.
Cell
2016
;
166
:
21
45
.
4.
Hay
ED
,
Zuk
A
.
Transformations between epithelium and mesenchyme: normal, pathological, and experimentally induced
.
Am J Kidney Dis
1995
;
26
:
678
90
.
5.
Ye
X
,
Weinberg
RA
.
Epithelial–mesenchymal plasticity: a central regulator of cancer progression
.
Trends Cell Biol
2015
;
25
:
675
86
.
6.
Pastushenko
I
,
Blanpain
C
.
EMT transition states during tumor progression and metastasis
.
Trends Cell Biol
2019
;
29
:
212
26
.
7.
Bierie
B
,
Pierce
SE
,
Kroeger
C
,
Stover
DG
,
Pattabiraman
DR
,
Thiru
P
, et al
.
Integrin-β4 identifies cancer stem cell–enriched populations of partially mesenchymal carcinoma cells
.
Proc Natl Acad Sci USA.
2017
;
114
:
E2337
E46
.
8.
Lu
M
,
Jolly
MK
,
Levine
H
,
Onuchic
JN
,
Ben-Jacob
E
.
MicroRNA-based regulation of epithelial-hybrid-mesenchymal fate determination
.
Proc Natl Acad Sci USA
2013
;
110
:
18144
9
.
9.
Tian
XJ
,
Zhang
H
,
Xing
J
.
Coupled reversible and irreversible bistable switches underlying TGFβ-induced epithelial to mesenchymal transition
.
Biophys J
2013
;
105
:
1079
89
.
10.
Zhang
J
,
Tian
XJ
,
Zhang
H
,
Teng
Y
,
Li
R
,
Bai
F
, et al
.
TGFβ-induced epithelial-to-mesenchymal transition proceeds through stepwise activation of multiple feedback loops
.
Sci Signal
2014
;
7
:
ra91
.
11.
Yu
M
,
Bardia
A
,
Wittner
BS
,
Stott
SL
,
Smas
ME
,
Ting
DT
, et al
.
Circulating breast tumor cells exhibit dynamic changes in epithelial and mesenchymal composition
.
Science
2013
;
339
:
580
4
.
12.
Udyavar
AR
,
Wooten
DJ
,
Hoeksema
M
,
Bansal
M
,
Califano
A
,
L
, et al
.
Correction: novel hybrid phenotype revealed in small cell lung cancer by a transcription factor network model that can explain tumor heterogeneity
.
Cancer Res
2019
;
79
:
1014
.
13.
Li
H
,
Xu
F
,
Li
S
,
Zhong
A
,
Meng
X
,
Lai
M
.
The tumor microenvironment: An irreplaceable element of tumor budding and epithelial–mesenchymal transition-mediated cancer metastasis
.
2016
;
10
:
434
46
.
14.
Tsai
JH
,
Yang
J
.
Epithelial–mesenchymal plasticity in carcinoma metastasis
.
Genes Dev
2013
;
27
:
2192
206
.
15.
Hao
Y
,
Baker
D
,
Ten Dijke
DP
.
TGFβ-mediated epithelial–mesenchymal transition and cancer metastasis
.
Int J Mol Sci
2019
;
20
:
2767
.
16.
Lamouille
S
,
Xu
J
,
Derynck
R
.
Molecular mechanisms of epithelial–mesenchymal transition
.
Nat Rev Mol Cell Biol
2014
;
15
:
178
96
.
17.
Shibue
T
,
Weinberg
RA
.
EMT, CSCs, and drug resistance: the mechanistic link and clinical implications
.
Nat Rev Clin Oncol
2017
;
14
:
611
29
.
18.
Szklarczyk
D
,
Franceschini
A
,
Wyder
S
,
Forslund
K
,
Heller
D
,
Huerta-Cepas
J
, et al
.
STRING v10: protein–protein interaction networks, integrated over the tree of life
.
Nucleic Acids Res
2015
;
43
:
D447
52
.
19.
Kanehisa
M
,
Goto
S
,
Sato
Y
,
Furumichi
M
,
Tanabe
M
.
KEGG for integration and interpretation of large-scale molecular data sets
.
Nucleic Acids Res
2012
;
40
:
D109
14
.
20.
Steinway
SN
,
Zanudo
JG
,
Ding
W
,
Rountree
CB
,
Feith
DJ
,
Loughran
TP
Jr
, et al
.
Network modeling of TGFβ signaling in hepatocellular carcinoma epithelial-to-mesenchymal transition reveals joint sonic hedgehog and Wnt pathway activation
.
Cancer Res
2014
;
74
:
5963
77
.
21.
Choi
M
,
Shi
J
,
Jung
SH
,
Chen
X
,
Cho
KH
.
Attractor landscape analysis reveals feedback loops in the p53 network that control the cellular response to DNA damage
.
Sci Signal
2012
;
5
:
ra83
.
22.
Choi
M
,
Shi
J
,
Zhu
Y
,
Yang
R
,
Cho
KH
.
Network dynamics–based cancer panel stratification for systemic prediction of anticancer drug response
.
Nat Commun
2017
;
8
:
1940
.
23.
Tripathi
S
,
Kessler
DA
,
Levine
H
.
Biological networks regulating cell-fate choice are minimally frustrated
.
Phys Rev Lett
2020
;
125
:
088101
.
24.
Shin
SY
,
Kim
T
,
Lee
HS
,
Kang
JH
,
Lee
JY
,
Cho
KH
, et al
.
The switching role of β-adrenergic receptor signaling in cell survival or death decision of cardiomyocytes
.
Nat Commun
2014
;
5
:
5777
.
25.
Wang
RS
,
A
,
Albert
R
.
Boolean modeling in systems biology: an overview of methodology and applications
.
Phys Biol
2012
;
9
:
055001
.
26.
Licata
L
,
Lo Surdo
P
,
Iannuccelli
M
,
Palma
A
,
Micarelli
E
,
Perfetto
L
, et al
.
SIGNOR 2.0, the signaling network open resource 2.0: 2019 update
.
Nucleic Acids Res
2020
;
48
:
D504
D10
.
27.
Silveira
DA
,
Mombach
JCM
.
Dynamics of the feedback loops required for the phenotypic stabilization in the epithelial–mesenchymal transition
.
FEBS J
2020
;
287
:
578
88
.
28.
Zhang
J
,
Tian
XJ
,
Xing
J
.
Signal transduction pathways of EMT induced by TGFbeta, SHH, and WNT and their cross-talks
.
J Clin Med
2016
;
5
:
41
.
29.
Yuan
X
,
Wu
H
,
Han
N
,
Xu
H
,
Chu
Q
,
Yu
S
, et al
.
Notch signaling and EMT in non–small cell lung cancer: biological significance and therapeutic application
.
J Hematol Oncol
2014
;
7
:
87
.
30.
Mulder
KM
.
Role of Ras and Mapks in TGFbeta signaling
.
Cytokine Growth Factor Rev
2000
;
11
:
23
35
.
31.
Chang
CJ
,
Chao
CH
,
Xia
W
,
Yang
JY
,
Xiong
Y
,
Li
CW
, et al
.
p53 regulates epithelial–mesenchymal transition and stem cell properties through modulating miRNAs
.
Nat Cell Biol
2011
;
13
:
317
23
.
32.
H
,
Quintanilla
M
,
Cano
A
.
Transforming growth factor beta-1 induces snail transcription factor in epithelial cell lines: mechanisms for epithelial–mesenchymal transitions
.
J Biol Chem
2003
;
278
:
21113
23
.
33.
Xu
J
,
Lamouille
S
,
Derynck
R
.
TGFβ-induced epithelial-to-mesenchymal transition
.
Cell Res
2009
;
19
:
156
72
.
34.
Smith
AP
,
Verrecchia
A
,
Faga
G
,
Doni
M
,
Perna
D
,
Martinato
F
, et al
.
A positive role for Myc in TGFbeta-induced Snail transcription and epithelial-to-mesenchymal transition
.
Oncogene
2009
;
28
:
422
30
.
35.
Eberlein
C
,
Rooney
C
,
Ross
SJ
,
Farren
M
,
Weir
HM
,
Barry
ST
.
E-cadherin and EpCAM expression by NSCLC tumor cells associate with normal fibroblast activation through a pathway initiated by integrin alphavbeta6 and maintained through TGFβ signaling
.
Oncogene
2015
;
34
:
704
16
.
36.
Powell
E
,
Piwnica-Worms
D
,
Piwnica-Worms
H
.
Contribution of p53 to metastasis
.
Cancer Discov
2014
;
4
:
405
14
.
37.
Font-Clos
F
,
Zapperi
S
,
La Porta
CAM
.
Topography of epithelial–mesenchymal plasticity
.
Proc Natl Acad Sci U S A
2018
;
115
:
5902
7
.
38.
Fischer
KR
,
Durrans
A
,
Lee
S
,
Sheng
J
,
Li
F
,
Wong
ST
, et al
.
Epithelial-to-mesenchymal transition is not required for lung metastasis but contributes to chemoresistance
.
Nature
2015
;
527
:
472
6
.
39.
Zheng
X
,
Carstens
JL
,
Kim
J
,
Scheible
M
,
Kaye
J
,
Sugimoto
H
, et al
.
Epithelial-to-mesenchymal transition is dispensable for metastasis but induces chemoresistance in pancreatic cancer
.
Nature
2015
;
527
:
525
30
.
40.
Deng
JJ
,
Zhang
W
,
Xu
XM
,
Zhang
F
,
Tao
WP
,
Ye
JJ
, et al
.
Twist mediates an aggressive phenotype in human colorectal cancer cells
.
Int J Oncol
2016
;
48
:
1117
24
.
41.
Imrich
S
,
Hachmeister
M
,
Gires
O
.
EpCAM and its potential role in tumor-initiating cells
.
2012
;
6
:
30
8
.
42.
Gao
J
,
Yan
Q
,
Liu
S
,
Yang
X
.
Knockdown of EpCAM enhances the chemosensitivity of breast cancer cells to 5-fluorouracil by downregulating the antiapoptotic factor Bcl-2
.
PLoS One
2014
;
9
:
e102590
.
43.
Richter
CE
,
Cocco
E
,
Bellone
S
,
Silasi
DA
,
Ruttinger
D
,
Azodi
M
, et al
.
High-grade, chemotherapy-resistant ovarian carcinomas overexpress epithelial cell adhesion molecule (EpCAM) and are highly sensitive to immunotherapy with MT201, a fully human monoclonal anti-EpCAM antibody
.
Am J Obstet Gynecol
2010
;
203
:
582
.
44.
KA
,
Yau
C
,
Hinoue
T
,
Wolf
DM
,
Lazar
AJ
,
Drill
E
, et al
.
Cell-of-origin patterns dominate the molecular classification of 10,000 tumors from 33 types of cancer
.
Cell
2018
;
173
:
291
304
.
45.
van der Meer
D
,
Barthorpe
S
,
Yang
W
,
Lightfoot
H
,
Hall
C
,
Gilbert
J
, et al
.
Cell model passports-a hub for clinical, genetic and functional datasets of preclinical cancer models
.
Nucleic Acids Res
2019
;
47
:
D923
D9
.
46.
Kudaravalli
S
,
den Hollander
P
,
Mani
SA
.
Role of p38 MAP kinase in cancer stem cells and metastasis
.
Oncogene
2022
;
41
:
3177
85
.
47.
Kurimoto
R
,
Ebata
T
,
Iwasawa
S
,
Ishiwata
T
,
Y
,
Tatsumi
K
, et al
.
Pirfenidone may revert the epithelial-to-mesenchymal transition in human lung adenocarcinoma
.
Oncol Lett
2017
;
14
:
944
50
.
48.
Grosse-Wilde
A
,
Fouquier d'Herouel
A
,
McIntosh
E
,
Ertaylan
G
,
Skupin
A
,
Kuestner
RE
, et al
.
Stemness of the hybrid epithelial/mesenchymal state in breast cancer and its association with poor survival
.
PLoS One
2015
;
10
:
e0126522
.
49.
Jolly
MK
,
Boareto
M
,
Huang
B
,
Jia
D
,
Lu
M
,
Ben-Jacob
E
, et al
.
Implications of the hybrid epithelial/mesenchymal phenotype in metastasis
.
Front Oncol
2015
;
5
:
155
.
This open access article is distributed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International (CC BY-NC-ND 4.0) license.