Chunhong
Gong
a,
Yi
Zhang‡
a,
Harish
Shankaran§
a and
Haluk
Resat¶
*abc
aComputational Biology and Bioinformatics Group, Pacific Northwest National Laboratory, Richland, WA 99352, USA
bSchool of Chemical Engineering and Bioengineering, Washington State University, Pullman, WA 99164, USA. E-mail: haluk.resat@wsu.edu; Tel: +1-509-335-6579
cSchool of Electrical Engineering and Computer Science, Washington State University, Pullman, WA 99164, USA
First published on 2nd October 2014
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.
Biomedical and biomolecular studies have shown that inhibition of HER dimerization and HER-mediated signaling can be an extremely effective therapeutic strategy.21 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.14 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 migration45 and enhances cell transformation and resistance to apoptosis.46 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.47 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 growth50 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.56 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.39 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.59 STAT signaling promotes RAS/RAF/ERK signaling indirectly.60 Activation of MAPKs leads to the activation of STAT-family transcription factors, which promotes their translocation to the nucleus32,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,13 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.
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.13 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,13 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.
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.72 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−1; Invitrogen, Carlsbad, CA) and/or puromycin (2 μg ml−1; 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 studies13,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.
Cell linesa | 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 |
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.
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.79 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 30000 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.13 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.
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.
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 rij is calculated from the measured Rij by solving , where the sum is over the nodes k (k ≠ i) that were perturbed in the experiments and rii = −1.81,82 Typically the total least squares estimation is used in Monte Carlo based simulations to estimate the coefficients.83 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 Rij 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 rij were calculated for 107 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 = 〈rij〉), 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.84 In BVSA, the constitutive MRA equation is changed to , where the added variable Aij represents if direct interaction between modules i and j are present (Aij = 1) or absent (Aij = 0). For noisy datasets, which is the usual case in biological experiments, the global responses Rij 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 Aij, and decided on the inferred network topology by defining a threshold equal to the average posterior edge probability.84 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 Aij are not zero. In the second MRA stage, we have set the local response matrix rij elements corresponding to the interactions that were found to be non-existing in the first round to zero and computed the values of rij 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 Rij 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 Rij 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 Rij(R2+3−) − Rij(R2−3−) and Rij(R2+3+) − Rij(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 Rij(R2−3+) − Rij(R2−3−) and Rij(R2+3+) − Rij(R2+3−) differences for every ligand combination where the main difference was the HER3 expression.
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.
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.85 (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).
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.
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).
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.18 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,18 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 switch96). 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.
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. |
This journal is © The Royal Society of Chemistry 2015 |