A Cell-Fate Reprogramming Strategy Reverses Epithelial-to-Mesenchymal Transition of Lung Cancer Cells While Avoiding Hybrid States

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.

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.

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

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.

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.

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 x_i 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, r_in(x_j), and outgoing edges, r_out(x_j), for node j (Supplementary Fig. S1C and S1D). Detailed methods for this simulation are in Supplementary Method S5.

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.

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

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

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

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% CO₂. Cells were periodically tested for Mycoplasma contamination (e-Myco plus Mycoplasma PCR Detection Kit, iNtRON Biotechnology).

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

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

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.

To analyze the effect of chemotherapeutic drugs on cell viability, transfected cells were seeded in 96-well plate (4 × 10³ 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.

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

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.

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.

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β.

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).

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.

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.

We performed double-node perturbation analysis and also computed the sum of Δr_out 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 Δr_out of the FBLs (Fig. 3B). Then top-ranked results were ordered by rEMT score combined with Δr_out 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.

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).

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 IC₅₀ 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.

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 Δr_out or Δr_in 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.

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 Δr_out 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).

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).

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).

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.

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.