Measuring phosphorylation and castration resistant survival in response to treatment
To obtain a diverse response across multiple phosphosites in LNCaP, PC3, and MDA-PCa-2b cells, the cells were treated with the ligands EGF (Epidermal growth factor), IGF1 (insulin-like growth factor 1), IL6, TNFα (Tumor necrosis factor alpha), dihydrotestosterone (DHT) which is an androgen receptor agonist, and the chemotherapeutic docetaxel. LNCaP cells were additionally treated with the targeted kinase inhibitors LY294002 (Phosphoinositide 3-kinase (PI3K) inhibitor), U0126 (Mitogen-activated protein kinase kinase (MEK) inhibitor), wedelactone (IκB Kinase (IKK) α/β inhibitor), temsirolimus (Mammalian target of rapamycin (mTOR) inhibitor), and SB202190 (p38 inhibitor), each in combination with the previously mentioned ligands (Additional file 1: Table S4). These ligands and drugs were selected because of their involvement in moderating prostate signaling pathways which have been implicated in castration resistant growth of prostate cancer, as well as their availability and characterized activity. Whole cell lysates were collected at 30 minutes, 4 hours, and 24 hours post treatment and assayed using 384 well plate phospho-ELISA assays to measure the response of phosphorylation sites in key pathways to treatment with these ligands and inhibitors. In the signaling pathways diagram, a simplistic representation of the interactions between the measured phosphoproteins, the pathways which contain those proteins, and the effect of the targeted inhibitors can be observed (Figure 1A). The phosphosites which were measured in response to treatment are listed (see Phosphosites measured). These particular phosphosites were selected based on an examination of the literature, and their potential to enable cell growth in androgen depleted conditions.
Phosphosites measured
Erk1: T202/Y204 [15]
Erk2: T185/Y187
Akt1: S473 [15, 16]
Akt2: S474
Akt3: S472
RPS6: S235/S236 [16]
GSK3α: S21 [17]
GSk3β: S9
p38δ: T180/Y182 [18]
JNK1 and JNK2: T183/Y185 [19]
JNK3: T221/Y223
HSP27: S78/S82 [20]
Stat3: Y705 [21]
After the phosphoprotein data was collected and normalized (see methods), hierarchical clustering analysis was applied across the phosphosites at the three time points as well as the treatment groups. This analysis measures the similarity between each observation using a Euclidean distance metric (Figure 1B). Across the y dimension of the X matrix, the treatments were found to cluster first by cell line and then by inhibitor treatment (for LNCaP cells only), with little clustering in the ligand treatment groups (Figure 1B and Additional file 2: Table S1). In the x dimension the phosphoprotein activation was generally found to cluster the three time points of each phosphoprotein together (Figure 1B and Additional file 3: Table S2). This clustering indicated that the cell line, and then inhibitor, and finally the ligand treatment imparted the most substantial changes in the cells in the y dimension (treatment conditions). In the × dimension (phosphoproteins), the data indicated that the change by time point tended to cause the most substantial response in phosphoprotein levels.
For each treatment, biological duplicates were measured and the absolute percentage difference between the two replicates was determined (Additional file 4: Table S3). A mean difference of 20.4% was observed across all cell lines which when compared to the finding that the phosphosites varied by approximately 670% on average over untreated controls, was considered an acceptable amount of error.
Regression analysis correlating phosphoprotein measurements to cell survival in androgen depleted conditions
In an attempt to understand how the alterations in signaling may lead to variations in survival outcomes in cells grown in androgen depleted conditions, we built a statistical model using PLS regression. The data was arranged so that the phosphoprotein data was regressed against the survival data using PLS regression on the complete data set of 8 phosphoproteins, at 3 time points, using 3 cell lines, with 6 treatments. After calculating the model parameters the leave-one-out cross validated R2 value was determined to be 0.616 with 3 latent variables, and the predicted versus measured survival values were plotted (Figure 2A). Additional latent variables beyond three had marginal explanatory power due to the fact that the majority of the variation in the X matrix could be described in terms of these latent variables, therefore three components were used for all analyses. When this calculated R2 value was compared to the mean R2 value calculated from randomized models (X matrix rows randomized against Y matrix rows) we observed that this model was 6.36 standard deviations above the mean randomized value of 0.1847 corresponding to a P-value less than 0.0001. This result indicates that this model can correlate to survival significantly better than by random chance.
Upon determining that this model was significantly more accurate than a randomized model, we examined the regression coefficients to determine weights calculated on the different phosphoproteins. Consistently positive coefficients for p-Erk (Extracellular signal-regulated kinases) were noted, as well as consistently increased p-RPS6 (Ribosomal protein S6) across all time points (Figure 2B). p-JNK regression coefficients were negative at all time points along with p-Akt and p-Stat3. p-GSK3 (Glycogen synthase kinase 3) additionally had minimal early and late time point regression coefficients, however had a substantially increased 4 hour regression coefficient.
In order to better assess the contribution of the regression coefficients to the model outcome the absolute value of the coefficients was taken for each time point and the mean plotted for each phosphoprotein in descending order (Figure 2C). From this, p-Erk was determined to most strongly contribute to the model, followed by p-RPS6 and p-JNK. We used this data to plot the R2 value of models built on increasing amounts of data, starting with p-Erk and adding phosphoproteins in order of their mean absolute value of regression coefficients. It can be seen that a model built solely on p-Erk, p-RPS6, and p-JNK resulted in R2 values of 0.4655 as compared to the complete model which gave us a R2 value of 0.616 (Figure 2D). Beyond these phosphoproteins, only the Akt phosphoprotein added substantial further information to the model, increasing the R2 from 0.484 to 0.570, indicating this data added substantial accuracy to the model without having a large regression coefficient. From these results it was concluded that the phosphorylation levels of Erk, RPS6, JNK, and Akt were able to explain the majority of variation in castration resistant survival across these three cell lines.
The amount of error between the predicted values from the model and the measured values were also grouped by treatment, cell line, and inhibitor (for LNCaP cells treated in combination with targeted inhibitors) (Additional file 5: Figure S1A, B, and C). The only significant difference that was observed between any conditions was a much higher docetaxel error (Additional file 5: Figure S1A). This is likely due to the fact that docetaxel is a chemotherapeutic which causes cell death, however little variation in the phosphoproteome as compared to controls was seen. Therefore a model of phosphoproteomic signaling was unable to predict docetaxel’s apoptotic effect.
The effect of androgen treatment on phosphoprotein signaling
The effect of DHT on phosphoprotein activation was examined across the different treatments conditions. Previous research indicates that the activated AR may act through growth factor pathways such as PI3K (Phosphoinositide 3-kinase), and by causing the transcription of genes which may directly activate the cell cycle [22]. Upon examining the DHT treatment group an increase in the 24 hour p-RPS6 and p-Akt levels as compared to controls was observed in LNCaP cells (Figure 3A). The effect of DHT on PC3 and MDA-PCa-2b cells was also examined. PC3 cells exhibited no substantial alterations in signaling which is consistent with previous reports where PC3 cells had minimal to no AR expression [23]. MDA-PCa-2b cells exhibited an increase in p-RPS6 and p-GSK3β at 4 hours which was not maintained through 24 hours, although DHT treatment of MDA-PCa-2b cells did not cause survival increases to the extent that EGF or IGF1 treatment did.
The survival of LNCaP cells in response to DHT treatment was examined and an increase of 38% was observed as compared to the control condition (Figure 3B). This survival advantage was completely abrogated when treated in combination with LY294002 (PI3Ki) which reduced p-Akt, p-GSk3, and p-RPS6 to below baseline levels at all time points. The combination of DHT plus LY294002 caused a non-significant increase in survival of 25% over the treatment of LY294002. There was little difference in phosphoprotein levels from LY294002 treatment alone, indicating direct activation of the cell cycle by AR or activation of other non-measured pathways by AR other than PI3K.
Based on these observations we propose a modification of the model originally proposed by Gosh et al. (Figure 3C) [24]. Here, the PI3K pathway can activate the AR which can activate the cell cycle. However, activation of the AR can also activate the PI3K pathway. Additionally, activation of the PI3K pathway can activate cell cycle through bypassing the AR via mTOR/RPS6.
Comparison of phosphoprotein alterations between LNCaP, MDA-PCa-2b, and PC3 cell lines
The differences between the signaling of the three different cell lines used were examined by taking the mean phosphoprotein level across all treatments, with the exception of inhibitor treatments in LNCaP cells. Several observations were noted in this data including the consistent trend across p-Akt, p-RPS6, and p-GSK3 of higher values in the LNCaP cells, somewhat reduced values in the PC3 cells, and the lowest amount of phosphoprotein in MDA-PCa-2b cells (Figure 4A). These phosphosites are part of the PI3K pathway which likely explains their similar levels of activation (the measured GSK3 phosphosites of GSK3α at S21 and GSK3β at S9 are activated by p-Akt). When p-Erk levels were measured in MDA-PCa-2b cells, consistently lower amounts of this phosphoprotein were found as compared to LNCaP and PC3 cells (10.7% of LNCaP levels and 11.3% of PC3 levels, Figure 4A). Based on the substantial weight placed on the p-Erk regression coefficient, this explains one of the major reasons for reduced castration resistance in MDA-PCa-2b cells.
A final observation made regarding the mean phosphoprotein levels across all treatments was the decreasing levels of phosphorylation in JNK from MDA-PCa-2b cells to LNCaPs and then PC3 cells (Figure 4B, C, and D). Initially, this was a counterintuitive observation due to the fact that this phosphosite has previously been described as an oncogene, and we have measured castration resistance in the cell lines inverse to the amount of p-JNK (Additional file 6: Figure S2) [25]. However, this observation corroborates recent work indicating that JNK acts as an oncogene in tumor development and a tumor suppressor in regards to castration resistant growth [19].
In order to better illustrate the activation of phosphoproteins between cell lines in response to treatments, graphs were created which plot the phosphoprotein response as a function of edge thickness (Figure 5A, B, and C). Upon examining these graphs substantial variation between the cell lines is observed with the most castration resistant cell line, PC3, having the weakest response generally to the various treatments, followed by moderate responses in LNCaP cells, and strong sensitivity to certain growth factors in MDA-PCa-2b cells. Furthermore, there were differences between the cell lines in response to the same growth factor. In PC3 and LNCaP cells EGF stimulates Erk to various extents, however in MDA-PCa-2b cells EGF had little effect on Erk and strongly increased p-RPS6 along with IGF1 which was not seen to have an effect LNCaP or PC3 cells.
Modeling the effect of treatments and targeted inhibitors
The effect of treatment with five targeted kinase inhibitors on protein phosphorylation and the LNCaP cell survival in androgen depleted media as compared to controls can be seen (Figure 6A, B, and C). Cells were treated with concentrations five times the published IC50 (half maximal inhibitory concentration) values of the target kinases which, assuming a hill coefficient of one, is equal to IC83. Some of the targeted kinase inhibitors did not reduce their target phosphoproteins to the anticipated levels, possibly due to degradation. Incomplete inhibition of targets should have no effect on model performance because the response is predicted according to actual measured phosphoprotein levels. We calculated a separate PLS regression model solely on all of the LNCaP data, including inhibitor treatments. A leave-one-out cross valuidated R2 value of 0.58 (Additional file 7: Figure S3) was observed across this data set indicating that the response from inhibitor treatment can predict the majority of the variation in cell survival.
The effect of complete PI3K inhibition with LY294002 versus mTor inhibition alone with temsirolimus was also examined. Based on the relative survival levels of LNCaP cells treated with LY294002 versus temsirolimus it was determined that the temsirolimus treated group had 31% increased cell survival over cells treated with LY294002. However, both treatments reduced the p-RPS6 to similar levels which were near complete inhibition from basal levels, while LY294002 also strongly reduced measured p-Akt and p-GSK3 levels (Figure 6D and 6E). Based on this observation it was concluded that signaling upstream of mTor (such as p-GSK3 which was observed to be highly correlated to p-Akt) accounted for the difference in survival between complete PI3K inhibition and inhibition of mTor alone.
Modeling the correlation between phosphosites’ activation
In order to better understand the correlation between different phosphoproteins’ activation under the same treatment we examined the Pearson correlation between them across the three separate cell lines (for LNCaP cells the inhibitor plus treatment data was excluded). The most consistent theme across the cell lines was the positive correlation between p-RPS6 and p-Akt, which occurs through mTor (Q-value of 0.0531, 0.0391, and 0.0160, for PC3, LNCaP, and MDA-PCa-2b cells, respectively, Figure 7). Additionally, there was a correlation between p-Akt and p-GSK3 present in LNCaP cells (Q-value of 0.00569) and MDA-PCa-2b cells (Q-value of 0.000216), but not PC3 cells (Q-value of 0.42972).