Signaling network state predicts twist-mediated effects on breast cell migration across diverse growth factor contexts.

Epithelial-mesenchymal transition (EMT), whether in developmental morphogenesis or malignant transformation, prominently involves modified cell motility behavior. Although major advances have transpired in understanding the molecular pathways regulating the process of EMT induction per se by certain environmental stimuli, an important outstanding question is how the activities of signaling pathways governing motility yield the diverse movement behaviors characteristic of pre-induction versus postinduction states across a broad landscape of growth factor contexts. For the particular case of EMT induction in human mammary cells by ectopic expression of the transcription factor Twist, we found the migration responses to a panel of growth factors (EGF, HRG, IGF, HGF) dramatically disparate between confluent pre-Twist epithelial cells and sparsely distributed post-Twist mesenchymal cells-but that a computational model quantitatively integrating multiple key signaling node activities could nonetheless account for this full range of behavior. Moreover, motility in both conditions was successfully predicted a priori for an additional growth factor (PDGF) treatment. Although this signaling network state model could comprehend motility behavior globally, modulation of the network interactions underlying the altered pathway activities was identified by ascertaining differences in quantitative topological influences among the nodes between the two conditions.

In the phenomenon of epithelial-mesenchymal transition (EMT), polarized epithelial cells loosen their cell-cell junctions and acquire the ability to migrate through extracellular matrices as single cells in a mesenchymal manner (1, 2). Although great progress has been made on identifying and understanding components and mechanisms involved in the process of EMT induction (e.g. (3, 4)), the “before” versus “after” consequences of this transition for signaling pathway control of cell migration has not yet been investigated from a multipathway, network-wide perspective. Cell migration results from a set of carefully orchestrated biophysical processes regulated by numerous key signaling pathways whose activities can be influenced downstream of a range of growth factor receptors. It is appreciated that these growth factor receptor-elicited signaling activities may be modulated in “before” versus “after” manner by EMT induction (5), whether by TGFβ or other developmental cues or inflammation-related stimuli (6, 7). However, a current challenge is to characterize this likely complex modulation from a multipathway network perspective and to establish an approach for predictive understanding of how the multiple pathway activities integrate to yield different migration behavior in postinduction compared with pre-induction conditions. This challenge is especially important for, among other motivations, gaining insights concerning how prospective targeted drug effects are influenced by whether tumor cells are in epithelial or mesenchymal state (8).

As one currently clinically urgent application example, the epidermal growth factor receptor (EGFR) is commonly overexpressed or mutated in epithelium-derived tumors, and its activation is linked to progression and poor prognosis (9). Therefore, EGFR has been the target of many small molecule inhibitors and monoclonal antibody antagonists, which have met with limited clinical success (10–12). Recent studies exploiting EMT markers and gene expression signatures suggest that cells with low levels of epithelial markers, such as E-cadherin, and high levels of mesenchymal protein expression, such as N-cadherin and vimentin, display resistance against these inhibitors (13, 14). Therefore, the decreased sensitivity of mesenchymal-like tumors to EGFR antagonists argues for an ability to bypass EGFR dependence to activate the downstream signaling pathways necessary for cell migration and survival (15). Cell activation through other receptors including the insulin-like growth factor-1 receptor (IGF-1R), fibroblast growth factor receptor (FGFR), and platelet-derived growth factor receptor (PDGFR), has been suggested to play a role in resistance to EGFR antagonists (14, 16). Thus, improved understanding of how EMT-mediated changes in multiple growth factor signaling networks contribute to cell invasion may necessarily shift investigational focus toward the design of novel therapeutics targeting tangential tyrosine kinase pathways or intracellular signaling nexi for use in treating EGFR inhibition-resistant carcinomas.

As a first multipathway network level study of how signaling pathway activities governing cell migration downstream of receptor tyrosine kinase stimulation differ between “before EMT” and “after EMT” conditions, we use here an established human mammary epithelial cell line (hMLE) immortalized and transformed via introduction of a minimal set of oncogenes (17) and focus on EMT induction by Twist1 (18), via its ectopic expression in hMLEs as previously characterized (19). Twist expression has been demonstrated in multiple studies in vitro, in mouse models, and in human patients, to be associated with breast tumor invasiveness, metastasis, and poor disease prognosis (e.g. (19–22)), and thus represents a pathophysiologically and clinically important system for analysis. It also may be as simple an induction process as can be examined, because other EMT inducers such as TGFβ and TNFα act via multiple transcription factors including Twist along with others (7), so our initial study here may indicate basic signaling network modulation insights that can be expanded upon in future analogous investigations of the more pleiotropic EMT inducers.

In this basic study, we quantitatively characterize the migration characteristics of hMLEs before and after Twist-mediated induction in both monolayer (indicative of epithelial mode) and single cell (indicative of mesenchymal mode) migration assays under stimulation by a panel of growth factors present in carcinoma environments including EGF, HRG, IGF, and HGF (16, 23–25). Across this broad landscape of extracellular treatment conditions, we measured phosphorylation states of 14 signaling pathway nodes to ascertain how Twist-mediated changes in numerous of these signals may be associated with consequent changes in the cell motility behaviors. Computational modeling with a partial least-squares regression (PLSR) framework demonstrated that quantitative combinations of multiple signals can account for the various motility behaviors across all growth factor treatments in both epithelial and mesenchymal migration modes—and, in fact, can successfully predict a priori the motility behavior for epithelial and mesenchymal modes in a new growth factor context, PDGF stimulation. We then constructed a complementary computational model, using a correlative topology framework, to identify influences among the signaling nodes that were modulated by the Twist-mediated EMT induction.

EXPERIMENTAL PROCEDURES

Antibody Reagents, Growth Factors, and Inhibitors

Bio-Plex bead-based phospho-protein ELISA kits were purchased from Bio-Rad Laboratories (Hercules, CA) and are against phospho(p)-EGF receptor (pan-p-Tyr), p-IGF-1 receptor (Tyr1131), p-Src (Tyr418), p-Erk1/2 (Thr202/Tyr204, Thr185/Tyr187), p-HSP27 (Ser78), p-JNK (Thr183, Tyr185), p-Akt (Ser473), p-IRS-1 (Ser636/Ser639), and p-GSK3α/β (Ser21/9). Antibodies against p-PKCδ (Thr505), p-PLCγ (Tyr771), p-Met (Tyr1234/1235), p-β-catenin (Ser33/37, Tyr41), EGFR, IGF-1R, HER2, and GAPDH were purchased from Cell Signaling Technologies (Danvers, MA). The antibody against p-FAK (Tyr397) was purchased from Invitrogen (Carlsbad, CA). Antibodies against E-cadherin, N-cadherin, and vimentin, as well as human recombinant IGF-1, HGF, and PDGF-BB were purchased from BD Biosciences (San Jose, CA). Human recombinant EGF and HRG-β1 were purchased from Peprotech (Rocky Hill, NJ) and Sigma-Aldrich (St. Louis, MO), respectively.

Cell Culture

Immortalized human mammary epithelial cells (hMLEs) expressing either the empty pBabe puro vector (pBp) or pBP-Twist1 (referred to as epithelial or mesenchymal hMLEs) were obtained from the Robert Weinberg laboratory at the Whitehead Institute for Biomedical Research (Cambridge, MA) and cultured as described previously (26) with complete medium consisting of 50% mammary epithelial basal medium (MEBM) (Lonza; Walkersville, MD), 25% Ham's F-12, 25% Dulbecco's modified Eagle's medium (Invitrogen; Carlsbad, CA) and an MEGM Bulletkit (Lonza) containing bovine pituitary extract, EGF, insulin, hydrocortisone, and gentamycin. Serum-free medium contained all components of the complete medium, except bovine pituitary extract, EGF, and insulin. MDA-MB-231, MDA-MB-453, BT549, and T47D human breast carcinoma cells were obtained from ATCC (Manassas, VA) and cultured in Dulbecco's modified Eagle's medium supplemented with 10% fetal bovine serum (Invitrogen).

Quantitative Cell Signaling Analysis

Epithelial and mesenchymal hMLE cells were seeded in six-well or 10 cm tissue culture plates at a density of 50,000 cells/cm2 or 6000 cells/cm2, respectively, in complete medium. Cells were grown overnight and then starved in serum-free medium for 24 h. Cells were stimulated with either 100 ng/ml EGF (Peprotech; Rocky Hill, NJ), 80 ng/ml HRG (Sigma), 100 ng/ml IGF-1, 50 ng/ml HGF, or 50 ng/ml PDGF-BB (BD Biosciences; San Jose, CA). After times indicated, cells were lysed in Bioplex Lysis Buffer (Bio-Rad Laboratories; Hercules, CA), and lysate protein concentrations were determined with a microplate BCA assay (Pierce Biotechnology; Rockford, IL). Lysates were aliquotted at equal protein concentrations and kept at −80 °C.

We determined the optimal loading mass of lysates for the Bioplex bead-based ELISA and Western blot assays by validating for their suitability for quantitation by ensuring the linearity of the assay signal versus protein concentration: 20 μg per well for p-Src, p-p38, p-JNK and all Western blots; 5 μg for p-HSP27, p-IRS-1, p-GSK3α/β, p-Erk1/2, p-IGF1-R; 1 μg for p-Akt and p-EGFR. A positive control lysate was loaded for signal intensity comparison between blots or assay plates. Manufacturer's instructions were followed for the Bioplex assay. Protein lysates were resolved on a 4–12% NuPage Novex Bis-Tris gels (Invitrogen) with Laemmli Running Buffer and transferred onto nitrocellulose membranes. Blots were probed with various primary antibodies listed above. Western blot bands were visualized via IRDye secondary antibodies (Rockland Immunochemicals Inc, Gilbertsville, PA) and the Li-Cor Odyssey Imaging system (Lincoln, Nebraska) and subsequently quantified by densitometry using the Odyssey software (Li-Cor). All raw signaling data was eventually normalized to the respective 0 min time point phosphorylation in epithelial cells.

Cell Migration Assays

For quantitation of migratory behavior, cells were incubated with 4 μm 5-chloromethyl-fluorescein diacetate (CMFDA) (Invitrogen) in serum free media for 20 min and washed twice with PBS before seeding overnight, at seeding densities described above. For the sparse migration assay, only labeled cells were seeded. For the monolayer migration assay, labeled cells and unlabeled cells were seeded at a 1:20 ratio. Epithelial cells were seeded at twice the density as mesenchymal cells to account for the cell size differential. After incubating overnight, hMLE cells were serum-starved for 24 h before growth factor stimulation. Other human breast cancer cells were incubated for 24 h in full serum media without serum-starvation. After 1 h of growth factor stimulation, cells were placed on an environment-controlled Nikon TE2000 microscope (Nikon Instruments; Melville, NY). Cells were imaged with 30-min time intervals on differential interference contrast microscopy (DIC) and 488 nm excitation. Fluorescence movies were analyzed with Bitplane Imaris software (Zurich, Switzerland) using the built-in “Spots” function. 12-hour tracks (6 h after stimulation) were generated using the “Brownian Motion” algorithm and manually validated to obtain cell speeds.

Three-Dimensional Collagen Migration Assay

As described previously (26), cells were mixed with pH-neutralized collagen I in serum-free medium (Inamed Biomaterials; Fremont, CA) at 500,000 cells/ml and 2.0 mg/ml final collagen concentration. The matrix-cell solution was placed on glass-bottom cell culture dishes (MatTek; Ashland, MA) and polymerized for 1 h at 37 °C. The culture was immediately serum-starved for 24 h and stimulated with 100 ng/ml EGF before imaging cells via brightfield microscopy at 10 × magnification using the environment-controlled Nikon TE2000 microscope.

Partial Least-Squares Regression Modeling

The relationship between the experimentally determined cell signaling and cell migration response data sets was assumed to be linear and was modeled using partial least squares regression (PLSR). For this purpose, two metrics were extracted from each signaling time course: (1) the initial phosphorylation (before growth factor stimulation, referred to as T = 0) and (2) integral of activation over the entire time course (int). All average signaling measurements was assembled into the signaling data matrix (X) consisting of the signaling metrics as columns (14 signals × 2 metrics) and of the growth factor treatments across the cell lines as rows (2 cell lines × 5 growth factor conditions). Epithelial monolayer migration speeds and Mesenchymal sparse migration speeds were assembled into a migration response vector (Y) consisting of 12 rows. X and Y were mean-centered and unit-variance scaled within each column.

The PLSR model was implemented using the noniterative partial least squares (NIPALS) algorithm in SIMCA-P 11.0 (Umetrics; Kinnelon, NJ). A detailed description can be found elsewhere (27). Briefly, PLSR attempts to solve the linear regression problem Y = f(X) = X●b, where b is a vector of regression coefficients containing information regarding each signaling metric's contribution to the cell migration response (loadings). This regression problem does not have a unique solution because X is rank deficient (i.e. X consists of more columns than rows). PLSR arrives at the solution via a linear regression of the data sets in a reduced-dimensionality principal component space with regression coefficients associated with principal components.

Models were generated using two to three principal components under the standard optimization criteria. Model calibration was performed using leave-one-out cross-validation and model uncertainties were calculated by jack-knifing (28). The quality of the models was evaluated by the goodness-of-fit parameters, R2X and R2Y, which are calculated as the fraction of sum of squares of all the X and Y variables, respectively, explained by the number of principal components. The predictive capability of the model was quantified by Q2, which is the fraction of the total variation in the response variable according to cross-validation. Model loadings and Y-scores were calculated using mean-centered coefficients w x c and ui, respectively, from the i-th principal component. The projected coordinates of the loadings and scores are complementary and indicate corresponding co-variations in the principal component space. To predict cell speeds a priori under growth factor treatments not included in the training data set (such as PDGF), we extracted the two signaling metrics for this growth factor treatment and used model regression coefficients to calculate the cell speeds. The two-component model was found adequate for comprehending the full range of experimental studies, so we focused our analysis on it.

Correlation Network Modeling

Pairwise Pearson correlation was used to quantify the relatedness between signaling nodes in the epithelial and mesenchymal cell states. First, the geometric mean of the phosphorylation fold-change relative to time zero from two to three biological replicates was calculated for each nonreceptor phosphosite time-course. Using only the four nonzero time points across five growth factor treatments, this gave 20 data points per phosphosite per cell state. Given 11 nonreceptor phosphosites, the Pearson correlation was then calculated between each pair of phosphosites using this 11 × 20 data matrix. The p values for nonzero correlation were calculated using a Student's t distribution for a transformation of the correlation. This provided 11 × 10/2 = 55 unique pairwise correlation coefficients and p values, neglecting self-correlations. Three separate methods were applied to account for multiple hypothesis testing: Bonferroni (29), Benjamini (30), and Storey (31). Bonferroni is the most conservative, and Storey the least conservative, of these alternative methods, with respect to assigning statistical significance to correlations.

RESULTS

Diverse Cell Motility Behavior and Growth Factor Treatment Responses in Epithelial versus Mesenchymal Mode

hMLEs ectopically expressing a vector control or Twist1, a transcription factor previously shown to induce EMT, were used as a model of EMT-induced phenotypic switch (called “epithelial” or “pre-Twist”, versus “mesenchymal” or “post-Twist” cells hereafter). Cells were cultured in serum-free medium upon seeding to assess growth factor-stimulated cell migration. The cells in epithelial and mesenchymal modes maintained their respectively appropriate EMT markers in this medium (supplemental Figs. S1).

Although invasive carcinomas and cells of mesenchymal developmental origins may invade as single cells, epithelial cells can also migrate but do so within established monolayers. To consider both types of migration, we seeded cells labeled with whole-cell tracking dye either sparsely to achieve single-cell migration or in a confluent monolayer with unlabeled cells for migration with cell-cell contact (Fig. 1A and 1B). Upon serum-starvation, cells were treated with saturating levels of EGF. As anticipated, sparse post-Twist cells migrated significantly, whereas pre-Twist cells that were maintained as single cells throughout the experiment exhibited little movement (Fig. 1A). Pre-Twist cells with intact cell-cell contacts (Fig. 1A) or in a confluent monolayer (Fig. 1B) displayed significant locomotion, consistent with previous reports of mammary epithelial cells (32, 33). In contrast, post-Twist cells exhibited a contact-mediated reduction in motility. Moreover, consistent with clinical observations (15), post-Twist cells displayed resistance to inhibition of invasion via inhibition of EGF signaling (Fig. 1C). Similar differences in motility behavior were observed with respect to invasion into a three-dimensional collagen I matrix (Fig. 2A); post-Twist cells invaded to a significant extent whereas pre-Twist cells did not.

Fig. 1.

EMT and growth factor-dependent human mammary epithelial cell migration is contingent on its context, which is recapitulated by other human breast cancer cell lines. A, B, DIC and epifluorescence overlay and cell tracks of epithelial (left) and mesenchymal (right) cells in the sparse (A) and monolayer (B) migration assay. Cells were labeled with a whole-cell dye CMFDA and either seeded sparsely (A) or mixed with unlabeled cells and seeded in confluence (B) before a 24-hour serum-starvation and treatment with saturating levels of EGF. After 1 h of stimulation, migration tracks over 18 h (red) were generated via semi-automated tracking of centroids (gray circles) of labeled cells. C, Cell speeds of epithelial cells in monolayer migration assay (black) and mesenchymal cells in sparse migration assay (red) in presence of varying levels of an EGF receptor kinase inhibitor AG1478. AG1478 was added simultaneously to EGF. Cell speeds are normalized to their respective no inhibitor control cell speeds. p < 0.0001 via two-way ANOVA between cell lines. D, E, Cell speeds of epithelial (black) or mesenchymal (red) cells under stimulation of various growth factors quantified from the sparse (D) or monolayer (E) migration tracks. Cell speeds were calculated from cell tracks after 7 to 19 h of stimulation. All data is shown as mean ± S.E. Box-and-whisker plots of individual cell speeds are shown in supplemental Fig. S5. n = 109–175 (C) and n = 269–390 (D, E) cells for monolayer migration and n = 31–66 (C) and n = 16–117 (D, E) and cells for sparse migration obtained from three independent biological replicates. * p < 0.05, ** p < 0.01, *** p < 0.0001 compared with serum-free condition (D, E) or sparse condition (C) (see Materials and Methods for details on statistical analyses).

Fig. 2.

Epithelial (top) and mesenchymal (bottom) cells were seeded in a neutralized 2.0 mg/ml collagen I solution. Upon gelation of collagen I, cells were serum-starved for 24 h and stimulated with saturating levels of EGF. Cells were imaged via phase-contrast microscopy over 6 h. Arrows indicate migratory mesenchymal cells in the three-dimensional collagen I matrix. Dashed lines are provided as a reference. Details of methodology can be found in (23). B, Cell speeds of human breast cancer cell lines migrating in complete medium in a sparse or monolayer migration assay. Cell lines are color-coded according to their widely accepted EMT state; red for mesenchymal and black for epithelial. All data are shown as mean ± S.E.

We also considered whether this differential behavior might be generalized to other breast tumor cell lines and similarly examined motility behavior of a panel of breast carcinoma cell lines in both confluence and sparse conditions. We found that lines representing the luminal subtype (T47D, MDA-MB-453) showed an epithelial pattern of migration, whereas lines representing the basal subtype (BT549, MDA-MB-231) diverged in their pattern of migration (Fig. 2B). However, the respective levels of Twist expression can explain the latter divergence: the MDA-MB-231 cells, which exhibited epithelial-like migration pattern similar to the T47D and MDA-MB-453 cells, likewise express Twist at only low levels whereas the BT549 cells which exhibited mesenchymal-like migration pattern express Twist at high level (34, 35). Taken all together, these findings indicate that motility behavior of the mammary epithelial cells is substantively altered by Twist expression and that insights gained in our model system may be relevant in at least some clinically relevant contexts.

To determine whether migration in response to carcinoma-related growth factors was altered upon EMT, we measured steady-state migration of epithelial and mesenchymal cells in response to EGF, HRG-β1, IGF-1, and HGF to activate the ErbB family, IGF1-R, and Met, respectively. Epithelial cells migrated very little as singular cells for all stimuli, whereas single mesenchymal cells migrated robustly in response to select growth factors, notably EGF (Fig. 1D, supplemental Fig. S3). Conversely, epithelial cells moved rapidly within monolayers, at or above the speeds attained by singular mesenchymal cells, even in the absence of exogenous stimuli, with only modest enhancement by some of the growth factors (Fig. 1E, supplemental Fig. S3). Within monolayers, mesenchymal cells exhibited very low cell speeds that were enhanced only slightly by growth factor treatments. These results suggest that the degree of motility in both migratory modes is highly growth factor- and EMT-dependent. Each type of cell responded differentially to growth factor treatments based on its phenotype, indicating distinct processing of growth factor-elicited signals in pre-Twist versus post-Twist cells.

Quantitative Analysis of Growth Factor-Elicited Multiple-Pathway Signaling Network Dynamics

We hypothesized that the changes in Twist-related gene expression could induce alteration of multiple pathways in the signaling network downstream of growth factor cues, leading to the observed EMT-dependent migratory responses. An exciting previous study has reported measurement of more than 1000 biomolecular species at the mRNA, protein, phosphopeptide, or phosphoprotein level in tumors with epithelial and mesenchymal phenotypes to generate annotated molecular network graphs (4), but our focus here is a quantitative analysis of changes in multipathway signaling network activities in comparative manner from before to after EMT induction in a particular cell line with the goal of constructing computational model-based prediction of signaling pathway relationships to motility behavior. To achieve this, we assessed the early activation kinetics of 14 proteins downstream of receptor tyrosine kinase activation (Fig. 3A) in confluent pre-Twist cells and sparse post-Twist cells (at 0, 5, 10, 30, or 60 min after growth factor activation). Measurements of 14 phospho-sites over five time points, five growth factor treatments, two cell lines, and two to three technical replicates resulted in greater than 1800 data-points (Fig. 3C).

Fig. 3.

Early activation profiles of key regulators of cell migration exhibit altered signaling pathway activities upon Twist-induced EMT. A, Simplified schematic of receptor tyrosine kinase-activated signaling network involved in cell migration and the candidates for which phosphorylation was assessed via quantitation. Arrows indicate direct binding of the growth factors used in this study to their respective receptors. Solid lines indicate direct interactions between proteins whereas dashed lines indicate demonstrated indirect interactions. Some candidates have been grouped based on their demonstrated involvement in various biophysical processes of cell migration. This figure is intended as an illustration of the complexity of the signaling network and does not fully account for all components or interactions assessed. Phosphorylation of candidates with green or red “P” symbols has been measured via a quantitative multiplexed bead-based ELISA assay or quantitative Western blot assay, respectively. B, Ratio of basal epithelial and mesenchymal phosphorylation, demonstrating changes in signaling before growth factor stimulation. Data is shown as mean ± S.E., based on error of measurements in each context. n = 2–3 biological replicates per context. C, Sixty-minute time-courses of phosphorylation after stimulation of hMECs with various growth factors. Confluent epithelial cells (solid circle) or sparsely seeded mesenchymal cells (dashed diamond) were lysed at various times after stimulation with EGF (red), HRG (black), IGF-1 (blue), and HGF (green), and subjected to a high-throughput multiplexable bead-based ELISA or quantitative Western blot using antibodies against various phosphorylation sites. Assay wells were loaded with equal mass of protein as assessed by a quantitative bichronic acid assay. Mean fluorescence intensities (MFI) for Western blots were obtained via densitometry. MFI values were normalized to the 0 min Epithelial value within each phosphosite. n = 2–3 biological replicates. Data is shown as mean ± S.E.

Interestingly, total receptor expression levels were not readily correlated with their activity in the context of EMT. For example, although EGFR activation was much greater in epithelial cells (Fig. 3B, 3C) and the total EGF receptor levels were comparable or only slightly higher in epithelial cells (supplemental Fig. S1), Mesenchymal cells were strikingly more responsive to EGF treatment. This is not necessarily surprising, because it is appreciated that receptor expression changes (whether at mRNA or protein level) alone are typically not predictive of associated activity or inhibitor effectiveness; a prominent instance of this is the lack of correlation of EGFR expression in patient tumors with anti-EGFR kinase inhibitor efficacy (e.g. (36)). Thus, assays that focus on receptor expression levels may not by themselves effectively identify key targets for therapeutic intervention.

The resulting activation profiles showed diverse kinetics across individual signals that were growth factor- and EMT state-dependent. Basal phosphorylation levels were dependent on the EMT state, with EGFR, Met, Erk, Src, β-catenin, HSP27, and IRS-1 displaying significantly higher initial phosphorylation in epithelial cells, but Akt, GSK3α/β, PKCδ, PLCγ, and JNK displaying higher phosphorylation levels in mesenchymal cells (Fig. 3B). Dynamic changes in JNK, IRS-1, Src, HSP27, GSK3α/β, and β-catenin phosphorylation after growth factor treatment were cell-state specific and correlated with their initial phosphorylation levels (Fig. 3C). However, activation of PKCδ and PLCγ along with EGFR canonical pathways Erk and Akt were relatively growth factor-dependent and insensitive to EMT state in most cases (Fig. 3C). Visual inspection of signal differences across the diverse treatments and contexts offered little insight into which signals contribute most significantly to the profoundly different EMT-dependent migratory responses. The consequent implication is that cells must quantitatively integrate the activities of multiple signaling pathways to generate robust decisions concerning context- and treatment-dependent migration responses.

Partial Least-Squares Regression Model Accounts for Diverse Motility Behavior Across Phenotypic Modes and Growth Factor Treatments

To understand this quantitative multipathway integration, we applied PLSR to the signaling data set to correlate the intracellular signaling activities to the phenotypic response of the cells and evaluate the ability of intracellular signaling nodes to predict the disparate EMT motility response (Fig. 4). Despite likely differences in which underlying biophysical processes (e.g. lamellipod protrusion, cytoskeletal contraction, and cell and substratum adhesion and de-adhesion) may be rate-limiting for the epithelial and mesenchymal modes of migration, we found that a single model comprised of two principal components (quantitative combinations of the key signaling node phosphorylation states, for both background activity {‘T0′} and growth factor stimulation-induced activity integrated over the 1-hour time period {‘Int’}) was able to account in a unified manner for migration across both EMT states and all four growth factor conditions. The Loadings and Scores plots shown in Fig. 4A for this model illustrate some interesting insights. Phosphorylation state of some signaling nodes (GSK3α/β, PLCγ, PKCδ, JNK, Akt) are more closely associated with mesenchymal migration, whereas some others (β-catenin, IRS-1, Src) are more closely associated with epithelial migration; for yet some other signaling nodes (Erk, HSP27) growth factor-stimulated activities are more closely associated with mesenchymal migration whereas their constitutive activities are more closely associated with epithelial migration. Thus, there is no individual signal that is uniquely associated with either migration mode, nor is solely stimulated signaling the crucial determining feature. Nonetheless, the quantitative combination of these signals is able to account for the entire scope of behavior comprehensively, as Fig. 4B demonstrates leave-one-out cross-validation of this two-component model. Basal phospho-signal levels were found to be important for prediction of motility in conjunction with the dynamic treatment-induced phospho-signal levels (R2/Q2 of 0.79/0.58 for the full model, 0.78/0.09 without initial time-points), demonstrating that the Twist-induced network modulation comprises both growth factor-dependent (stimulus-induced) and growth factor-independent (background, before stimulation) effects; both types of effects are integrated seamlessly into the PLSR model, as has been similarly found in previous work (27).

Fig. 4.

A multivariate partial least squares regression model captures signaling metrics contributing most to the prediction of both epithelial and mesenchymal cells. A PLSR model has been constructed using the initial phosphorylation levels and those integrated over 60 min of the 14 signals described in Fig. 3 across serum-free, EGF, HRG, IGF, and HGF treatments. A, Projection of loadings (left) and scores (right) onto the first two principal components. Loadings of individual signaling metrics (Int = integral of phosphorylation; T0 = initial phosphorylation) are plotted in black. Loading of cell speed metric is plotted in red. Scores of each growth factor treatment are plotted black for epithelial and red for mesenchymal cells. B, Leave-one-out cross-validation of the PLSR model with cell speeds predicted by the two principal component model versus experimentally measured cell speeds.

To test the predictive capability of our PLSR model with respect to its central multivariate “signal-response” relationship, we undertook assessment in a new growth factor condition, PDGF-BB, for both epithelial monolayer and mesenchymal sparse conditions. PDGF is increasingly appreciated as an important growth factor in many breast cancer situations (37). As before, the transient signaling dynamics were measured (Fig. 5A); now they were simply inserted into the previously trained model for a priori prediction of consequent migration behavior in the new PDGF pre- and post-Twist cases. This model successfully predicted the migratory behavior (Figs. 5B, 5C). This test is especially stringent, as migration speeds in both the epithelial and mesenchymal contexts were outside those previously observed within each phenotypic context, requiring the model to generate extrapolative predictions. We also tested whether inhibition of PKCδ would compromise mesenchymal motility more than epithelial motility, because of its stronger predicted association with the former (Fig. 4A) as well as our previous findings of its involvement in fibroblast motility (38). Treatment of pre-Twist and post-Twist cells with the small molecule inhibitor rottlerin showed, as predicted, that the motility behavior of the latter was more strongly reduced in dose-response manner (supplemental Fig. S4); similar effects were found for the BT549 breast tumor cell line, which also exhibits high expression of Twist (supplemental Fig. S4). Although this pharmacological intervention was successful, technical limitations precluded us from confirming the result using an independent approach. Prediction tests such as this one using small molecule inhibitors, or RNA interference methods to knock down individual signaling nodes can be undertaken, as we emphasize that our PLSR model is correlative rather than causal in nature. However, the prevalence of off-target effects and/or transcriptional compensation mechanisms renders the effectiveness of a single-target perspective problematic without narrowly restricting study conditions (see Discussion).

Fig. 5.

PDGF-stimulated epithelial and mesenchymal cell speeds are predicted A, Sixty-min time courses of epithelial cells in monolayer (solid purple) and sparse mesenchymal cells (dashed purple) stimulated with 50 ng/ml PDGF-BB for the 14 signals. Already presented time courses are provided in gray for comparison. B, Cell speeds of epithelial and mesenchymal cells under stimulation of PDGF (purple) quantified from the sparse (top) or monolayer (bottom) migration tracks as described previously. Cell speeds under stimulation of other growth factors (gray) are plotted as comparison. All data is shown as mean ± S.E. n = 233–310 for monolayer migration and n = 15–70 for sparse migration. C, Prediction of epithelial monolayer and mesenchymal sparse cell speeds stimulated by PDGF by the PLSR model trained only on the unstimulated, EGF, HRG, IGF, and HGF data set. Predicted PDGF cell speeds (purple) are plotted against experimentally measured cell speeds.

Node-to-Node Correlation Topology Model Reveals Quantitatively Different “Information Flow” Between Epithelial And Mesenchymal Modes

Based on the inability of receptor expression to explain changes in growth factor responsiveness, striking changes in observed signaling, and the retained ability of cells to migrate in all contexts, we hypothesized that differences in downstream signaling might arise from quantitatively different multipathway integration, despite retention of the qualitative network topology. This is in addition to likely signaling-independent changes. In order to investigate Twist-induced differences in node communication downstream of receptor signaling, correlative topology modeling was performed for the epithelial and mesenchymal contexts as described in the Experimental Procedures section. Separate network topology inference models were constructed for each EMT cell state.

Results using the Storey method are shown in Fig. 6, whereas the Bonferroni and Benjamini methods' results are shown in supplemental Fig. S5. In the context of network inference, each significant correlation value represents one undirected edge in the inferred network. In the epithelial state, using the Storey method with a false discovery rate of 0.08 (p < 0.046) provides for 12 significant correlation values, resulting in an estimated 0.08 × 12 edges ≈ 1 false positive edge. In the mesenchymal state, using the Storey method with a false discovery rate of 0.11 (p < 0.028) provides for 10 significant correlation values, resulting in an estimated 0.11 × 10 edges ≈ 1 false positive edge. Not surprisingly, the greatest number of node-to-node influence arcs were found in the Storey models, and the smallest number in the Bonferroni models, given the more conservative nature of the latter algorithm for assigning significance. Accordingly, all arcs found in the Bonferroni models were found in the corresponding Benjamini models and all arcs found in the Benjamini models were found in the corresponding Storey models.

Fig. 6.

Correlative topological modeling, comparing epithelial ( Edges between phosphorylation sites indicate statistically significant positive (black) or negative (red) Pearson correlation (Storey multiple hypothesis correction, ∼1 false positive edge). Models using more stringent multiple hypothesis testing are contained in supplemental Fig. S5. Edge end annotation indicates literature evidence for detected correlation (listed in supplemental Table S1), including direct phosphosite-specific or pathway-level evidence (arrowheads), protein-level evidence (diamond ends) and complex-level evidence (dot ends). Nodes with red font indicate phosphosites that signal for inhibition or degradation of the protein. Dashed edges have the highest first-order partial correlation p value within each three-node triplet.

First-order partial correlation determines if the correlation between two nodes may be explained because of their mutual correlation with a third node (39). First-order Pearson partial correlations were calculated for all three-node triplets present in the networks. A three-node triplet occurs when three nodes, A, B, and C form a complete subnetwork, whereby significant pairwise correlation exists between nodes A and B, A and C, and B and C. Using the Storey method's resultant networks, in the epithelial network a three-node triplet exists between JNK, Erk1/2, and IRS-1. In the mesenchymal network, two three-node triplets exist: Akt, GSK3α/β, and HSP27; and PKCδ, GSK3α/β, and IRS-1. For each triplet, the partial correlation edge with the highest p value is shown as a dashed edge in Fig. 6. These results suggest that the dashed edge in each triplet may exist because of mutual correlation with the third node.

Striking differences in the set of significant edges suggests network modulation occurs upon Twist-induced EMT, wherein information is processed via changes in the influence signaling pathway nodes have upon one another. Key similarities and differences will be noted, in context of available literature information, in the Discussion section.

DISCUSSION

Our objective in this report has been to investigate how activities in multiple signaling pathways downstream of a range of receptor tyrosine kinases are changed between pre-EMT condition and post-EMT condition, especially with respect to their contributions to regulation of cell motility. We emphasize that we have not aimed to investigate signaling pathway activities involved in or responsible for the act of EMT induction per se, for that question has been addressed with great effectiveness by several laboratories during the past decade (e.g. (1–6)). We also note that our analysis focuses on EMT induced by ectopic expression of the Twist1transcription factor, for two reasons: first, Twist expression is strongly implicated in clinical tumor biology; and, second, it may represent a relatively tightly defined dysregulation, compared with extracellular inducers such as TGFβ and TNFα, which alter expression of multiple EMT-associated transcription factors along with Twist (7).

We have found although both epithelial and mesenchymal-like cells possess migratory potential, their growth factor-elicited behavior is substantively distinct with respect to contexts under which vigorous motility is exhibited. Epithelial cells are predominantly motile only within confluent monolayers in which cell-cell contacts are maintained (consistent with previous findings (32, 33), whereas mesenchymal-like cells are motile mainly as individual cells and exhibit this best when sparsely distributed (Fig. 1). The responsiveness to any particular growth factor depends on whether the cells are in an epithelial or mesenchymal-like state. With respect to EMT-associated alterations in growth factor-induced signaling network activities downstream of the stimuli/cues that might be critically involved in disparate motility responses, we showed that quantitative and dynamic properties of numerous phosphoprotein signaling nodes were comparatively modulated from pre- to post-Twist conditions across the different growth factor treatments (Fig. 3). PLSR analysis successfully demonstrated that multipathway signaling information can be quantitatively integrated to account for motility behavior across all observed contexts (Fig. 4)—and can even predict a priori the motility responses in both epithelial and mesenchymal situations to treatment by an additional growth factor, PDGF (Fig. 5). Finally, we analyzed each EMT-condition separately in order to identify differences in signaling. Correlative topological modeling suggested Twist-dependant differences in terms of a network-level explanation for disparate motility responses (Fig. 6). Specifically, we propose a concept of “operational rewiring,” in which the dominance of particular nodes on motility is altered by quantitative modulation of node-to-node influences.

We can provide validation of this correlation network model results by comparison to previous literature reports, toward biological relevance of our inference insights, illustrated by some particular examples here along with a larger set of information in supplemental Table S1. For example, three arcs are ascertained to be invariant between epithelial and mesenchymal states: Akt-GSK3, Akt-HSP27, and Erk-IRS1. Although the role of GSK3 as a substrate downstream of Akt is strongly established (e.g. (40)), emerging evidence indicates that HSP27 helps enhance Akt activity (41, 42); similarly, recent reports suggest that serine phosphorylation on IRS1 is mediated by p70S6K in a manner regulated by ERK (43, 44). As another example, in the epithelial cell state Src phosphorylation on tyrosine 416 is closely associated with β-catenin phosphorylation, whereas in the mesenchymal cell state the same Src phosphosite is more strongly associated with PLCγ; the predominance of an activity correlation between Src and β-catenin in the epithelial state is recognized (e.g. (45)), and mechanisms for positive influence of Src activity on PLCγ-mediated membrane protrusion in disperse, Mesenchymal state cells are being elucidated (46, 47). One final example worth emphasizing here is that the influence of PKCδ changes more dramatically between pre-Twist and post-Twist conditions than that for any of the other nodes in our network: it loses a negative influence on Akt and gains positive influences on β-catenin, GSK3α/β, and IRS-1, and thus has a stronger presence in network dynamics for the mesenchymal state relative to the epithelial state. This change is consistent with our findings that inhibiting PKCδ more sensitively affects mesenchymal motility than epithelial (supplemental Fig. S4). Previous reports on the role of this kinase on tumor cell motility have been conflicting, with both migration-enhancing and migration-diminishing effects observed depending on the cell type and environmental context (48, 49). Direct comparison of how EMT alters influences of PKCδ on downstream activities has not been previously investigated, although a recent report implicates PKCδ in Wnt-induced invasive motility of colorectal tumor cells via interaction with a β-catenin/Dishevelled complex (50).

Because of the highly multivariate nature of the relationship demonstrated here (Figs. 3 and 4) between signaling activities and migration behavior, and the likelihood that molecular perturbation of any node (whether by small molecule inhibitor or RNA interference) will generate substantial cross-talk effects among pathways, aspiring to narrowly test specific model predictions for individual signals across the diverse conditions is fraught with intricate complication—not merely technical but conceptual. Although under narrow circumstances qualitative trends may be found consistent with a model prediction for the effect of modulating a particular individual signal (e.g. (51–53)), successful prediction of intervention effects over a range of environmental conditions typically requires capture of key cross-talk consequences across multiple pathways (e.g. (54)).