Epithelial-to-mesenchymal transition (EMT) has been associated with cancer cell heterogeneity, plasticity, and metastasis. However, the extrinsic signals supervising these phenotypic transitions remain elusive. To assess how selected microenvironmental signals control cancer-associated phenotypes along the EMT continuum, we defined a logical model of the EMT cellular network that yields qualitative degrees of cell adhesions by adherens junctions and focal adhesions, two features affected during EMT. The model attractors recovered epithelial, mesenchymal, and hybrid phenotypes. Simulations showed that hybrid phenotypes may arise through independent molecular paths involving stringent extrinsic signals. Of particular interest, model predictions and their experimental validations indicated that: (i) stiffening of the extracellular matrix was a prerequisite for cells overactivating FAK_SRC to upregulate SNAIL and acquire a mesenchymal phenotype and (ii) FAK_SRC inhibition of cell–cell contacts through the receptor-type tyrosine-protein phosphatases kappa led to acquisition of a full mesenchymal, rather than a hybrid, phenotype. Altogether, these computational and experimental approaches allow assessment of critical microenvironmental signals controlling hybrid EMT phenotypes and indicate that EMT involves multiple molecular programs.
A multidisciplinary study sheds light on microenvironmental signals controlling cancer cell plasticity along EMT and suggests that hybrid and mesenchymal phenotypes arise through independent molecular paths.
Metastasis is a hallmark of cancer and the leading cause of mortality among patients with cancer. Despite intensive effort in basic and clinical research, metastatic cancers still present a major barrier to favorable clinical outcomes (1). Therefore, the fight against cancer calls for a better understanding of the involved cellular mechanisms.
The progression from carcinoma to metastatic cancer has been proposed to involve a shift from an epithelial to a mesenchymal phenotype, in a highly plastic and dynamic process referred to as epithelial-to-mesenchymal transition (EMT). During this transition, epithelial cells downregulate epithelial markers, lose their connections with neighboring cells through the breakdown of E-cadherin (ECad)-mediated adherens junction (AJ), upregulate mesenchymal markers, and acquire a marked migratory capacity mediated by the dynamic remodeling of focal adhesions (FA). Colonization at distant sites and metastatic outgrowth may require that disseminated cancer cells lose their migratory capacities and reestablish AJs through a mesenchymal-to-epithelial transition (MET). However, cancer cells are rarely purely mesenchymal or purely epithelial. They often exhibit both epithelial and mesenchymal features. Cells with such hybrid phenotypes appear to reside at intermediate states along the epithelial-to-mesenchymal continuum. Unlike mesenchymal phenotypes, these hybrid phenotypes may bear multiple advantages to cancer cells, including drug resistance and tumor-initiating potential (2).
The cancer-associated EMT program is hardly activated in a cell autonomous manner. Indeed, the tumor microenvironment induces EMT in carcinoma cells by releasing paracrine cell–cell signaling molecules, (WNT, the NOTCH ligand DELTA), growth factors (EGF and hepatocyte growth factor, HGF), and inflammatory signals (IL6; reactive oxygen species, ROS; and TGFβ; ref. 2). In addition, recent evidence suggests that direct physical interactions of tumor cells with their neighbors or with the extracellular matrix (ECM) greatly affect the EMT/MET program. In particular, ECM stiffening can induce EMT and promote tumor invasion and metastasis (3). In contrast, the cell adhesion molecule FAT4, which heterophilically interacts with its ligand on adjacent cells, prevents EMT and metastasis of gastric cancer cells (4). Other cell adhesion molecules, such as receptor-type tyrosine-protein phosphatases (RPTP) of the R2B subfamily, which display homophilic cell–cell adhesion capabilities, could also have major impact on the EMT/MET program. Indeed, they are regulators of AJs, are frequently mutated in solid cancers, and appear to display tumor suppressive capabilities (5). These microenvironmental signals cooperate to induce a group of EMT-inducing transcription factors (EMT-TF), including the zinc-finger proteins SNAIL (SNAIL and SLUG) and the E-box-binding protein ZEB1. These EMT-TFs regulate each other and, in different combinations, control the expression of genes associated with the epithelial and mesenchymal phenotypes, respectively (2). However, it remains unclear which combinations of external signals can stabilize carcinoma cells into hybrid phenotypes.
The complexity of the molecular networks involved in EMT has prompted numerous modeling studies, recently reviewed in ref. 6. Ordinary differential equation models focused on core regulatory circuitries and demonstrated the existence of stable, hybrid phenotypes (7, 8). In contrast to these continuous models, discrete logical models were developed that encompass numerous players and signaling pathways (9–12). Steinway and colleagues' model was built around the TGFβ pathway in the context of hepatocellular carcinoma, hybrid phenotypes being obtained through model perturbations (10). Other models considered the regulatory control of cell fate decisions between cell-cycle arrest, apoptosis, EMT, invasion and migration (11), or epithelial, mesenchymal, and senescent states (12). Overall, these studies demonstrated that, despite a coarse-grain representation, logical models are suitable to account for large and intricate networks, recovering physiologic cell phenotypes, and providing relevant predictions.
Here, we used computational modeling and experimental validations to investigate the role of selected microenvironmental signals on two features largely affected during EMT: cell–cell and cell–ECM adhesion properties (workflow in Supplementary Fig. S1). In addition to purely epithelial and mesenchymal phenotypes, our logical model recovers hybrid phenotypes, which could represent terminal phenotypical steps and whose stability depends on stringent inputs from the microenvironment. Among these, we have validated experimentally the role of cell–cell contacts through RPTP-kappa and of ECM stiffening. Altogether, our work identified key microenvironmental cues involved in maintaining hybrid EMT phenotypes that could be used as targets for developing therapeutic strategies against cancer cells with these phenotypes.
Materials and Methods
The logical formalism allows to handle the complexity of networks and the lack of quantitative data for most regulatory mechanisms (13). In brief, it considers an influence graph where the nodes embody regulatory components (e.g., genes, proteins, or phenomenological properties), and signed edges represent regulatory interactions activating or inhibiting the target node. Each regulatory component is associated with a discrete variable (Boolean or multivalued). This variable represents the component qualitative functional level (or state), that is, depending on the nature of the component, its level of expression, of activity, the effective complex formation or phenotype acquisition etc. Finally, a logical regulatory rule defines the values of this variable depending on the states of the components influencing it (i.e., its regulators). A model state is defined by a tuple of the variable values.
The model evolution is defined by the logical rules that, for any model state, indicate whether each component level should be updated or not. Here, we used an asynchronous updating scheme, meaning that all concurrent transitions are generated, each transition corresponding to the update of a single variable. Hence, any model state has as many successors as the number of components that are called to update in that state. Note that such an asynchronous updating scheme leads to nondeterministic dynamics (13).
The asymptotic behavior of the model is defined by its attractors. These are sets of mutually reachable states in which the model dynamics is trapped (no outgoing transitions). These attractors are defined by: a single stable state, with no successors, that is, all the model components are stable in that state, or by several states, corresponding to a maintained oscillatory behavior. The model presented here has only stable states and no cyclical attractors.
Computational tools and methods
The logical model of the EMT regulatory network (Fig. 1) was built using GINsim (version 3.0.0b, http://ginsim.org/), a software dedicated to the logical formalism (14). GINsim implements the determination of the stable states, export facilities, and other relevant features.
To validate the classification of the stable states into phenotypic categories based on the model readouts, we used the multiple correspondence analysis (MCA), which extends the principal component analysis to categorical (e.g., Boolean or discrete) variables and show how phenotypes cluster with respect to one another (15). This analysis was performed with R (version 3.5.0) using the two packages FactoMineR (16) and factoextra (https://rpkgs.datanovia.com/factoextra).
To confirm the absence of cyclical attractors, the model was reduced and exported to BoolSim (https://www.vital-it.ch/software/genYsis), a tool that efficiently constructs the dynamics of Boolean models (17).
In nondeterministic, asynchronous dynamics, a reachability probability can be calculated for the model attractors, assuming that concurrent transitions are equiprobable. Such a quantification was performed using the built-in GINsim functionality implementing Avatar, a modified Monte Carlo simulation (18), with a number of runs set to 105, ensuring the convergence of estimated probabilities.
The capacity of a system to adapt to environmental variation (plasticity) can be explored through model-checking analyses (19). We used the model checker NuSMV-ARCTL (20) with action restricted computation tree logic (ARCTL) temporal operators to study switches between phenotypes, while specifying input values (see Supplementary Materials and Methods; Supplementary Figs. S7 and S8). To ease this analysis, we used the reduction method available in GINsim, which allowed to propagate the values of the inputs FAT4L = 1, WNT = 0, and DELTA = 0 (21).
More quantitative view of the dynamics is provided by stochastic simulations as performed by MaBoSS (https://maboss.curie.fr/; ref. 22). The model was exported into MaBoSS format (this export feature is available in GINsim). MaBoSS computes stochastic trajectories and provides the time evolution of probabilities of the component values. We considered equal transition rates, a time step of 0.1, and a maximum time of 50. This allowed to plot the evolution of phenotype probabilities (Figs. 5B and 6B and C), which can be interpreted as the dynamical constitution of a cell population (each run corresponding to a cell).
Cell lines, culture conditions, and drug treatment
The MCF10A-ER-SRC cell line was kindly provided by K. Struhl (Harvard University, Boston, MA; ref. 23) and was cultured as described in ref. 24. To treat cells with 4OH-TAM (tamoxifen) or EtOH, 50% confluent cells were plated and allowed to adhere for at least 16 hours before being treated with 1 μmol/L 4OH-TAM (Sigma-Aldrich; H7904) or with identical volume of EtOH for the time period indicated. MDCK-pTR cSRCY527F-GFP cell line was kindly provided by C. Hogan (European Cancer Stem Cell Research Institute, Cardiff University, United Kingdom). Cells were treated with 2 μg/mL tetracycline and cultured as described in ref. 25. All experiments were performed with passage number between 8 and 25.
Collagen gel preparation
Collagen gels were prepared using Corning Collagen I High Concentrated, Rat Tail (Corning; 354249), according to the manufacturer's instruction. High-concentrated collagen (0.15 mL of 1 mg/mL or of 5 mg/mL) was used to coat the most central 14-mm diameter of 35 mm petri dishes (MatTek Corporation; P35G-1.5-14-C), which were precoated at the periphery with 1.2 % Poly-HEMA (Sigma-Aldrich; 3932).
Generation of stable cell line overexpressing PTPRK
MCF10A-ER-SRC cells were plated at 50% confluence and allowed to adhere for 16 hours. Cells were then infected with the PTPRK guide RNA or control lentiviral activation particles (Santa Cruz Biotechnology; sc-402287 LAC and sc-437282) in complete growth medium containing 8 μg/mL Polybrene (Millipore; TR-1003-6). Infected cells were then selected using media containing 400 μg/mL of Hygromycin B (Gibco/Thermo Fisher Scientific; 10687-010), 5 μg/mL of Blasticidin S HCl (Santa Cruz Biotechnology; sc-204655A), and 2 μg/mL of puromycin dihydrochloride (Calbiochem; 540411).
Phase-contrast microscope images
MCF10A-ER-SRC cells plated on collagen gels were imaged before or after treatments with EtOH or tamoxifen for 48 hours. Phase-contrast images were obtained on Leica DMi1 inverted microscope equipped with a Leica MC170 HD Camera (Leica Microsystems), using a 10 × 0.22 PH1 objective. Numbers of isolated cells per well were evaluated by applying the following settings: subtract background; Gaussian blur - sigma = 4. Single cells (size = 500–20,000 pixels and circularity = 0.80–1.00) and total cell (size = 500–inf pixels and circularity = 0.00–1.00) were then measured using the Analyse Particles Tool from Fiji.
A total of 20,000 cells, resuspended in 200 μL of complete cell culture media, were plated on top of a 0.7% agarose solidified layer. Images were acquired 24 hours later on a commercial Leica High Content Screening Microscope Leica DMI6000 equipped with a Hamamatsu Flash Orca 4.0 sCMOS Camera (Leica Microsystems), using a 10 × 0.30 NA objective. Numbers of isolated cells per well were evaluated on single microscopic fields using the Analyse Particles Tool from Fiji. Cell area between 500 and 2,000 pixels or between 2,000 and 10,000 pixels were labeled as isolated or aggregated cells, respectively.
Primary antibody used was rabbit anti-ECad monoclonal antibody 24E10 (1:50; Cell Signaling Technology, 3195S). For extended immunofluorescence analysis methods, see Supplementary Materials and Methods.
Primary antibodies used were rabbit anti-pSRC (pY419) (1:1,000; Invitrogen; 44–660G), rabbit anti-pFAK (pY397) (1:1,000; Cell Signaling Technology; 3283S), mouse anti-total FAK (1:500; BD Transduction, 610088), mouse anti-ECad (1:1,000, BD Transduction, 610182), rabbit monoclonal anti-RPTPK (1:1,000, clone H4; ref. 26), or rabbit anti-GAPDH (1:5,000, Sigma-Aldrich; G9545). For extended immunoblot analysis methods, see Supplementary Materials and Methods (26).
Atomic force microscopy
A PicoPlus 5500 Atomic Force Microscope (Keysight Technologies), coupled to an Inverted Fluorescence Microscope (Observer.Z1, Zeiss), was used for force mapping acquisition, in force spectroscopy mode. A MLCT-BIO-DC-D Cantilever (Bruker) was used. Before measuring collagen stiffness, the tip calibration was performed using the thermal K method obtaining a k = 0.048 N/m. Briefly, a laser beam was focused on the back of the cantilever, deflection of it was transduced in changing of the beam position in a photodiode. The calculation of the applied force was determined according to the Hooke's law. Once the laser aligned, three samples of 1 mg/mL or 5 mg/mL collagen I gels were measured in PBS 1×. Force maps were acquired using a grid of 6 × 6 points in an area of 5 × 5 μm2. The maximum force applied was 4.9 nN, with a loading rate cycle of 13 μm/second. The force curve maps were analyzed using the AtomicJ software 1.8 version, using the following parameters: blunt pyramidal model, Poisson ratio 0.5, and log-normal distribution. The collected approach force–displacement curves were fitted with the Hertz model to obtain the Apparent Young moduli. The topography images of each collagen concentration nanosurfaces were acquired using the Tapping mode in PBS 1×, with a scan speed of 1 Hz, in an area of 5.1 × 5.1 μm2, displayed at 2D and 3D images using the WSxM5.0 software (27).
Sample sizes and statistics
Quantifications of PTPRK and PTPRM mRNA levels in cells treated with EtOH or tamoxifen for 36 hours were from six biological replicates in triplicates. Quantifications of PTPRU or ECad mRNA were from four or three biological replicates, respectively, in triplicates. Quantifications of PTPRK mRNA over time and of SNAIL, SLUG, ZEB1, and ZEB2 mRNA levels were from three biological replicates in triplicates. Quantifications of PTPRK mRNA levels in cells expressing the mock or PTPRK guide RNA (PTPRK+), treated for 24 or 36 hours were performed in three or five biological replicates in triplicate, respectively. Quantifications of isolated cells by aggregation assays were from five biological replicates in triplicates. Quantifications of the number of isolated MCF10A-ER-SRC cells grown on stiff or soft collagen gels were from four biological replicates in triplicates. All analyses were performed using Prism 7.0 Software (GraphPad Inc.). For statistical comparison of two or more than two independent groups unpaired t tests or one-way ANOVA tests were used, respectively. An unpaired t test with Welch correction was used to compare the stiffness between soft and stiff collagen gels by an atomic force microscope. A two-tailed significance level of 5% was considered as statistically significant (P < 0.05).
Data and software availability
The computational model is available in the GINsim repository at http://ginsim.org/model/EMT_Selvaggio_etal. GINsim and other software tools used for the model definition and analysis are freely available (see section Materials and Methods). Scripts used to facilitate the model-checking analysis (Supplementary Materials and Methods) are available upon request.
Construction of a microenvironment-dependent EMT regulatory network revealing cell adhesion properties
We built a minimal regulatory network encompassing EMT-TF, epithelial (ECad and miR200, green in Fig. 1), and mesenchymal (SNAIL, SLUG, ZEB, TCF_LEF, and BCat, dark brown in Fig. 1) markers, as well as known EMT signaling pathways (RAS, NOTCH, WNT, TGFB, JAK/STAT, Hippo, and AKT), controlled by inputs from the tumor microenvironment (Supplementary Table S1). We focused on molecular interactions validated in experimental carcinoma models, prioritizing when possible human data. The regulatory rule of each internal component was defined on the basis of experimental evidence when available. If not specified otherwise, a component activation requires the presence of at least one of its activators combined with the absence of all its inhibitors (Supplementary Table S1).
Among the plethora of noncell autonomous signals involved in EMT, we included HGF, EGF, TGFB, IL6, ROS, WNT, and DELTA. In addition, to evaluate the impact of cell–cell and cell–ECM interactions, we incorporated the FAT4 ligand (FAT4L), the RPTP ligand of the R2B subfamily (RPTPL) and the ECM. While the ECM levels indicate the elastic property of the substrate, either stiff or soft (ECM = 1 or 0, respectively), all other input levels convey the presence (1) or absence (0) of these signals in the microenvironment. These 10 signals denoting environmental conditions constitute the input components of the model (gray nodes in Fig. 1).
Internal components (intracellular components) were associated with binary levels (present/active = 1, absent/inactive = 0), with the exception of FAK_SRC and ECad_AJ. The former was assigned a three-valued level based on the observations that increasing levels of SRC activity correlate with the acquisition of premalignant and malignant carcinoma features, respectively, in vivo and in vitro (24, 28). Because microarray profiling and RNA-sequencing indicate that oncogenic SRC downregulates RPTPK, the gene encoding for RPTP-kappa (29, 30), a predicted inhibitory interaction between FAK_SRC and RPTP was included, such that this interaction takes place only at high FAK_SRC activity (FAK_SRC = 2; blue interaction in Fig. 1).
ECad_AJ, the other multivalued component, conveys the state of ECad at the membrane: ECad_AJ = 0 indicates that the protein is not located at the membrane; ECad_AJ = 1 stands for the correct localization and binding to BCat_AJ; and finally, ECad_AJ = 2 represents the complete formation of the stable complex ECad, BCat, and p120.
Two multivalued readouts define the cell commitment in assembling AJs, an attribute of an epithelial state, and in remodeling FAs, typical of migrating properties displayed by mesenchymal cells (2). AJ values convey the capacity of ECad and BCat in assembling stable AJs, and FA values denote the status of FA recycling, and consequently the cellular migratory capacities (Supplementary Table S1).
Modeling cell adhesion properties accounts for phenotypical repertoire observed in EMT
The model attractors were identified: there were 1,452 stable states and no cyclic attractor. For conciseness, we focused on the 136 stable internal states obtained by considering only the internal components and ignoring the input components. Furthermore, we compressed these states by using wildcards (*, ?, and !), resulting in 31 patterns, each gathering multiple stable internal states (Table 1, see also Supplementary Materials and Methods). These patterns were then grouped together into eight phenotypes, each characterized by specific values of the readout components AJ and FA:
(i) States in which ECad localized with BCat at AJs, while displaying basal levels of FA recycling (AJ = 2 and FA = 0), were defined as epithelial phenotypes (E1).
(ii) States in which epithelial and mesenchymal markers were present were classified as hybrid phenotypes (H1, H2, and H3). Among these, H1 assembled AJs and showed the presence of both TCF_LEF and BCat and a weak ability to recycle FAs (AJ = 2 and FA = 1). The H2 phenotype failed to assemble AJs but exhibited the presence of ECad, miR200, TCF_LEF, and BCat with an intermediate ability to recycle FAs (AJ = 1 and FA = 2). Finally, H3 assembled AJs, while activating TCF_LEF and BCat and showing a high ability to recycle FAs (AJ = 2 and FA = 3).
(iii) States having only mesenchymal markers (and no epithelial markers) were classified as mesenchymal phenotypes (M1, M2, and M3). These phenotypes differed in their increasing capacity of recycling FAs (AJ = 0 and FA = 1, 2, or 3).
(iv) Four states referred to an undefined phenotype (UN) with the presence of mesenchymal markers only, and basal levels of FA recycling (AJ = 0 and FA = 0).
Note: Stable states are grouped into patterns (the table rows) classified into four phenotypes according to the model readouts AJ and FA. E, epithelial; H, hybrid; UN, undetermined; M, mesenchymal. Cell colors denote the levels of the components: 0 (red), 1 (green), 2 (dark green), 3 (olive). The wildcard * (white) indicates any admissible level; the coupled wildcard ? (light gray) indicates any same admissible level for the components marked with this symbol (the same applies for the coupled wildcard ! in dark gray). The penultimate column indicates the different phenotypes and their respective color codes. For each phenotype, the total number of stable states is given in the last column; this is calculated as the sum of the number of states represented in each pattern that depends on the wildcards: 2⁁(number of *) if there are no coupled wildcard; 2 × 2⁁(number of *) if there are only coupled wildcards ?; 2 × 2 × 2⁁(number of *) if there are both coupled wildcards ? and ! (see Supplementary Materials and Methods).
This phenotypical classification based on cell adhesion properties (model readouts) allowed to draw a parallel with the continuum of phenotypes along the EMT axis including hybrid phenotypes, such as H3, typical of a collective cell migration behavior (2).
Interestingly, the model pointed toward the epithelial phenotype as a reference phenotype in the absence of specific external signal: simulations starting from any state with all inputs set to 0 evolve toward an epithelial stable state.
To validate this classification, we performed a MCA, displaying the 136 stable internal states in a low dimension space (Fig. 2). This analysis showed that states adequately cluster together within each phenotype, with a clear separation between those clusters. We also tested the robustness of the model by setting the functional levels of some components to their maximum (i.e., X = 1 or 2) or to their minimum (i.e., X = 0), mimicking gain or loss of function conditions, respectively. The phenotypes reached by the perturbed model reproduced experimental observations reported in the literature (Supplementary Fig. S2). Interestingly, for some perturbations, the phenotypic landscapes were reshaped, leading to novel phenotypes and/or impairing the stability of others. We then performed a systematic analysis of single and double mutants. As hybrid phenotypes may provide pluripotent abilities to cancer cells (2), we searched for perturbations impeding these phenotypes (Supplementary Table S2). Three single mutants (TGFBR E1, BCad KO, and TCF_LEF KO) promoted a shift toward mesenchymal phenotypes or toward the two extreme phenotypes in the EMT continuum. In contrast, when excluding the effect of these single mutants, 43 double mutants triggered a shift toward mesenchymal only, and 20 toward mesenchymal and epithelial phenotypes (Supplementary Table S2).
Taken together, the qualitative assessment of cell adhesion properties defined in our computational model is capable of accounting for the phenotypic repertoire observed in the EMT continuum and permits to formulate predictions discussed hereafter.
SRC downregulates PTPRK
FAK and SRC integrate signaling from numerous microenvironmental inputs and play a central role in tumor-associated EMT (31). A predicted inhibitory interaction between FAK_SRC and RPTP has been included in the model (Fig. 1). Because of the latency period between SRC activation and the downregulation of PTPRK by microarray and RNA-sequencing (29, 30), we assumed that this interaction takes place exclusively at high levels of FAK_SRC activity (FAK_SRC = 2). Simulations of the overexpression of FAK_SRC (FAK_SRC = 2) indicated that when starting from an E1 phenotype with input configurations set to [ECM, GF, IF, CC, WNT, DELTA] =  and in the presence of an inhibitory interaction of FAK_SRC on RPTP, the model reached a full M3 mesenchymal phenotype in 100% of simulation runs (Fig. 3A). In contrast, when forcing RPTP maintenance, we observed a switch to the H3 phenotype in 50% of cases (Fig. 3B), suggesting that the repression of RPTP by FAK_SRC promotes a full mesenchymal phenotype.
To validate experimentally the relevance of RPTP repression by FAK_SRC on EMT, we took advantage of the MCF10A-ER-SRC cell line with conditional SRC activation. This normal cell line contains a fusion of v-SRC and the ligand-binding domain of the estrogen receptor, inducible with tamoxifen treatment. Upon induction, these cells acquire transformed features within 24–36 hours and invade in three-dimensional Matrigel cultures (23, 24). Consistent with an interdependency between SRC and FAK, tamoxifen-treated MCF10A-ER-SRC grown on plastic, displayed stepwise increases of phosphorylated ER-SRC (ER-pSRC) and endogenous SRC (pSRC; Fig. 3C), and of phosphorylated FAK (pFAK; Fig. 3D). Moreover, these cells acquire EMT features, including a significant downregulation of ECad mRNA levels, which was reduced by 48% 24 hours after tamoxifen treatment, and dropped by 72% at 36 hours (Fig. 3E). In addition, tamoxifen-treated MCF10A-ER-SRC cells suffered a morphologic shift from epithelial to mesenchymal starting 24 hours after treatment (Fig. 3F). This was in contrast to control MCF10A-ER-SRC cells, which maintained an epithelial morphology with ECad localized at cell–cell contacts during the 36 hours of EtOH treatment (Fig. 3G). Finally, the set of genes deregulated in tamoxifen-treated MCF10A-ER-SRC cells showed a significant 4.47 enrichment (13%, 33/246 genes; P < 0.0001; Hypergeometric test) for core EMT signature identified in mammary cells (32).
We then analyzed whether conditional SRC activation affected the expression of the four RPTP members of the R2B subfamily, encoded by the PTPRM, PTPRK, PTPRU, and PTPRT genes. Consistent with microarray and RNA-sequencing data (29, 30), the ratio of PTPRK mRNA levels between cells treated with tamoxifen and EtOH indicated that PTPRK was significantly downregulated by 33% 8 hours after SRC activation, and dropped by 70% at 36 hours (Fig. 3H; Supplementary Fig. S3A). Accordingly, the levels of the P subunit of RPTPK, generated from one precursor protein, displayed a progressive reduction during the 36 hours of transformation in tamoxifen-treated cells (Fig. 3I). In contrast, SRC activation had no significant effect on PTPRU or PTPRM mRNA levels (Supplementary Fig. S3B and S3C), while PTPRT mRNA could not be detected using two independent sets of primers. Although PTPRK expression is known to increase with cell confluency (33), the reduction in PTPRK levels by SRC is unlikely due to a density-dependent regulation of PTPRK expression, as it precedes the loss of cell–cell adhesion (Fig. 3D). Thus, we conclude that PTPRK downregulation results from SRC activation.
Forcing PTPRK expression in SRC-overactivated cells restores cell–cell adhesion
To test the model prediction by which maintaining PTPRK expression in SRC-overactivated cells favors the acquisition of an H3 phenotype (Fig. 3B), we forced PTPRK expression in MCF10A-ER-SRC cells using the clustered regularly interspaced short palindromic repeats–based activation system. At 24 hours after EtOH or tamoxifen treatment, MCF10A-ER-SRC cells stably expressing the PTPRK guide RNA (PTPRK+) displayed significantly higher PTPRK levels, compared with those expressing the mock guide RNA (mock; Fig. 4A). However, forcing PTPRK expression did not prevent the repressive effect of SRC on PTPRK expression 24 (Fig. 4A) and 36 hours (Fig. 4B) after tamoxifen treatment. Nevertheless, it was sufficient to restore PTPRK expression in tamoxifen-treated cells to levels similar to those of control EtOH-treated cells expressing the mock (Fig. 4A and B). Accordingly, the P subunit of RPTPK accumulated at higher levels in MCF10A-ER-SRC cells expressing PTPRK+ and treated with EtOH or tamoxifen for 24 or 36 hours (Fig. 4C). Then, using cell aggregation assays we tested whether restoring PTPRK expression would impact on the loss of cell–cell adhesion induced by SRC activation. MCF10A-ER-SRC cells expressing the mock and treated with EtOH for 36 hours were mainly found aggregated, with only 13.6% of isolated cells (Fig. 4D). As expected, treating these cells with tamoxifen prevented the formation of aggregates and increased the number of isolated cells to 54%. Strikingly, restoring PTPRK expression in tamoxifen-treated MCF10A-ER-SRC cells significantly restored cell aggregation with the number of isolated cells dropping to 36.6% (Fig. 4D). This effect did not result from a negative feedback of RPTPK on SRC or FAK activity, as MCF10A-ER-SRC cells overexpressing PTPRK still maintained high levels of ER-pSRC (Fig. 4E) and pFAK (Fig. 4F) levels, 24 or 36 hours after tamoxifen treatment. Taken together, these observations support the predicted requirement of the downregulation of PTPRK by SRC to promote the emergence of a full mesenchymal phenotype.
ECM stiffening synergizes with SRC to induce a full mesenchymal phenotype
Simulations indicate that stable states reached upon FAK_SRC overactivation would depend on the ECM status. In a soft ECM, and with an input configuration set to [ECM, GF, IF, CC, WNT, DELTA] = , FAK_SRC overactivation (FAK_SRC = 2) failed to generate an M3 phenotype when starting from an E1 state. Instead, the M2 and H2 phenotypes were reached in 22% and 78% of simulation runs, respectively (Fig. 5A). In addition, the model indicated that SNAIL (Fig. 5B), SLUG, and ZEB (Supplementary Fig. S4A and S4B) were expressed in 100% simulation runs with a stiff ECM (ECM = 1). In contrast, with a soft ECM (ECM = 0), SNAIL (Fig. 5B) was never expressed, while SLUG and ZEB were expressed in 20% of cases only (Supplementary Fig. S4A and S4B). Altogether, this would indicate a synergistic effect between ECM stiffening and FAK_SRC to generate a full mesenchymal phenotype, and to upregulate EMT-TF.
To test this prediction experimentally, we analyzed the effect of type-I collagen-coated rigid plastic of distinct concentrations on the behavior of tamoxifen-treated MCF10A-ER-SRC cells. Collagen gels (1 mg/mL) with average Young moduli of 1.641 kPa ± 0.691 were used to mirror soft ECM surrounding normal mammary epithelial cells (34). While 5 mg/mL collagen gels of average Young moduli of 3.084 kPa ± 1.458 (Supplementary Fig. S5A–S5C) were used to resemble stiffer matrix reported for stroma adjacent to transformed cells (34). Consistent with the model prediction, tamoxifen-treated MCF10A-ER-SRC cells expressed significantly higher levels of SNAIL when grown on stiff gels, compared with those plated on soft gels (Fig. 5C). However, the presence of a stiff gel was not sufficient to significantly induce SLUG, ZEB1, or ZEB2 expression in tamoxifen-treated MCF10A-ER-SRC cells (Supplementary Fig. S6A–S6C). The ECM stiffness-dependent effect on SNAIL expression in this cell line was also associated with distinct cell morphologies. While on stiff gel, tamoxifen-treated cells appeared rounder and isolated, when grown on soft gels, these cells maintained cell–cell contacts at “tip-like” junctions (Fig. 5D). These differential effects were not because of increased SRC activity following matrix stiffening, as the levels of ER-pSRC did not significantly vary in tamoxifen-treated MCF10A-ER-SRC between those grown on stiff and soft gels (Fig. 5E). Moreover, on both, soft and stiff gels, SRC activation could downregulate PTPRK (Fig. 5F). To confirm the synergistic effect between ECM stiffness and SRC, we analyzed the effect of soft and stiff collagen gels on the Madin–Darby Canine Kidney (MDCK) epithelial cell line expressing cSRCY527F in a tetracycline-inducible manner (MDCK-pTR cSRCY527F-GFP). In contrast to tetracycline-treated MDCK-pTR cSRCY527F-GFP cells grown on soft gels, those plated on stiff gels poorly accumulated ECad at the cell membrane (Fig. 5G). Moreover, in this cell line, a stiff gel synergized with cSRCY527F to significantly upregulate SNAIL (Fig. 5H) and SLUG (Supplementary Fig. S6D). ZEB1 and ZEB2 could not be tested in this cell line using two independent sets of primers.
Taken together, our experimental assessments of the effect of SRC and ECM stiffness are consistent with the model predictions by which ECM stiffening synergizes with SRC to potentiate SNAIL and SLUG expression and to promote a mesenchymal phenotype.
Phenotypes are plastic upon changes in environmental conditions
To further analyze the impact of microenvironment signals on EMT dynamics, we performed systematic reachability analyses of the model stable states upon switching input combinations. For concision, EGF and HGF and ROS on the one hand, IL6 and TGFβ on the other hand were gathered to, respectively, denote the presence or absence of growth factors (GF) and of inflammation (IF). The RPTPL value denoted the state of cell–cell contact (CC) and the other inputs were assigned fixed values with DELTA and WNT, both considered absent, and FAT4 considered present (FAT4 = 1). We were thus left with four input configurations (ECM, GF, IF, and CC). The model was then reduced by propagating the values of the fixed inputs, and stable states reached by the model were classified into distinct phenotypes, depending on the values of the output nodes AJ and FA, as shown previously (Table 1). Note that setting WNT to value 0 discarded the UN phenotype, which requires the presence of the WNT signal.
Considering each input configuration and for each phenotype, we checked the following properties (Supplementary Materials and Methods):
(i) Does it exist as a state in the set of states defining the phenotype such that there is a trajectory leading directly to a target phenotype (i.e., without visiting any other phenotype)? This resulted into a “weak” reprogramming graph (Supplementary Fig. S7);
(ii) For every state of the phenotype, is there a trajectory leading directly to a target phenotype (i.e., without visiting any other phenotype)? This resulted into a “strong” reprogramming graph (Supplementary Fig. S8).
When property (i) does not hold for an input configuration, it means that all the states of the phenotype are stable for that configuration. This “strong” stability property is reported as a self-loop in the “strong” reprogramming graph (Supplementary Fig. S8). From this graph, we discarded the phenotypes M1 and H1, as no input configurations ensure the stability of all the states of these phenotypes. This permitted to obtain a simplified version of the strong reprogramming graph, which displays striking behaviors (Fig. 6A). First, while the H2, H3, and M2 phenotypes can evolve into all other phenotypes, this is not the case for the E1 and M3 phenotypes. Moreover, none of the inputs considered (ECM, GF, IF, and CC) permit to attain H3 neither from E1 nor from M3, or to reach H2 from M3. However, some specific states in the other phenotypes can lead to H2 and H3 (Supplementary Fig. S7). In contrast, while the E1 and M3 phenotypes show strong stability, a large variety of environmental inputs allow to leave the H2, H3, or M2 phenotypes. Another relevant observation is that switching off ECM is necessary and sufficient to convert the M3 phenotype into M2 and/or E1 phenotypes (Fig. 6A). Again, the “weak” reprogramming graph shows that other phenotypes are possibly reached from restricted sets of M3 states (Supplementary Fig. S7).
Relying on the “strong” reachability graph, we then explored possible scenarios by which microenvironmental inputs could trigger EMT. We assumed that (i) alterations in microenvironmental inputs are more likely to arise sequentially (only one input altered at a time) and; (ii) ECM stiffening is a late event, as it has been mainly associated with later stages of cancer progression (35). Starting from an E1 phenotype with input configurations set to [ECM, GF, IF, CC] = [0001 or 0000], four possible scenarios were achievable.
When the input IF was first switched on in the absence of CC (configuration ), this led to transient occurrences of the H1, H2, and M1 phenotypes, that all ultimately reached the M2 phenotype. Then, switching the ECM to 1 (configuration ) led to the disappearance of the M2 phenotype, which was replaced by M3 (Fig. 6B). In this scenario, the presence of CC (configuration ) had only subtle effects on the trajectory. The only divergence was the transient emergence of the M3 rather than the M1 phenotype (Fig. 6B). Moreover, the kinetics and probabilities of acquiring these phenotypes were not altered by switching GF to 1 after stabilization of the M2 or M3 phenotypes by IF or ECM, respectively. In contrast, when GF was first switched on in the absence of CC (configuration ), this led to a transient appearance of an H1 phenotype, and ultimately a mixed stable population with about 5% and 95% of cells with M2 and H2 phenotypes, respectively (Fig. 6C). Then, switching IF to 1 (configuration ) and then ECM to 1 (configuration ) converts H2 into a M2 and then a M3 phenotype in 100% of the cells (Fig. 6C). Strikingly, in this scenario, the CC status drastically influences the response of epithelial cells, as in the presence of CC all cells maintained an E1 phenotype upon GF stimulation (configuration ; Fig. 6C). Several predictions arise from this analysis: (i) IF hampers the maintenance of the hybrid H2 and H3 phenotypes while favoring the mesenchymal M2 and M3 phenotypes; (ii) ECM stiffening promotes the emergence of the H3 and M3 phenotypes; and (iii) CC (illustrated by RPTPL) prevents EMT in the presence of GF.
Our in silico model uses cell adhesions properties as readouts for the acquisition of EMT traits. It robustly recapitulates cell behaviors observed experimentally. First, the model accounts for the phenotypic repertoire identified using experimental and mathematical modeling assessments of EMT (7, 9, 36, 37), including pure epithelial (E1) and mesenchymal (M) phenotypes, as well as hybrid phenotypes (H). Second, the ease by which these hybrid phenotypes evolve upon microenvironmental stimulations could reflect the observed pluripotent abilities of cancer cells in these hybrid states (2). Third, we experimentally show that the downregulation of PTPRK by SRC is required to promote the emergence of a full mesenchymal phenotype. Fourth, our experimental assessment of the effect of ECM stiffening corroborates the model prediction by which a stiff ECM and FAK_SRC synergize to potentiate SNAIL and SLUG expression, therein triggering a full mesenchymal phenotype. Our observations that SLUG, ZEB1, and ZEB2 did not respond to increase gel stiffness in tamoxifen-treated MCF10A-ER-SRC could reflect cell type–specific effects of these EMT-TFs, as reported in other context (38, 39). In any case, as wisely underscored by Jolly and colleagues, “no one-size-fits-all,” that is, no model, either biological or mathematical, can account for all biological contexts (37). Fifth, consistent with our simulations indicating that ECM stiffening boosts cell movement, irrespectively of the mode of migration (collective through H3 vs. mesenchymal through M3), a stiffer microenvironment promotes metastatic transition (34, 40), induces faster migration (41), and does not influence the mode of cell migration (42). Sixth, as observed in malignant tissues and cultured cancer cell lines (43), our model suggests the coexistence of cell populations with distinct phenotypes, revealing tumor heterogeneity.
Our model also allows to formulate diverse predictions. Hybrid phenotypes might not constitute mere transient steps along an EMT continuum, but could represent independent trajectories, as the acquisition and maintenance of the H2 and H3 phenotypes require more stringent inputs than those required for gaining mesenchymal traits. In addition, simulations of single or double mutants and of microenvironmental signals on EMT phenotypes predict key players that would facilitate or hinder the maintenance of these hybrid phenotypes. As these phenotypes have been associated with worse progression-free survival in patients with cancer (44), our work point to actionable molecular targets that could serve to design novel therapeutic interventions. In particular, inhibitors against molecular targets that would drag cancer cells into mesenchymal phenotypes only could be beneficial for patients with cancer, as they could also halt metastasis seeding. Among involved microenvironmental signals, inflammation triggers the transient appearance of these phenotypes, but does not allow their maintenance. Instead, inflammation favors the acquisition of mesenchymal phenotypes. Accordingly, the autocrine production of inflammatory signals, such as TGFβ, stabilizes a mesenchymal phenotype in vitro (45). However, our model suggests that transient versus sustained inflammatory signals would have different consequences on tumor cell behavior and likely aggressiveness.
In contrast to inflammatory signals, the model predicts that growth factors support the acquisition and maintenance of the hybrid H2 phenotype, if and only if cell–cell contacts through RPTP are absent. Thus, RPTPs could fulfill tumor suppressor functions in the presence of growth factors. Consistent with this hypothesis, the downregulation of PTPRK in solid cancers correlates with a poor disease free-survival time (5). Moreover, knocking down PTPRK increases invasiveness in breast cancer cells, whereas expressing PTPRK exogenously in melanoma cells reduces cell migration (46, 47). Furthermore, we provide evidence that PTPRK is downregulated by SRC. The model also predicts that RPTP is absolutely required to support an H3 phenotype. In agreement with this prediction, R2B RPTPs have been proposed to promote the collective movement of tumor cells along nerve bundles (48). Moreover, they have also been classified as oncogenes (5). Furthermore, we show that forcing PTPRK expression in tamoxifen-treated MCF10A-ER-SRC cells restores cell–cell adhesion. Although PTPR downregulation is not an obligate requirement for a mesenchymal phenotype, our computational and experimental approaches suggest that PTPRK downregulation prevents collective cell migration to allow full mesenchymal transition. Thus, RPTPK could have opposite effects during EMT: early, it could prevent the induction of an EMT program, in particular in the presence of growth factors, whereas later, PTPRK could determine the mode of cell migration (collective vs. single cell migration).
Our in silico analysis suggests that, as observed experimentally (36), a soft ECM is not required to maintain an epithelial phenotype. It further predicts that ECM softening is sufficient to revert both mesenchymal M3 and hybrid H3 phenotypes into epithelial ones. Thus, alterations in ECM rigidity could be a critical determinant for seeding metastasis. Accordingly, diverse observations support a role of ECM remodeling during premetastatic niche formation (49).
Further model analyses would permit to identify noncell autonomous role of WNT and DELTA, as well as individual contributions of the growth factors EGF and HGF and of the inflammatory inputs IL6, TGFβ, and ROS. Moreover, as EMT alters many cell behaviors in addition to adhesion properties (2), the model could be extended by including relevant pathways and associated readouts, such as survival, proliferation, and other cell fates. As future work, we also aim at embedding the cellular model into a multicellular context, which will allow to properly simulate collective dynamics, accounting for cell–cell interactions and their synergistic shaping of the microenvironment.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Conception and design: F. Janody, C. Chaouiya
Development of methodology: S. Canato, M.M. Brás, F. Janody, C. Chaouiya
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Pawar, P.S. Guerreiro, M.M. Brás
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): G. Selvaggio, S. Canato, A. Pawar, P.T. Monteiro, P.S. Guerreiro, F. Janody, C. Chaouiya
Writing, review, and/or revision of the manuscript: G. Selvaggio, S. Canato, A. Pawar, P.T. Monteiro, P.S. Guerreiro, M.M. Brás, F. Janody, C. Chaouiya
Study supervision: F. Janody, C. Chaouiya
Other (funding acquisition): F. Janody, C. Chaouiya
We specially thank K. Struhl and C. Hogan for providing the MCF10A-ER-SRC and MDCK-Src cell lines, respectively, P. Hermanowicz for his assistance on the use of the AtomicJ program and for AFM data analysis, H.J. Shape for sharing their anti-RPTPK antibody, and R. Pais for defining a first version of the computational model. We acknowledge the Imaging and Cytometry and Genomics facilities at IGC, and the Cell Culture and Genotyping, the Genomics, the Bioimaging and the Advanced Light Microscopy platform, both members of the national infrastructure PPBI-Portuguese Platform of BioImaging (POCI-01-0145-FEDER-022122), as well as the Biointerfaces and Nanotechnology platforms at i3S. Funds from Fundação para a Ciência e Tecnologia (FCT; PTDC/BEX-BCB/0772/2014 to C. Chaouiya) supported this work and A. Pawar and S. Canato. Funds from FCT, co-financed by Fundo Europeu de Desenvolvimento Regional (FEDER) through Programa Operacional Competitividade e Internacionalização (POCI-01-0145-FEDER-016390 to C. Chaouiya and F. Janody), supported this work and G. Selvaggio. G. Selvaggio and C. Chaouiya were supported by the Gulbenkian Foundation. F. Janody was the recipient of IF/01031/2012.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.