Integrated analysis reveals that STAT3 is central to the crosstalk between HER/ErbB receptor signaling pathways in human mammary epithelial cells

Abstract

Human epidermal growth factor receptors (HER, also known as ErbB) drive cellular proliferation, pro-survival and stress responses by activating several downstream kinases, in particular ERK, p38 MAPK, JNK (SAPK), the PI3K/AKT, as well as various transcriptional regulators such as STAT3. When co-expressed, the first three members of HER family (HER1–3) can form homo- and hetero-dimers, and there is considerable evidence suggesting that the receptor dimers differentially activate intracellular signaling pathways. To better understand the interactions in this system, we pursued multi-factorial experiments where HER dimerization patterns and signaling pathways were rationally perturbed. We measured the activation of HER1–3 receptors and of the sentinel signaling proteins ERK, AKT, p38 MAPK, JNK, STAT3 as a function of time in a panel of human mammary epithelial (HME) cells expressing different levels of HER1–3 stimulated with various ligand combinations. We hypothesized that the HER dimerization pattern is a better predictor of downstream signaling than the total receptor activation levels. We validated this hypothesis using a combination of model-based analysis to quantify the HER dimerization patterns, and by clustering the activation data in multiple ways to confirm that the HER receptor dimer is a better predictor of the signaling through p38 MAPK, ERK and AKT pathways than the total HER receptor expression and activation levels. We then pursued combinatorial inhibition studies to identify the causal regulatory interactions between sentinel signaling proteins. Quantitative analysis of the collected data using the modular response analysis (MRA) and its Bayesian Variable Selection Algorithm (BVSA) version allowed us to obtain a consensus regulatory interaction model, which revealed that STAT3 occupies a central role in the crosstalk between the studied pathways in HME cells. Results of the BVSA/MRA and cluster analysis were in agreement with each other.

---

Introduction

Because of their fundamental regulatory role in physiological processes such as cell proliferation, survival, and migration,¹ Human Epidermal growth factor Receptor (HER, also known as ErbB) signaling has been extensively studied over the years. HER receptors belong to the receptor tyrosine kinase superfamily, and the first three members of this receptor family (EGFR/HER1, HER2 and HER3) are often co-expressed, and they play important roles in carcinogenesis.^2,3 HER receptors are highly homologous and their activation occurs through similar biochemical steps, although the individual receptors possess certain distinct properties. For example, HER2 receptors have no known ligand, and HER3 lacks intrinsic tyrosine kinase activity.^4,5 Upon ligand binding, receptors undergo conformation changes that favor the formation or stabilization of receptor hetero- and homo-dimers,^6,7 and almost every possible pairwise combination of HER dimers has been reported.⁸ Dimer formation patterns are affected by ligand availability,⁹ receptor expression levels,¹⁰ and kinetic dimerization hierarchy.^11–13 Each HER receptor appears to have a characteristic repertoire of adaptor proteins depending on its dimerization partner.¹⁴ Which receptor–adaptor protein complexes are formed determines differential downstream signal transduction, and can trigger diverse biological responses.^15,16 Thus, the types of receptor homo- and hetero-dimers that are formed, and the signal transduction pathways that are activated depend upon the specifics of the system and treatment conditions. This severely complicates investigation of the HER signaling pathways and the HER-initiated cellular responses.^17–20

Biomedical and biomolecular studies have shown that inhibition of HER dimerization and HER-mediated signaling can be an extremely effective therapeutic strategy.²¹ The role of HER receptors seems to be different in different cell lines and under different conditions. In HER2 positive breast cancers, heterodimer interactions between HER2 and its partners are often constitutively activated and their disruption has proven to be an effective means for inhibiting HER2-mediated aberrant responses.^1,22 HER3 overexpression has also been shown to correlate with poor prognosis in epithelial cancers.^23,24 Activation of pro-survival responses through the HER3 pathway can also lead to drug resistance in cancer treatment.^25,26 Further, interactions between members of the HER family and signaling redundancies can also contribute to drug resistance when specific inhibitors are used against a single HER receptor and the importance of HER3–HER1 dimers in drug resistance is increasingly recognized.^21,25,27 Thus, the use of combination therapeutics that target multiple members of the HER family of receptors is advocated as an effective treatment strategy for a variety of tumors.^25,26,28,29

Receptor tyrosine kinases drive cellular responses via the activation of several distinct intracellular signaling pathways.^30,31 HER receptors activate pathways of several downstream kinases, in particular ERK, p38 MAPK, JNK (SAPK), the PI3K/AKT signaling pathways, as well as various transcriptional regulators such as STAT3.^1,15,32–42 These are the prominent sentinel proteins of the cellular processes associated with cell proliferation, pro-survival and stress-response pathways of cells. These inter-related signaling pathways can be stimulated to varying degrees by different receptor–adaptor protein complexes.¹⁴ The abundance of these complexes in turn is a function of the relative abundances of the various HER dimer types in the cell.

Many stimuli, including growth factors, cytokines, transforming agents, and carcinogens activate the ERK (extracellular signal-regulated kinase) pathway. Sustained ERK activation can have significant functional consequences.^43,44 For example, constitutive or prolonged ERK activation through HER receptors contributes to speedy cell migration⁴⁵ and enhances cell transformation and resistance to apoptosis.⁴⁶ The serine/threonine kinase AKT promotes cellular survival by blocking apoptosis, which it achieves by binding and regulating many downstream effectors; e.g., NF-κB, Bcl-2 family proteins, and MDM2.⁴⁷ AKT also induces protein synthesis and therefore plays a role in tissue growth. AKT has been implicated in many types of cancers because of its possible role in angiogenesis and tumor development.^48,49 PI3K is a well-studied regulator of cellular growth⁵⁰ and is known to be activated by HER3, especially in HER2 over-expressing cancer cells.^51,52 HER3 in particular has a unique ability to efficiently engage the pro-survival PI3K/AKT pathway.^53–55

ERK, p38 MAPK and JNK are the three major subgroups of the MAPK family.⁵⁶ In contrast to the ERKs, which are mainly mitogen activated proliferation/differentiation factors, JNK and p38 MAPK are stress-activated proteins involved in apoptosis/survival. The p38 MAPK kinase regulates EGFR down-regulation, survival, and cellular migration,^38,57,58 and it can suppress cell proliferation by antagonizing the JNK/c-JUN pathway.³⁹ Although all MAPKs can be activated by HER receptors in different cellular systems, p38 MAPK and JNK phosphorylation may also occur in an EGFR-independent manner.⁵⁹ STAT signaling promotes RAS/RAF/ERK signaling indirectly.⁶⁰ Activation of MAPKs leads to the activation of STAT-family transcription factors, which promotes their translocation to the nucleus^32,61 where they elevate the transcription of cell proliferation related genes.^32,60

Upon HER receptor stimulated activation of cell proliferation, pro-survival and stress-response signaling pathways, cellular responses are regulated collectively by these prominent sentinel proteins.^30,31 Identification of how the regulatory interactions between sentinel proteins are affected by the unavoidable simultaneous activation of multiple signaling pathways, i.e., pathway crosstalk, has been very challenging and requires multi-factorial experiments where HER dimerization patterns and crosstalk between signaling pathways are altered. Towards this goal, we have measured HER1–3 activation, and the activation of key signaling proteins ERK, AKT, p38 MAPK, JNK, STAT3 as a function of time in a panel of human mammary epithelial (HME) cells expressing different levels of HER1–3 stimulated with various ligand combinations, e.g., when the signaling pathways were activated differentially in cells with same cellular background but with varying HER expression levels. Conditions covered in these studies resulted in large differences in receptor dimerization,¹³ and in the activation levels of investigated signaling proteins. These time-series measurements were supplemented with combinatorial measurements where we determined how the inhibition of one of the sentinel proteins affected the activation levels of the other sentinels. This data was used to identify the regulatory interactions between the studied sentinel proteins using the modular response analysis method.

Experimental materials and methods

General considerations

There is considerable evidence that the HER dimerization pattern is an important determinant of the consequences of HER signaling,^62,63 which suggests that activation of the downstream signaling pathways would be functions of receptor dimer formation and activation.¹⁸ In this regard, we hypothesized that the HER dimerization pattern is a better predictor of downstream signaling than total receptor activation levels. Validation of this hypothesis and quantification of the relative potency of the HER dimer types for activating distinct downstream pathways would enable rational perturbation of the signaling patterns within a cell.

Due to the redundancy in downstream pathway activation, and the extensive crosstalk between signaling pathways, establishing the relationship between HER dimers and signaling requires: (1) the ability to quantify the HER dimerization patterns, and (2) multi-factorial experiments where HER dimerization patterns and signaling crosstalk is altered. Direct isolation and measurement of the receptor dimers is challenging, and one is largely restricted to the measurement of the total receptor activation levels. However, as we have recently shown that the dimer activation levels (point #1 above) can be derived from receptor activation data using an integrated experimental-computational approach.¹³ To satisfy the second requirement, we measured HER1–3 activation, and the activation of key signaling proteins ERK, AKT, p38 MAPK, JNK, STAT3 as a function of time in a panel of HME cells expressing different levels of HER1–3 stimulated with various ligand combinations. Conditions in these studies result in large differences in receptor dimerization,¹³ and in the activation levels of investigated signaling proteins (cf., Results section below). Ability to obtain internally consistent and complementary data sets in a panel of HME cells with the same cellular background is a novel aspect of our study and it was critical for deriving statistically significant regulatory interactions from the experimental data.

Cell lines

The HME cell line 184A1L5 used in this study was originally provided by Martha Stampfer (Lawrence Berkeley Laboratory, Berkeley, CA). It expresses approximately 200 000 molecules of EGFR/ErbB1/HER1.^64,65 Cells were maintained at 37 °C in 5% CO₂ air in DFCI-1 medium supplemented with 12.5 ng ml⁻¹ EGF (PeproTech, Rocky Hill, NJ) as described previously.⁶⁶ At about 70–80% cell confluency, DFCI-1 medium was replaced with bicarbonate-free DFHB minimal medium lacking EGF, plus 0.1% bovine serum albumin. Cells were then brought to quiescence for 12–18 hours before treatment. We note that, unlike many model cell lines with limited physiological responses to growth factors, HME cells are an excellent model system for studying the properties of the HER system because, like many epithelium derived cell types, HME cells require EGFR activation for proper proliferation and migration responses.⁶⁷ The 184A1 strain of HME cells that was used to develop our clone library retains the EGFR-dependent regulatory machinery of the primary cell type from which it was derived.⁶⁸ This parental cell line expresses high levels of EGFR and relatively low copy numbers of HER2 and HER3. It is labeled the HER2−/3− cell line, and it was used to create the HER2 and HER3 overexpressing cell lines by retroviral transduction.⁶⁵

Briefly, the retroviral MFG vector containing HER2 was constructed and transfected into the ψCRIP packaging cell line.^69–71 The transfected clones with highest level of HER2 expression were expanded to produce viral supernatant. This supernatant was then used to infect the HME cells. Individual clones of retrovirus transduced HME cells were isolated using cloning rings. The degree of HER2 overexpression was determined by immunofluorescence, flow cytometry, and equilibrium-binding studies using labeled Fab fragments. HER3 expressing cells were derived in a similar retrovirus-based strategy. Retroviral pBM-IRESpuro vector containing HER3 and the puromycin resistance gene was transfected into the Phoenix packaging cells.⁷² Viral supernatant was collected from transiently transfected Phoenix cells. The supernatant was then used to infect the parental HME cells and the HER2 expressing cells (discussed in the previous paragraph) to create the HER3 expressing and HER2–HER3 co-expressing cell lines, respectively (Table 1). Growth medium for the HER2 and/or HER3 expressing cells was the same as that of the parental HME cells except for addition of the antibiotic G418 (250 μg ml⁻¹; Invitrogen, Carlsbad, CA) and/or puromycin (2 μg ml⁻¹; Sigma, St. Louis, MO) to ensure selection for HER2 and HER3, respectively. Utilized cell lines are labeled and referred to as R2a3b where a or b = +/− indicates that the cell line expresses/do not express (i.e., very low endogenous expression) the type 2 and type 3 HER receptors, respectively. We note that these cell lines were labeled with their clone names in our earlier studies^13,65 as R2−3− (parental, or Par for short), R2+3− (24H), R2−3+ (B5), and R2+3+ (D20) and these clone names are used interchangeably in this report as well.

Table 1 Human mammary epithelial cell library

Cell lines^a			Receptor abundance (1000 receptors per cell)
Label	Clone tag	Description (HER1/2/3)	EGFR	HER2	HER3
a Cell labels and clone tags are used interchangeably in the text.
HER2−3−	Par	Parental, very low HER2/3 (+/−/−)	200	30	2
HER2+3− or HER2+	24H	HER2 expressing (+/++/−)	137	600	2
HER2−3+ or HER3+	B5	HER3 expressing (+/−/+)	189	30	28
HER2+3+	D20	HER2/3 expressing (+/++/+)	86	644	29

---

Signal initiation

To activate the cells, ligands were added to the culture medium, and cells were incubated at 37 °C for fixed amounts of time. Stimulation of EGFR/HER1 with 12 ng ml⁻¹ EGF, of HER3 with 40 ng ml⁻¹ HRG, as well as the co-stimulation of both receptors via simultaneous addition of both ligands was studied. For simplicity, in the remainder of the paper, these ligand treatment conditions will be labeled and referred to with +E, +H, and +EH where E and H stand for EGF and HRG, respectively. No-ligand measurements provided the negative control as 0 min time point. We note that same ligand doses were used in the ELISA, Western blot and inhibition experiments described below to obtain consistent data sets.

Receptor activation measurements using ELISA

Phosphorylation levels of the receptors were measured at 5, 10, 20, 30, 60, 90 and 120 min after ligand stimulation. Data for the phosphorylation levels of HER1–3 receptors were collected by us earlier.^13,65 These measurements used the R&D DuoSet IC ELISA kits (R&D Systems Inc., Minneapolis, MN), and detailed descriptions can be found in our earlier studies.^13,17,65 Briefly, cells were washed 3× with ice cold PBS and incubated at room temperature for one minute to allow surface receptor dephosphorylation. After another round of cold PBS washing, cells were solubilized with ice cold lysis buffer (1% NP-40, 20 mM pH 8.0 Tris buffer, 137 mM NaCl, 10% glycerol, 2 mM EDTA, supplemented with 1 mM heat activated sodium orthovanadate and 1% protease inhibitor cocktail III; Calbiochem, La Jolla, CA) for 20 min. Lysates were centrifuged at 13 000 rpm for 10 min at 4 °C, and the supernatants were transferred into fresh microtubes. Obtained cell lysates were either analyzed immediately or stored at a −80 °C freezer until needed.

ERK and AKT activation measurements using ELISA

Data for the phosphorylation levels of the kinases ERK and AKT were quantified in ELISA experiments by us earlier and has already been communicated.¹⁸ Briefly, R&D DuoSet IC ELISA kits were used and the manufacture's protocol was followed in these measurements. Before each ELISA assay, protein concentration of cell lysate was measured using the Bicinchoninic Acid protein quantitation kit (Sigma, St. Louis, MO). ERK and AKT phosphorylation levels were measured at 5, 10, 20, 30, 45, 60, 90 and 120 min after ligand stimulation. Further details of these experiments can be found in Zhang et al.¹⁸

Western blot analysis of p38 MAPK, JNK and STAT3

Cells were grown to near confluency and were starved by replacing regular DFCI-1 culture medium with DFCI-1 medium lacking all supplements but 0.1% bovine serum albumin for 12–18 hours. Cells were stimulated at 37 °C by adding the appropriate concentration of ligands for the desired durations. Afterwards, cells were washed twice with cold PBS and cell lysates were collected with a scraper. Cells were lysed using RIPA buffer containing protein protease and phosphatase inhibitors on ice for 10 to 30 min. Lysates were centrifuged at 12 000 rpm for 10 min at 4 °C, and only the supernatants were kept in new microtubes. Cell lysates were either analyzed immediately or stored at −80 °C until needed. Total protein concentration of cell lysate was determined by the Bicinchoninic Acid protein quantitation kit (Sigma, St. Louis, MO), and same total protein amount (25 μg) was loaded onto the SDS-PAGE gel. NuPAGE Bis-Tris precast 4–12% gradient gels with 15 wells were used. Proteins were transferred from SDS-PAGE gel to PVDF membrane using XCell SureLock™ Mini-Cell Electrophoresis System (Invitrogen). PVDF membrane were developed using ECL reagents and Western blot images were captured by Lumi-Imager F1 device and analyzed using the LumiAnalyst 3.1 software.

JNK, STAT3, and p38 MAPK phosphorylation levels were measured at 0, 15, 30, 60, 120, 240 and 480 min after ligand stimulation. For each experiment, at least two independent measurements were performed for every treatment condition. To be able to scale the data collected in different gels, we allocated two of the lanes to “standard” measurements, which allowed us to normalize the measurements across gels. As having common positive controls is necessary for reliable quantification, results for the stimulation with both EGF and HRG ligands at 15 min time point for the parental (R2−3−) and B5 (R2−3+) cells were measured in each gel and then used as the standard data for the gel-to-gel normalization. Additionally, to further improve the comparison of the results for the four cell-lines, each gel was designed to measure the activation of two cell lines. This allowed us to do quality check across experimental time series datasets.

Inhibition studies

To further investigate the crosstalk between the signaling pathways, we selectively inhibited the studied proteins one at a time and measured how the particular inhibitions affect the activities of the studied sentinels. Cell lysates for inhibitor-treated cell assays were collected as described above for the no-inhibition studies except that during cell starvation inhibitors were added to the medium and incubated for certain time. U0126 monoethanolate (ERK1/2 inhibitor), SB 203580 (p38 MAPK inhibitor), and AKT inhibitor VIII were added to starvation medium 2 hours before ligand stimulation all at 20 μM concentration. LY294002 (hydrochloride; PI3K inhibitor) was added to starvation medium 1 hour before ligand stimulation at a concentration of 10 μM. JNK inhibitor SP600125 and STAT3 inhibitor VIII were added to starvation medium 18 hours before ligand stimulation at 40 μM concentration. Phosphorylation levels of the studied proteins were measured 15 min after ligand stimulation in these inhibition studies using Western blots for p38 MAPK, JNK and STAT3, and ELISAs for ERK and AKT as described in the previous sections.

Reagents

Epidermal growth factor (EGF, human recombinant) and heregulin-β1 (HRG, human recombinant) were purchased from Peprotech (Rocky Hill, NJ). Primary Western blot antibodies against total p38 MAPK (#9212), SAPK/JNK (#9258), and STAT3 (#9139), and against phosphorylated p38 MAPK (Thr180, Tyr182, #9215), SAPK/JNK (Thr183/Tyr185, #9255), and STAT3 (Tyr705, #9145) were purchased from Cell Signaling Technology (Danvers, MA). Primary antibody for β-actin (A2228) was purchased from Sigma. Anti-mouse IgG-peroxidase secondary antibody (A4416) was purchased from Sigma (St. Louis, MO), and anti-rabbit IgG HRP-linked secondary antibody (7074) was purchased from Cell Signaling Technology (Danvers, MA). U0126 monoethanolate (ERK1/2 inhibitor), SB 203850 (p38 MAPK inhibitor), AKT inhibitor VIII, and STAT3 inhibitor VIII were purchased from Calbiochem (EMD Millipore, Billerica, MA). PI3K inhibitor LY294002 (L9908) hydrochloride was purchased from Sigma (St. Louis, MO). JNK inhibitor SP600125 was purchased from LC Laboratories (Woburn, MA). ELISA assay kits were purchased from R&D Systems Inc. (Minneapolis, MN). G418 was purchased from Invitrogen (Carlsbad, CA), and puromycin was purchased from Sigma (St. Louis, MO). Cell culture supplements, SDS-PAGE gels, and PVDF membranes were purchased from Invitrogen (Carlsbad, CA). The high sensitivity peroxidase-based Electrochemiluminescence (ECL) immunoassay substrate for Western blot was purchased from Thermo Scientific (Pittsburg, PA). All the other reagents were purchased from Sigma (St. Louis, MO) unless otherwise indicated.

Computational methods

HER1–3 dimerization

One of the aims of this study is to show that cellular signaling initiated by the HER receptors can be better characterized by the dimerization patterns between the HER family members. In other words, not the receptor phosphorylation per se but which receptor dimers are formed and activated is the determinant of the information flow through HER signaling pathways. This requires the quantification of the HER receptor dimer activation profiles. Although this information can in principle be obtained in co-IP type experiments, such experimental studies are technically challenging and do not provide the level of quantitation required for our analysis. We therefore measured the receptor phosphorylation instead, and used a model-based analysis to extract the HER receptor dimerization profiles from the collected experimental data. Our modeling approach is described in detail elsewhere.^13,73 Therefore, it will only be briefly summarized here.

Model used to analyze the receptor activation considers ligand binding to receptors, receptor dimerization and trafficking through endocytosis and degradation. All possible dimer combinations are included except HER3–HER3 pair, which is excluded because this dimer has very low kinase activity and is not formed at significant levels.^5,26,74–78 Although HER3 has been shown to have some kinase activity, as a kinase it is ∼1000 times less efficient than EGFR/HER1.⁷⁹ Therefore, when other HER types are present, HER3 can be considered as kinase deficient for practical reasons. Additionally, even at its highest expression, HER3 level in the cells is about 30 000 receptors, which is about an order of magnitude less than HER1 and HER2 in our HME cells (Table 1) as well as in many other epithelial cell types.¹³ Thus, even if they are formed, the number of HER3–HER3 receptor pairs will be low due to the availability of large number of HER1 or HER2 receptors for partnering and the contribution of HER3–HER3 pair to HER3 phosphorylation would be minute because of its low kinase activity. For these reasons, HER3–HER3 pair was not included in our model to avoid addition of extra parameters unnecessarily, i.e., to avoid parameter overestimation.

Model parameters were fitted to the experimental data to obtain its kinetic parameters. Parameter sets that explain the data the best were then used to compute the receptor dimerization and activation profiles. One key feature of our model is that total phosphorylation of an HER receptor type is expressed as a sum of the contributions from all possible dimer types. For example, phosphorylation of EGFR/HER1, pR1, can be expressed as the sum of contribution of phosphorylated HER1 to cellular signaling when it is in dimer complexes with HER1 (pR1r1), with HER2 (pR1r2), and with HER3 (pR1r3): pR1 = pR1r1 + pR1r2 + pR1r3. We note that our analysis assumes that the phosphorylated forms of the proteins are active. Therefore, in the remainder of the paper, we use the terms “active” and “phosphorylated” receptors interchangeably.

Correlation analysis

Above described experimental studies provided us the phosphorylation data as time series. Our analysis combined our earlier data (for HER1–3, ERK and AKT) with the new data (for p38 MAPK, JNK, STAT3) that were collected as part of this study. Because of the limited number of lanes in SDS-PAGE gels, we were not able to measure the protein activation levels at all the time points that were used in our earlier ELISA experiments (cf. Experimental methods section). We therefore used the 15 min time point in the Western blot experiments and interpolated the ERK and AKT ELISA results for this time point by averaging their measured activation levels at 10 and 20 minutes.

A set of 12 time series data were collected for every protein; 3 ligand treatment conditions (+E, +H, and +EH) in 4 cell lines (R2−3−, R2−3+, R2+3−, and R2+3+). Obtained time series data was used to identify the correlations between the activation patterns of the proteins at different levels, and the cluster analysis was used to group the proteins. The analysis was performed using MATLAB's clustering software. In the cluster analysis, 1 − |Corr(X,Y)| values were used as the distance measure where Corr(X,Y) is the correlation between protein X and protein Y over the included conditions. We preferred to use the absolute value of the correlation in the distance measure to allow for the possibility that the protein activation patterns may be negatively correlated, which too imply functional relationships. However, significant correlations in the investigated data sets were all positive so reached conclusions would have been the same if a 1 − Corr(X,Y) distance criteria was used. Hierarchical cluster analysis utilized the average distance criteria in assigning the clusters and the connection values. Corr computation only included the common time points, and correlations were calculated using the protein activation values at the same time points.

Availability of data for different cell lines and different ligand treatment conditions allowed us to pursue the cluster analysis at multiple levels: first, we combined the data for four cell lines and three ligand treatment cases and created an “all-data” set for the proteins. This set made it possible to investigate the overall correlations between the investigated proteins at a higher level. We then divided up the data to perform the analysis on a “cell-type” and “treatment-type” manner. In the former set, data about a protein for all ligand treatment cases were combined for each cell line to form cell-line based datasets. Similarly, in the latter set, data for the four cell lines for a particular ligand treatment condition was combined.

Additionally, the analysis was pursued by investigating the correlation between the downstream sentinels and receptors by (1) using the activation data for the receptors directly (i.e., without considering dimers), and (2) using the dimer contributions to the receptor phosphorylations. Comparison of these two approaches allowed us to validate our hypothesis. We note that, as it consists of closely related cell types and overlapping time spans, data collected and analyzed in this study was internally cohesive, and this has allowed for a systematic investigation.

Modular response analysis

Correlation analysis described above does not provide information about the causal regulatory interactions between studied proteins. Modular Response Analysis^80–82 (MRA) is a powerful method to reverse engineer a system to infer the existing regulatory interactions and their causality in a quantitative manner. In MRA, one starts with perturbing the system to measure how varying a module impacts the response of the other modules in the system. Observed variations are then used to derive the local response matrix elements r_ij, which quantifies the sensitivity of module i (i.e., activation of sentinel protein i in our case) to changes in module j of the network system, r_ij = ∂x_i/∂x_j provided that activity of all other nodes are kept constant where x_i denotes the activity of module i.^81,82 When multiple interactions may be involved, changes in the activities can be characterized with the global response matrix with elements R_ij, which corresponds to the change in component i, x_i, in response to a perturbation in component j, P_j, as R_ij = ∂ ln x_i/∂ ln P_j. For small enough variations R_ij may be approximated as R_ij ≅ 2(x^(j)_i − x⁽⁰⁾_i)/(x^(j)_i + x⁽⁰⁾_i) where x^(j)_i and x⁽⁰⁾_i are the activity levels of module i after and before the perturbation to module j.⁸¹

In its original formulation, MRA was developed to estimate the local matrix elements that represent the strength of the involved interactions between modules. The local response matrix r_ij is calculated from the measured R_ij by solving , where the sum is over the nodes k (k ≠ i) that were perturbed in the experiments and r_ii = −1.^81,82 Typically the total least squares estimation is used in Monte Carlo based simulations to estimate the coefficients.⁸³ Utilizing this approach, we have computed the estimated probability distributions of the response coefficients by forming random realizations of the node activities that were drawn from a normal distribution with a mean equal to those of the measured values. Standard deviation in the distribution was assumed to be 30% of the mean for each of the measured R_ij values. This deviation was larger than the deviations observed in our Western blot experiments but we preferred to use a conservative standard deviation value in our analysis. Response matrices r_ij were calculated for 10⁷ realizations of the randomly sampled data sets, and statistical distribution data obtained for the resulting local response coefficients were used to calculate the mean (μ_ij = 〈r_ij〉), standard deviation (σ_ij), and coefficient of variation (CV = σ_ij/μ_ij) of the computed response coefficients. These MRA simulations were performed using the in-house written MATLAB codes.

Recently, Santra et al. have developed a variant of MRA, the Bayesian Variable Selection Algorithm (BVSA) that integrates Bayesian selection into MRA to infer regulatory interactions in biochemical networks.⁸⁴ In BVSA, the constitutive MRA equation is changed to , where the added variable A_ij represents if direct interaction between modules i and j are present (A_ij = 1) or absent (A_ij = 0). For noisy datasets, which is the usual case in biological experiments, the global responses R_ij contain uncertainties and the above equality usually does not hold exactly. The BVSA algorithm estimates the probability of module i being directly influenced by module j without having to estimate the connection coefficients as in the MRA algorithm. As a statistical inference algorithm, BVSA uses Bayesian statistics: starting with a prior distribution, distributions are updated based on the experimental data using Bayes' theorem to obtain the posterior distributions for the probabilities of direct connections between network modules. Santra et al. utilized Markov Chain Monte Carlo sampling to obtain the posterior distributions for A_ij, and decided on the inferred network topology by defining a threshold equal to the average posterior edge probability.⁸⁴ We have used the Gibbs sampler MATLAB program provided as ESI† to their publication in our simulations. Simulations were run for 1000 Gibbs scans (parameter: noit) and 5000 iterations (parameter: times). 50% of the early samples were assumed to be burn-ins and were not considered in the calculation of the posterior edge probabilities. Prior and posterior distributions were computed using the default values for the other parameters of the program.

In this study, we have combined these two algorithms in a two-stage approach to first determine which regulatory interaction links are likely to exist and then to determine the strength and type (activation/repression) of interactions. First, we have used BVSA to identify which interaction links may be present in the investigated biological network, i.e., find which A_ij are not zero. In the second MRA stage, we have set the local response matrix r_ij elements corresponding to the interactions that were found to be non-existing in the first round to zero and computed the values of r_ij for the interactions that were found to exist in our network. The second stage allowed us to rank order the strengths of the interactions and to characterize if an interaction is activation or repression relation.

Since we have obtained consistent datasets in four cell lines with same cellular background, we also investigated if the detected interactions depended on the expression of either HER2 or HER3 receptors. This was achieved by grouping the perturbation response data in three different ways. First, we defined an “all dataset” by considering all four cell lines and all three ligand-treatment cases. For this set R_ij data consisted of 12 separate conditions in which each sentinel protein was inhibited as perturbation for each condition and the activity (phosphorylation) changes of all proteins, i.e., x^(j)_i, were measured to determine the R_ij values. We also defined “HER2 and HER3 specific datasets” by comparing the responses in our cell lines. Cell lines in our library (Table 1) can be contrasted to each other to signify a particular HER receptor: for example, the principal difference between the R2−3− (Par) and R2+3− (24H) and between the R2−3+ (B5) and R2+3+ (D20) cell lines is the HER2 expression when HER1/EGFR and HER3 levels are roughly the same. Therefore, to form a “HER2 specific perturbation response dataset” from our measurements, we computed the differences between the measured responses in pairs of cell lines. i.e., the response matrix for HER2 specific data set was obtained by computing the R_ij(R2+3−) − R_ij(R2−3−) and R_ij(R2+3+) − R_ij(R2−3+) difference for every ligand combination. This resulted in a dataset for 6 conditions (2 cellular × 3 ligand treatments) where the main difference was the HER2 expression. “HER3 specific perturbation response dataset” for 6 different conditions was formed in a similar manner by computing R_ij(R2−3+) − R_ij(R2−3−) and R_ij(R2+3+) − R_ij(R2+3−) differences for every ligand combination where the main difference was the HER3 expression.

Results

Correlations of JNK and STAT3 with other proteins

Our cluster analysis indicated that JNK and STAT3 activation patterns only weakly correlated to HER, p38 MAPK, ERK, and AKT activations (Fig. 1A and 2A). Both STAT3 and JNK were more associated with the p38 MAPK/ERK/HER1 cluster than the AKT/HER3 cluster (Fig. 1A and 2A). However, the linkage distances were large indicating that auxiliary JNK and STAT3 pathways have rather weak correlative associations with the p38 MAPK, AKT, and ERK pathways. Correlation between JNK and STAT3 and other sentinels were also weak when the “treatment-type” and “cell-type” based datasets (cf., Computational methods section) were used in the cluster analysis (data not shown). Modular response analysis results discussed below further indicated that the activation patterns of STAT3 and JNK were complex and they may be differentially regulated by both ERK and AKT pathways. We therefore exclude JNK and STAT3 from the detailed discussion of the activation pattern correlations, and discuss the involvement of JNK and STAT3 in the regulations of other proteins later below.

Fig. 1 Cluster analysis of the time-dependent activation patterns of the investigated sentinel proteins and the HER1–3 receptors: (A) analysis using the all data set, i.e., data combined for 4 cell lines and 3 ligand-treatments; (B) ligand-treatment based analysis using the data grouped for 4 cell lines for EGF (E), HRG (H), and both ligands (EH) treatment conditions; and (C) cell-type based analysis using the data grouped for the 3 ligand treatments for each cell line, i.e., for the cells R2−3− (−−), R2+3− (+−), R2−3+ (−+), and R2+3+ (++). The value of y-axis presents the score of 1 − |Corr| where Corr is the correlation between the activities of the entities as discussed in the text.

Fig. 2 Cluster analysis of the time-dependent activation patterns of the investigated sentinel proteins and the HER receptor dimers. (A) Analysis using the all data set, i.e., data combined for 4 cell lines and 3 ligand-treatments; (B) ligand-treatment based analysis using the data grouped for 4 cell lines for EGF (E), HRG (H), and both ligands (EH) treatment conditions; and (C) cell-type based analysis using the data grouped for the 3 ligand treatments for each cell line, i.e., for the cells R2−3− (−−), R2+3− (+−), R2−3+ (−+), and R2+3+ (++). The value of y-axis presents the score of 1 − |Corr| where Corr is the correlation between the activities of the entities as discussed in the text. In the plots, Rarb indicates the activation of HERa receptors in HERa:HERb dimer complexes.

Correlation analysis without consideration of receptor dimers

We first analyzed the data using the receptor phosphorylation information without considering receptor dimers. In this “receptor based analysis”, we address the question of whether ERK, AKT and p38 MAPK activation can be explained with the total HER receptor phosphorylation levels. This is the typical approach used in a large majority of studies in the literature.

Fig. 1 reports the clusters when analysis was pursued using the all, treatment-type, and cell-type based datasets. All data set results (Fig. 1A) showed nicely clustered correlations between (i) AKT and HER3, and (ii) p38 MAPK, ERK and HER1/EGFR. Even though HER2 is believed to be the mediator of HER3 signaling and its tumorigenic potency,^1,10 interestingly, its activation did not correlate well with the other receptors and sentinel proteins. Ability to partition the data to cell-type and treatment-type subsets allowed us to investigate the dependence of these high level correlations on receptor expression and treatment conditions. In the analysis based on ligand-treatment (Fig. 1B), AKT correlated more with HER1, ERK and p38 MAPK instead, while HER2 and HER3 clustered together. Whereas, the correlation between AKT and HER3 still held on a cell-type basis (Fig. 1C) but the tight grouping between p38 MAPK, ERK and HER1 seemed to disintegrate at the cell-type level and HER1 correlated more with HER2.

Even though it may capture and reveal the correlations between EGFR/HER receptor activation and their signaling pathways, these results clearly illustrate that the interpretation of the data can depend on the way the receptor based analysis is performed. Most alarmingly, even though the analyzed data set was internally consistent, different analyses were found to lead to conflicting conclusions.

Correlation analysis using receptor dimers

We have repeated the cluster analysis using HER dimer activation levels instead of the experimentally measured receptor phosphorylation data. Contributions of dimer types to the overall receptor phosphorylations were inferred from the experimental data in a model-centric analysis (cf., Computational methods section).

Fig. 2 reports the clustering results for the receptor dimer-based approach. All data set results (Fig. 2A) paralleled the results in Fig. 1A: AKT clustered with HER3, and p38 MAPK, ERK and HER1 clustered together. Fig. 2A provides further details about which dimers participate in the correlations: AKT was most correlated with the HER3 signaling from the HER3:HER1 (R3r1) dimers. In contrast, p38 MAPK and ERK were correlated the most with HER1 signaling from its homo-dimers. As in the receptor based analysis discussed above, obtained clusters showed that activation patterns of HER2 correlated poorly with the activation patterns of the sentinel proteins for the signaling pathways. Fig. 2B and C report the results when the cluster analysis was pursued using the ligand treatment-type and cell-type data subsets, respectively. Unlike the receptor-based analysis reported in the previous subsection, where different analyses led to conflicting correlations, ligand treatment type (Fig. 2B) and cell-type (Fig. 2C) dependent clustering agreed to a large degree with the results of the analysis with all data set. Most notably: (i) AKT activation correlated the most with HER3 signaling from HER3:HER1 dimers in both the ligand treatment- and cell-type dependent analyses (Fig. 2B and C). Ligand treatment-based analysis indicated additional correlation between AKT and HER1 signaling from HER1:HER3 dimers. Interestingly, AKT activation when the receptors were stimulated with EGF better correlates with p38 MAPK/ERK/HER1 activation (AKT(E) in Fig. 2B; cluster on the left) while HRG stimulation by itself, AKT(H), or together with EGF, AKT(EH), leads to a better correlation with signaling from the HER3:HER1 dimer (Fig. 2B; cluster on the right). This could be a reflection of the known multi-functionality of the AKT pathway where its contribution to cell signaling depends on the received cues.⁸⁵ (ii) In agreement with the all-data set results, analyses using data subsets did not indicate any strong correlation between HER2 signaling and sentinel proteins. (iii) In the cell-line dependent analysis (Fig. 2C), p38 MAPK showed a stronger correlation with HER1 signaling from HER1:HER3 dimer and ERK was better correlated with HER1 signaling from HER1:HER1 homo- and HER1:HER2 hetero-dimers. However, p38 MAPK and ERK clusters merged at a score ∼0.35, i.e., a correlation of 0.65, which was still a significant correlation. Thus, these results were in good agreement with the higher level analysis results (Fig. 2A). A notable difference in the analyses results for the p38 MAPK/ERK/HER1 cluster was that HER2 containing dimers became part of this cluster in cell-line based analysis (Fig. 2C).

Modular response analysis of the inhibition studies

Investigation of the correlations between protein activation levels can reveal the associations between signaling pathways but they do not provide direct information regarding causal effects. We have therefore also pursued a systematic study where we have inhibited one of the proteins using specific inhibitors and measured the activation levels of the other studied proteins. This data (Fig. S1–S6, ESI†) would indicate the impact of the inhibition of an “effector”-protein on its “target”-protein. We have analyzed the obtained data using modular response analysis to infer a model for the interactions between ERK, AKT, PI3K, p38 MAPK, JNK and STAT3 in our HME cells.

Our quantitative systems approach does not focus on a particular response condition but uses the general trends observed across all of the studied conditions to determine robust crosstalk interactions to obtain a consensus interaction model (Fig. 3). A unique aspect of our study is that ability to obtain consistent datasets in four cell lines with same cellular background but different HER2/3 expression and ability to differentially activate particular HER receptor types using different ligand regimens allowed us to infer the interaction model in a receptor type specific manner. This was achieved by grouping the perturbation response data in three different ways, as described in the Computational methods section.

Fig. 3 Consensus regulatory interaction diagram between the studied sentinel proteins as derived from the inhibition experiments described in the text. Arrows with a pointed end indicate activation interaction and with a blunt end indicate repression interaction. Tags next to the arrows indicate if the interaction was observed regardless of HER2 or HER3 expression (all); or particularly in HER2 (R2) or HER3 (R3) expressing cells. Further details are explained in the main text.

The probabilities of the regulatory reaction interactions as determined in the BVSA analysis are tabulated in Table S1 (ESI†). To decide on the regulatory interactions, Santra et al. used the mean of the computed response matrix element values as the threshold. Mean value of the three datasets μ_r was 0.34 in our case with a standard deviation σ_r of 0.27, which is quite lower than a mean of 0.5 that would be for random interactions. We have therefore classified the response matrix elements into 3 classes (Table S1, ESI†): interactions with r values that are (i) r < μ_r, (ii) μ_r < r < μ_r + σ_r, and (iii) r > μ_r + σ_r. Interactions in the first category are unlikely to be present and the ones in the last category are most likely to be present in our system. Interactions with r probabilities in the intermediate second category are probable but they are hard to estimate with our data and their prediction is uncertain.

We have then pursued MRA in the second stage of our analysis (Table S2, ESI†). We have first computed the response matrix element values by including all of the possible links in the network, i.e., all 5 × 4 = 20 binary possibilities were included. Second, MRA was pursued by including the interactions detected in the BVSA (Table S2, ESI†). Further details of this analysis are described in the Methods section above.

The consensus model shown in Fig. 3 was obtained by combining the BVSA and MRA results, and it revealed that STAT3 occupies a central role in the crosstalk between studied proteins. Some of the interactions were consistently observed across different analysis. Most interestingly, as we have expected, our analysis has established that some of the regulatory interactions were more apparent in HER2 or HER3 expressing cell lines. We next discuss the combined analysis results for the identified interactions (Fig. 3).

Interaction between AKT and STAT3. BVSA analysis (Table S1, ESI†) indicated that AKT is regulated by STAT3 and this regulation was enhanced with HER3 expression. MRA computations (Table S2, ESI†) supported this result and further revealed that AKT is activated by STAT3 and this activation was more pertinent in HER3 expressing cells. AKT in turn was found to activate STAT3 in all cases strongly (Tables S1 and S2, ESI†).

Interactions between ERK, STAT3 and JNK are bi-directional and HER2 dependent. BVSA analysis has indicated that ERK, STAT3 and JNK interact with each other in a bi-directional manner (Table S1, ESI†). MRA computations (Table S2, ESI†) supported this result and further revealed that ERK–JNK and STAT3–JNK are activating interactions while ERK and STAT3 repress each other. All of these interactions were predicted to be more dominant when HER2 is expressed with the exceptions that STAT3 was found to be activated by JNK in all cases and repressed by ERK when either HER2 or HER3 is expressed. Based on the derived interaction strengths, mutually repressive interaction between STAT3 and ERK is the most questionable of all the binary interactions between sentinels. This was reflected as large standard deviations in the interaction values derived for the regulation of ERK (Table S2, ESI†).

Interactions between p38 MAPK, STAT3 and JNK are triangular and HER receptor expression dependent. Similar to the ERK–STAT3–JNK case, p38 MAPK was found to interact with STAT3 and JNK in a triangular manner. BVSA analysis indicated that these interactions may be HER expression dependent (Table S1, ESI†) and its findings were supported by the MRA analysis (Table S2, ESI†). STAT3–JNK interaction has been discussed in the previous paragraph. STAT3 and p38 MAPK were found to mutually repress each other in all cells (Tables S1 and S2, ESI†). JNK and p38 MAPK interaction was mutually activating in HER3 expressing cells (Tables S1 and S2, ESI†).

Comparison of the consensus model with cluster analysis results

The obtained consensus model can be compared and contrasted to the cluster analysis results to retrospectively understand the sources of the observed signaling correlations:

ERK and p38 MAPK correlation is indirect. Cluster analysis of the phosphorylation patterns discussed above indicated that ERK and p38 MAPK activations were closely correlated. Interestingly, modular response analysis predicted that there was no apparent direct impact of inhibiting ERK or p38 MAPK on the activation level of the companion protein. Derived model however illustrated that p38 MAPK and ERK activities are coordinated through STAT3 and JNK. Cluster analysis indicated that p38 MAPK/ERK/HER1 activation correlated better with the HER2 containing dimers (Fig. 2C). Obtained consensus model has validated this finding by showing that the triangulated interactions between ERK/JNK/STAT3 and JNK/STAT3/p38 MAPK are HER2 expression dependent. According to the model, ERK and p38 MAPK activity may be correlated through STAT3 (ERK is high, STAT3 is low, p38 MAPK is high; or ERK is low, STAT3 is high, p38 MAPK is low) or through JNK.

STAT3 is activated by AKT and suppressed by ERK. Our analysis has shown that AKT is an activator of STAT3 and ERK is a repressor of STAT3 in either HER2 or HER3 expressing cells (Fig. 3, Table S2, ESI†). These results indicate contrasting roles for ERK and AKT in the regulation of STAT3 in HME cells when responding to stimulation with different growth factors. Therefore, STAT3's correlative association with either ERK or AKT pathways is only partial, which gets reflected as low correlations in the cluster analysis.

JNK is mainly stimulated through ERK pathway and it is HER2 dependent. Obtained model indicated an interaction between ERK and JNK, particularly in the HER2 expressing cells (Fig. 3). In contrast, AKT stimulation was found to have low impact on JNK activation (Table S1, ESI†). These results supported the correlation analysis results (Fig. 1 and 2) that signaling through JNK pathway is more associated with the ERK pathway and likely to involve the EGFR/HER1 and HER2 receptors and not the HER3/AKT pathway.

STAT3 and JNK regulate p38 MAPK in a context dependent manner. We have found that STAT3 repressed p38 MAPK activation (Fig. 3). JNK in general tends to contribute to this suppression by activating STAT3. However, in HER3 expressing cells JNK relieves p38 MAPK's suppression by activating it directly. Thus, the p38 MAPK/JNK/STAT3 interactions are complex and context dependent in HME cells. When cell signaling is initiated through HER3/AKT pathway, activity of p38 MAPK would be repressed by STAT3, and enhanced through JNK. For the signaling via HER1/HER2/ERK pathway, p38 MAPK is suppressed by the concerted effort of STAT3 and JNK.

AKT activation. We have also investigated how AKT phosphorylation was affected by PI3K inhibition, and found that inhibition of PI3K almost diminishes the AKT activation (Fig. S6, ESI†). Activation of AKT by PI3K is well established in the literature,^86–88 and our results confirmed it for our cell lines. Our results also established that AKT and ERK activations were not impacted by the inhibition of each other (Fig. S1B and S2A, ESI†) and there were not any detected interactions (Tables S1 and S2, ESI†). These results establish that the crosstalk between the two signaling pathways³³ is not direct in our HME cells.

Discussion

Having access to a closely related cell library in which multiple signaling pathways could be activated differentially via different HER receptors allowed us to pursue a systematic study of crosstalk between proliferation, pro-survival and stress response pathways in HME cells under growth factor stimulated signaling conditions. Our first aim was to establish the possible relationships between the HER family member receptors and several key proteins that are believed to be sentinels for the investigated signaling pathways. Towards this aim, we have collected the experimental data for four HME cell lines and for three ligand treatment conditions. Analysis of the collected data using correlation based cluster analysis has allowed us to establish that: (i) AKT activation correlates the most with HER3 phosphorylation, in particular with the HER3 signaling from the HER3:HER1 dimer. (ii) ERK and p38 MAPK activations correlate the most with HER1 phosphorylation, in particular HER1 signaling from its homo-dimers, and with the HER1:HER2 hetero-dimer signaling. (iii) Even though HER2 is believed to be the mediator of HER3 signaling and its increased potency,^3,89–95 the HER2:HER3 dimers correlated only weakly with the activation of signaling proteins. (iv) Activations of JNK and STAT3 did not correlate strongly with the activation patterns of the other studied proteins in HME cells under the investigated conditions. STAT3 and JNK were more associated with p38 MAPK/ERK/HER1, albeit weakly, than with AKT/HER3. Modular response analysis results implied that the observed weak correlation was most likely due to the complex crosslink relationship between AKT and ERK pathways through STAT3.

Most importantly, our analysis led to the conclusion that a consistent analysis of the HER signaling pathways has to take dimer-based receptor signaling into account. There is considerable evidence that the HER dimerization pattern is an important determinant of the consequences of HER signaling,^62,63 which suggests that activation of the downstream signaling pathways would be functions of receptor dimer formation and activation.¹⁸ In this regard, we have hypothesized that the HER dimerization pattern is a better predictor of downstream signaling than total receptor activation levels. Correlation analysis using subsets of the experimental data led to cohesive interpretations only when the protein activation patterns were investigated in relationship to receptor signaling from dimer complexes. Analysis based on receptor measurements that did not consider dimer-dependent signaling resulted in conflicting interpretations. This finding can also be seen as a cause of alarm for the interpretation of experimental results when the dynamics of cells are stimulated through HER. Studies in the literature often investigate the activation profiles of the signaling proteins for a particular condition that is relevant for the addressed problem. Our results however clearly illustrated that the interpretation can depend on the way the analysis is performed.

Our correlation analysis did not take into account that there may be some time delay between the receptor activation and activation of the downstream elements. As we have shown earlier,¹⁸ the delay between the receptor and ERK/AKT activations is on the order of 3 minutes in the HME cells. So use of equal-time correlation in the reported analysis was reasonable.

To extend our analysis to identify possible causal regulatory relationships, we have collected the experimental data for how inhibition of a protein affects the activation levels of the other investigated proteins. Our ability to generate such data for multiple cell lines and ligand stimulation conditions made it possible to derive causal relationships that are consensus based instead of being applicable to particular conditions. This analysis led us to establish the regulatory interaction model between the studied proteins ERK, AKT, p38 MAPK, STAT3, and JNK in our HME cells.

Derived consensus model, which is shown in Fig. 3 and quantitatively described in Tables S1 and S2 (ESI†), illustrates that the proliferation and pro-survival cell signaling pathways are tightly cross-linked in HME cells. MRA of the collected inhibition data led us to identify a central switch-like regulatory role for STAT3 in HME cells with contrasting mechanisms for ERK and AKT in the regulation of STAT3 when responding to stimulation with different growth factors where STAT3 is suppressed by the ERK pathway and activated by the AKT pathway. MRA also identified triangular regulatory patterns among sentinels, which can be considered as evidence of the possible multi-functionality of the studied pathways and how central signaling pathways may be utilized by cells for differentially regulating their responses as a function of their physiological state.

Our analysis indicated HER2 or HER3 receptor driven interactions as well. Triangular interactions between ERK, STAT3 and JNK are mainly operational in HER2 expressing cells while the interactions between p38 MAPK, STAT3 and JNK have a more complex HER expression dependence and can be context dependent in HME cells. For the signaling via HER1/HER2/ERK pathway, p38 MAPK is suppressed by both STAT3 and JNK. When cell signaling is initiated through HER3/AKT pathway, activity of p38 MAPK would be repressed by STAT3 and enhanced through JNK. As JNK is regulated by ERK in HER2 expressing cells, ERK also contributes to STAT3's regulation indirectly through JNK. Thus, STAT3 can couple the ERK and AKT pathways in either direction and its regulatory role depends on pathway stimulation through activation of HER2. Overall, our analysis has established that the crosstalk between the canonical ERK and AKT pathways can occur at multiple levels and it can be dependent on the involvement of particular HER receptors.

As seen in the consensus model (Fig. 3), almost all of the identified interactions are bi-directional where the interactions are either both activating or repressing interactions in both directions. A bi-directional interaction that is activating in both directions is a positive feedback loop where an activating stimulation leads to the high activity of both proteins (i.e., both on or both off). A bi-directional interaction that is repressing in both directions is a feedback loop that leads to bi-stable dynamics where the activity of one of the proteins is high and of the companion protein is low (i.e., one is on and one is off as in a genetic toggle switch⁹⁶). Dynamic response of bi-stable circuits depends on the initial conditions of the system when the kinetic rates, i.e., interaction strength, are comparable in forward and reverse directions. Protein with higher activity level initially would continue to stay active while the companion protein's activity would be suppressed until a signal flips the states of the proteins. With these properties in mind, obtained likelihood of interaction to exist (Table S1, ESI†) and interaction strength and type (Table S2, ESI†) data may also be used to interpret the hierarchy of regulatory order in HME cells. For example, our analysis results indicated that STAT3 regulation by AKT to be more dominant than its regulation by p38 MAPK or JNK (Table S2, ESI†).

The strong correlation between ERK and p38 MAPK observed in the cluster analysis, and the lack of direct interaction between them as found in modular response analysis may appear as contradictory. However, based on the model shown diagrammatically in Fig. 3, ERK activation represses STAT3 and suppressed STAT3 activities lift the repression on p38 MAPK. Thus, ERK activation can lead to p38 MAPK activation through STAT3 in a positively correlated manner. Such a positive correlation is also possible for the activation of p38 MAPK via JNK upon signaling through ERK pathway. Therefore, the results of the cluster analysis and the modular response analysis are not inconsistent with each other.

Involvement of p38 MAPK, JNK, and STAT3 in the interplay between ERK and AKT signaling has been studied in the literature, where the conclusions seem to depend on the cell type and treatment conditions. To avoid such complications and possible contradictions, our study was designed to collect and analyze data that is internally consistent and abundant enough to allow more definitive conclusions. Model that we have constructed using the results of our analysis indicated that STAT3 is the focal point of the crosstalk between ERK and AKT pathways in the HME cells. STAT3 is aided by JNK and p38 MAPK in modulating the cellular responses to stimulation with growth factors. Our results indicated a less prominent role for JNK in coupling p38 MAPK activation to ERK. Overall structure of the signaling network connecting ERK and AKT to downstream elements consists of several positive and negative feedback loops, thus illustrating the complexity of the cellular signaling networks and that pathways do not function in isolation and wider profiling of pathways is needed to decipher the cellular responses reliably. We are in the process of extending these studies to (in)validate the derived network in other cell lines and to include the role of other related sentinel proteins.

Acknowledgements

The research described in this paper was funded by the National Institutes of Health Grant 5R01GM072821-07 to H.R. Pacific Northwest National Laboratory is a multiprogram national laboratory operated by Battelle for the U.S. Department of Energy under Contract DE-AC06-76RL01830.

References

1. Y. Yarden and M. X. Sliwkowski, Nat. Rev. Mol. Cell Biol., 2001, 2, 127–137.
2. N. Normanno, A. De Luca, C. Bianco, L. Strizzi, M. Mancino, M. R. Maiello, A. Carotenuto, G. De Feo, F. Caponigro and D. S. Salomon, Gene, 2006, 366, 2–16.
3. T. Holbro, G. Civenni and N. E. Hynes, Exp. Cell Res., 2003, 284, 99–110.
4. R. C. Harris, E. Chung and R. J. Coffey, Exp. Cell Res., 2003, 284, 2–13.
5. P. M. Guy, J. V. Platko, L. C. Cantley, R. A. Cerione and K. L. Carraway, 3rd, Proc. Natl. Acad. Sci. U. S. A., 1994, 91, 8132–8136.
6. J. Schlessinger, Cell, 2002, 110, 669–672.
7. J. P. Dawson, M. B. Berger, C. C. Lin, J. Schlessinger, M. A. Lemmon and K. M. Ferguson, Mol. Cell. Biol., 2005, 25, 7734–7742.
8. M. X. Sliwkowski, G. Schaefer, R. W. Akita, J. A. Lofgren, V. D. Fitzpatrick, A. Nuijens, B. M. Fendly, R. A. Cerione, R. L. Vandlen and K. L. Carraway, 3rd, J. Biol. Chem., 1994, 269, 14661–14665.
9. R. R. Beerli and N. E. Hynes, J. Biol. Chem., 1996, 271, 6071–6076.
10. D. Graus-Porta, R. R. Beerli, J. M. Daly and N. E. Hynes, EMBO J., 1997, 16, 1647–1655.
11. E. Tzahar, H. Waterman, X. Chen, G. Levkowitz, D. Karunagaran, S. Lavi, B. J. Ratzkin and Y. Yarden, Mol. Cell. Biol., 1996, 16, 5276–5287.
12. J. P. Duneau, A. P. Vegh and J. N. Sturgis, Biochemistry, 2007, 46, 2010–2019.
13. H. Shankaran, Y. Zhang, Y. B. Tan and H. Resat, PLoS Comput. Biol., 2013, 9, e1003201,  DOI:10.1371/journal.pcbi.1003201.
14. W. X. Schulze, L. Deng and M. Mann, Mol. Syst. Biol., 2005, 1, 2005.0008.
15. R. N. Jorissen, F. Walker, N. Pouliot, T. P. Garrett, C. W. Ward and A. W. Burgess, Exp. Cell Res., 2003, 284, 31–53.
16. Y. Yarden and B. Z. Shilo, Cell, 2007, 131, 1018.
17. H. Shankaran, Y. Zhang, L. Opresko and H. Resat, Biochem. Biophys. Res. Commun., 2008, 371, 220–224.
18. Y. Zhang, H. Shankaran, L. Opresko and H. Resat, IET Syst. Biol., 2008, 2, 273–284.
19. M. R. Birtwistle, M. Hatakeyama, N. Yumoto, B. A. Ogunnaike, J. B. Hoek and B. N. Kholodenko, Mol. Syst. Biol., 2007, 3, 144.
20. W. W. Chen, B. Schoeberl, P. J. Jasper, M. Niepel, U. B. Nielsen, D. A. Lauffenburger and P. K. Sorger, Mol. Syst. Biol., 2009, 5, 239.
21. C. D. Britten, Mol. Cancer Ther., 2004, 3, 1335–1342.
22. K. Sakai, H. Yokote, K. Murakami-Murofushi, T. Tamura, N. Saijo and K. Nishio, Cancer Sci., 2007, 98, 1498–1503.
23. B. Tanner, D. Hasenclever, K. Stern, W. Schormann, M. Bezler, M. Hermes, M. Brulport, A. Bauer, I. B. Schiffer, S. Gebhard, M. Schmidt, E. Steiner, J. Sehouli, J. Edelmann, J. Lauter, R. Lessig, K. Krishnamurthi, A. Ullrich and J. G. Hengstler, J. Clin. Oncol., 2006, 24, 4317–4323.
24. C. J. Witton, J. R. Reeves, J. J. Going, T. G. Cooke and J. M. S. Bartlett, J. Pathol., 2003, 200, 290–297.
25. J. A. Engelman, K. Zejnullahu, T. Mitsudomi, Y. Song, C. Hyland, J. O. Park, N. Lindeman, C. M. Gale, X. Zhao, J. Christensen, T. Kosaka, A. J. Holmes, A. M. Rogers, F. Cappuzzo, T. Mok, C. Lee, B. E. Johnson, L. C. Cantley and P. A. Janne, Science, 2007, 316, 1039–1043.
26. N. V. Sergina, M. Rausch, D. Wang, J. Blair, B. Hann, K. M. Shokat and M. M. Moasser, Nature, 2007, 445, 437–441.
27. E. K. Rowinsky, Annu. Rev. Med., 2004, 55, 433–457.
28. C. L. Arteaga, Exp. Cell Res., 2003, 284, 122–130.
29. N. Normanno, M. Campiglio, L. A. De, G. Somenzi, M. Maiello, F. Ciardiello, L. Gianni, D. S. Salomon and S. Menard, Ann. Oncol., 2002, 13, 65–72.
30. P. K. Kreeger and D. A. Lauffenburger, Carcinogenesis, 2010, 31, 2–8.
31. N. Kumar, R. Afeyan, H. D. Kim and D. A. Lauffenburger, Mol. Pharmacol., 2008, 73, 1668–1678.
32. M. Katz, I. Amit and Y. Yarden, Biochim. Biophys. Acta, Mol. Cell Res., 2007, 1773, 1161–1176.
33. E. Aksamitiene, A. Kiyatkin and B. N. Kholodenko, Biochem. Soc. Trans., 2012, 40, 139–146.
34. G. Sheng, J. Guo and B. W. Warner, Am. J. Physiol.: Gastrointest. Liver Physiol., 2007, 293, G599–G606.
35. I. D. Barrantes and A. R. Nebreda, Biochem. Soc. Trans., 2012, 40, 79–84.
36. T. Pawson and N. Warner, Oncogene, 2007, 26, 1268–1275.
37. F. Chiacchiera, V. Grossi, M. Cappellari, A. Peserico, M. Simonatto, A. Germani, S. Russo, M. P. Moyer, N. Resta, S. Murzilli and C. Simone, Cancer Lett., 2012, 324, 98–108.
38. M. R. Frey, R. S. Dise, K. L. Edelblum and D. B. Polk, EMBO J., 2006, 25, 5683–5692.
39. L. J. Hui, L. Bakiri, E. Stepniak and E. F. Wagner, Cell Cycle, 2007, 6, 2429–2433.
40. J. M. Nelson and D. W. Fry, J. Biol. Chem., 2001, 276, 14842–14847.
41. R. M. Neve, T. Holbro and N. E. Hynes, Oncogene, 2002, 21, 4567–4576.
42. E. F. Wagner and A. R. Nebreda, Nat. Rev. Cancer, 2009, 9, 537–549.
43. M. Raman, W. Chen and M. H. Cobb, Oncogene, 2007, 26, 3100–3112.
44. C. J. Marshall, Cell, 1995, 80, 179–185.
45. E. J. Joslin, L. K. Opresko, A. Wells, H. S. Wiley and D. A. Lauffenburger, J. Cell Sci., 2007, 120, 3688–3699.
46. C. J. Wu, X. Qian and D. M. O'Rourke, DNA Cell Biol., 1999, 18, 731–741.
47. G. Song, G. Ouyang and S. Bao, J. Cell. Mol. Med., 2005, 9, 59–71.
48. Y. Zhang, S. Banerjee, Z. Wang, H. Xu, L. Zhang, R. Mohammad, A. Aboukameel, N. V. Adsay, M. Che, J. L. Abbruzzese, A. P. Majumdar and F. H. Sarkar, Cancer Res., 2006, 66, 1025–1032.
49. J. Chen, P. R. Somanath, O. Razorenova, W. S. Chen, N. Hay, P. Bornstein and T. V. Byzova, Nat. Med., 2005, 11, 1188–1196.
50. J. A. Engelman, J. Luo and L. C. Cantley, Nat. Rev. Genet., 2006, 7, 606–619.
51. A. C. Hsieh and M. M. Moasser, Br. J. Cancer, 2007, 97, 453–457.
52. B. Schoeberl, A. C. Faber, D. Li, M. C. Liang, K. Crosby, M. Onsum, O. Burenkova, E. Pace, Z. Walton, L. Nie, A. Fulgham, Y. Song, U. B. Nielsen, J. A. Engelman and K. K. Wong, Cancer Res., 2010, 70, 2485–2494.
53. H. H. Kim, S. L. Sierke and J. G. Koland, J. Biol. Chem., 1994, 269, 24747–24755.
54. S. P. Soltoff, K. L. Carraway, 3rd, S. A. Prigent, W. G. Gullick and L. C. Cantley, Mol. Cell. Biol., 1994, 14, 3550–3558.
55. B. Schoeberl, E. A. Pace, J. B. Fitzgerald, B. D. Harms, L. Xu, L. Nie, B. Linggi, A. Kalra, V. Paragas, R. Bukhalid, V. Grantcharova, N. Kohli, K. A. West, M. Leszczyniecka, M. J. Feldhaus, A. J. Kudla and U. B. Nielsen, Sci. Signaling, 2009, 2, ra31.
56. M. H. Cobb, Prog. Biophys. Mol. Biol., 1999, 71, 479–500.
57. M. R. Junttila, S. P. Li and J. Westermarck, FASEB J., 2008, 22, 954–965.
58. K. L. Mueller, K. Powell, J. M. Madden, S. T. Eblen and J. L. Boerner, Transl. Oncol., 2012, 5, 327–334.
59. J. Dong, S. Ramachandiran, K. Tikoo, Z. Jia, S. S. Lau and T. J. Monks, Am. J. Physiol.: Renal Physiol., 2004, 287, F1049–F1058.
60. J. S. Rawlings, K. M. Rosler and D. A. Harrison, J. Cell Sci., 2004, 117, 1281–1283.
61. M. Sen, S. Joyce, M. Panahandeh, C. Y. Li, S. M. Thomas, J. Maxwell, L. Wang, W. E. Gooding, D. E. Johnson and J. R. Grandis, Clin. Cancer Res., 2012, 18, 4986–4996.
62. M. A. Olayioye, D. Graus-Porta, R. R. Beerli, J. Rohrer, B. Gay and N. E. Hynes, Mol. Cell. Biol., 1998, 18, 5042–5051.
63. L. Yen, N. Benlimame, Z. R. Nie, D. Xiao, T. Wang, A. E. Al Moustafa, H. Esumi, J. Milanini, N. E. Hynes, G. Pages and M. A. Alaoui-Jamali, Mol. Biol. Cell, 2002, 13, 4029–4044.
64. B. S. Hendriks, L. K. Opresko, H. S. Wiley and D. Lauffenburger, Cancer Res., 2003, 63, 1130–1137.
65. Y. Zhang, L. Opresko, H. Shankaran, W. B. Chrisler, H. S. Wiley and H. Resat, BMC Cell Biol., 2009, 10, 78.
66. V. Band and R. Sager, Proc. Natl. Acad. Sci. U. S. A., 1989, 86, 1249–1253.
67. J. Dong, L. K. Opresko, P. J. Dempsey, D. A. Lauffenburger, R. J. Coffey and H. S. Wiley, Proc. Natl. Acad. Sci. U. S. A., 1999, 96, 6235–6240.
68. M. R. Stampfer, C. H. Pan, J. Hosoda, J. Bartholomew, J. Mendelsohn and P. Yaswen, Exp. Cell Res., 1993, 208, 175–188.
69. S. A. Eming, J. Lee, R. G. Snow, R. G. Tompkins, M. L. Yarmush and J. R. Morgan, J. Invest. Dermatol., 1995, 105, 756–763.
70. P. P. Di Fiore, J. H. Pierce, M. H. Kraus, O. Segatto, C. R. King and S. A. Aaronson, Science, 1987, 237, 178–182.
71. O. Danos and R. C. Mulligan, Proc. Natl. Acad. Sci. U. S. A., 1988, 85, 6460–6464.
72. K. J. Garton, N. Ferri and E. W. Raines, Biotechniques, 2002, 32, 830.
73. H. Shankaran, Y. Zhang, W. B. Chrisler, J. A. Ewald, H. S. Wiley and H. Resat, Mol. BioSyst., 2012, 8, 2868–2882.
74. A. Citri, K. B. Skaria and Y. Yarden, Exp. Cell Res., 2003, 284, 54–65.
75. M. B. Berger, J. M. Mendrola and M. A. Lemmon, FEBS Lett., 2004, 569, 332–336.
76. N. Jura, Y. B. Shan, X. X. Cao, D. E. Shaw and J. Kuriyan, Proc. Natl. Acad. Sci. U. S. A., 2009, 106, 21608–21613.
77. E. M. Bublil, G. Pines, G. Patel, G. Fruhwirth, T. Ng and Y. Yarden, FASEB J., 2010, 24, 4744–4755.
78. S. E. Telesco, A. J. Shih, F. Jia and R. Radhakrishnan, Mol. BioSyst., 2011, 7, 2066–2080.
79. F. M. Shi, S. E. Telesco, Y. T. Liu, R. Radhakrishnan and M. A. Lemmon, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 7692–7697.
80. F. J. Bruggeman, H. V. Westerhoff, J. B. Hoek and B. N. Kholodenko, J. Theor. Biol., 2002, 218, 507–520.
81. B. N. Kholodenko, A. Kiyatkin, F. J. Bruggeman, E. Sontag, H. V. Westerhoff and J. B. Hoek, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 15245.
82. M. Andrec, B. N. Kholodenko, R. M. Levy and E. Sontag, J. Theor. Biol., 2005, 232, 427–441.
83. S. D. M. Santos, P. J. Verveer and P. I. H. Bastiaens, Nat. Cell Biol., 2007, 9, 324–U139.
84. T. Santra, W. Kolch and B. N. Kholodenko, BMC Syst. Biol., 2013, 7, 57.
85. N. Grabinski, K. Bartkowiak, K. Grupp, B. Brandt, K. Pantel and M. Jucker, Cell. Signalling, 2011, 23, 1952–1960.
86. K. D. Rodland, N. Bollinger, D. Ippolito, L. K. Opresko, R. J. Coffey, R. Zangar and H. S. Wiley, J. Biol. Chem., 2008, 283, 31477–31487.
87. H. Li, G. Schmid-Bindert, D. Wang, Y. Zhao, X. Yang, B. Su and C. Zhou, Adv. Med. Sci., 2011, 56, 275–284.
88. G. W. Krystal, G. Sulanke and J. Litz, Mol. Cancer Ther., 2002, 1, 913–922.
89. T. Holbro, R. R. Beerli, F. Maurer, M. Koziczak, C. F. Barbas, 3rd and N. E. Hynes, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 8933–8938.
90. K. Zhang, J. Sun, N. Liu, D. Wen, D. Chang, A. Thomason and S. K. Yoshinaga, J. Biol. Chem., 1996, 271, 3884–3890.
91. K. L. Carraway, 3rd, S. P. Soltoff, A. J. Diamonti and L. C. Cantley, J. Biol. Chem., 1995, 270, 7111–7116.
92. R. Pinkas-Kramarski, L. Soussan, H. Waterman, G. Levkowitz, I. Alroy, L. Klapper, S. Lavi, R. Seger, B. J. Ratzkin, M. Sela and Y. Yarden, EMBO J., 1996, 15, 2452–2467.
93. M. Alimandi, A. Romano, M. C. Curia, R. Muraro, P. Fedi, S. A. Aaronson, P. P. Di Fiore and M. H. Kraus, Oncogene, 1995, 10, 1813–1821.
94. C. A. Maurer, H. Friess, B. Kretschmann, A. Zimmermann, A. Stauffer, H. U. Baer, M. Korc and M. W. Buchler, Hum. Pathol., 1998, 29, 771–777.
95. P. M. Siegel, E. D. Ryan, R. D. Cardiff and W. J. Muller, EMBO J., 1999, 18, 2149–2164.
96. T. S. Gardner, C. R. Cantor and J. J. Collins, Nature, 2000, 403, 339–342.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4mb00471j

‡ Current address: Department of Medicine, University of Tennessee Health Science Center, Memphis, TN, 38103, USA.

§ Current address: Modeling and Simulation Group, AstraZeneca Pharmaceuticals, Waltham, MA, 02451, USA.

¶ School of Chemical Engineering and Bioengineering, Washington State University, P.O. Box 646515, Pullman, WA 99164-6515, USA.

---