Analysis of the molecular networks in androgen dependent and independent prostate cancer revealed fragile and robust subsystems

Abstract

Androgen ablation therapy is currently the primary treatment for metastatic prostate cancer. Unfortunately, in nearly all cases, androgen ablation fails to permanently arrest cancer progression. As androgens like testosterone are withdrawn, prostate cancer cells lose their androgen sensitivity and begin to proliferate without hormone growth factors. In this study, we constructed and analyzed a mathematical model of the integration between hormone growth factor signaling, androgen receptor activation, and the expression of cyclin D and Prostate-Specific Antigen in human LNCaP prostate adenocarcinoma cells. The objective of the study was to investigate which signaling systems were important in the loss of androgen dependence. The model was formulated as a set of ordinary differential equations which described 212 species and 384 interactions, including both the mRNA and protein levels for key species. An ensemble approach was chosen to constrain model parameters and to estimate the impact of parametric uncertainty on model predictions. Model parameters were identified using 14 steady-state and dynamic LNCaP data sets taken from literature sources. Alterations in the rate of Prostatic Acid Phosphatase expression was sufficient to capture varying levels of androgen dependence. Analysis of the model provided insight into the importance of network components as a function of androgen dependence. The importance of androgen receptor availability and the MAPK/Akt signaling axes was independent of androgen status. Interestingly, androgen receptor availability was important even in androgen-independent LNCaP cells. Translation became progressively more important in androgen-independent LNCaP cells. Further analysis suggested a positive synergy between the MAPK and Akt signaling axes and the translation of key proliferative markers like cyclin D in androgen-independent cells. Taken together, the results support the targeting of both the Akt and MAPK pathways. Moreover, the analysis suggested that direct targeting of the translational machinery, specifically eIF4E, could be efficacious in androgen-independent prostate cancers.

Figure 1. Schematic overview of the interaction network used in modeling the androgen response in prostate epithelial cells.

The model architecture was formulated by aggregating molecular modules into a single network (see insert for high level details). The model describes growth factor and hormone induced expression of cyclin D, PSA and the two forms of PAcP. The complete list of molecular interactions that comprise the model (along with kinetic parameter values) are given in Table S1.

Figure 2. Identification and properties of the prostate model ensemble.

A: Steady state PSA level as a function of cPAcP and sPAcP expression. The circles represent the values used to model the C-51 and C-81 LNCaP clones. All values are relative to C-33. B: Coefficient of Variation (CV; standard deviation of a parameter relative to its mean value) for the parameter ensemble used in this study. A small CV suggested a parameter was tightly constrained by the training data used for model identification. The parameters with the three smallest CVs are listed. C: Parameter identification strategy. Multiple monte-carlo trajectories were used to randomly explore parameter space. The simulation error and the correlation between parameter sets was used to generate the family of parameter sets used in the simulation study.

Figure 3. Simulation results for the addition of 10nm DHT at 1 hour to C-33 and C-81 LNCaP clones.

A: Her2 phosphoralation (circles) and cPAcP expression (squares) for C-33 cells following the addition of DHT. Experimental data reproduced from Meng et al. . B: PSA expression following the addition of DHT to C-81 (squares) and C-33 (circles) LNCaP clones. Experimental data reproduced from Lee et al. . The shaded region in each plot denotes one standard deviation centered about the ensemble mean (line).

Figure 4. Simulated PSA mRNA levels in C-33 cells with and without Her2 overexpression.

Her2 overexpression was modeled as a 50% increase in the expression rate of Her2. Bars denote the mean PSA mRNA level over the parameter ensemble while error bars denote one ensemble standard deviation. The experimental PSA mRNA data was adapted (replotted) from .

Figure 5. Simulation results for key species under androgen free conditions.

A: Effect of HER2 and MEK overexpression on LNCaP C-33 steady state PSA levels. The inhibition of MEK blocks the effect HER2 overexpression. Experimental data adapted from Lee et al. . B: Effect of HER2 and MEK inhibition on LNCaP C-33 steady state PSA levels. The inhibition of either HER2 or MEK blocks high AIPC PSA levels. Experimental data adapted from Lee et al. . C: Effect of PAcP isoforms on LNCaP steady state cyclin D levels. Experimental data adapted from Lingappa and coworkers (Prosetta Corporation, unpublished data). D: Transient activation of ERK via ligand dependent EGF signaling (8nM EGF at t = 60s) in HeLa cells. The HeLa data was reproduced from . Inset: Simulated phosphorylated ETS (ETSp) levels following the addition of 8nM EGF in the presence and absence of Her2. Her2 activation drives a sustained MAPK signal which in turns sustained ETS activation. The shaded region denotes one standard deviation centered about the ensemble mean (line).

Figure 6. Independent model predictions versus experimental observations.

A Ensemble prediction of cyclin D expression following the addition of DHT at 1 hour to C-33 clones. The ensemble predicted a dose dependent increase of cyclin D at 24 hours after DHT addition. Experimental data was adapted from Barnes-Ellerbe et al. . B Predicted effect of an AR knockdown on PSA expression following the addition of androgen at 1 hour to C-33 wild-type and C-33 AR knock-down clones. The ensemble predicted an approximate 50% decrease in androgen stimulated PSA expression due to AR knock-down 72 hours after treatment. Experimental data was reported by Eder et al. . The error bar denotes one standard deviation centered about the ensemble mean.

Figure 7. Sensitivity analysis of the model parameters.

A: Comparison of the mean OSSC parameter ranks for the C-33 and C-81 LNCaP models. Large ranks indicate fragility. Points left of the 45 line are more important in C-33, while shifts to the right show increased importance in C-81. Points are organized by biological function. B: Comparison of the mean OSSC parameter ranks for translation mechanisms (including the role of Akt signaling in translation initiation) in C-33 versus C-81 LNCaP clones. The error bars indicate one standard deviation centered about the mean ensemble value. C: The final mechanism in PSA transcription becomes increasingly more robust w.r.t cancer aggressiveness, as indicated by a significant reduction in mean OSSC Rank. D: The final mechanism in PSA translation (translation termination) was increasingly fragile w.r.t cancer aggressiveness, as indicated by a significant increase in mean OSSC rank. The results indicate a shift in the bottle neck for generation of PSA from transcription to translation as prostate cancer cells lose their androgen dependence. The top and bottom of each box denote the 25th and 75th percentile of the OSSC rank over the parameter ensemble. The center line denotes the median value. Whiskers show the furthest observations and black crosses indicate outliers.

Figure 8. Robustness analysis of functional protein markers.

The expression level of seven key proteins was altered by a factor of 10, .1 or 0 (knock-in, knock-down or knock-out) and robustness coefficients (area under the curve for the perturbed versus wild-type simulation) were calculated for cyclin D and PSA expression levels along with ERK and AR activation levels. Simulations were run for C-81, with the indicated perturbation, to approximate steady-state and 10nM of DHT was added for 72 hours. Ensemble mean values are reported.

Figure 9. Synergy analysis between the ERK and Akt signaling axes in LNCaP C-81 cells. The double ERK and Akt knock-out was used as the control.

A: The difference in steady state cyclin D expression (compared to the control) with the knock-in of Akt (left), ERK (center) and both (right). The predicted cyclin D levels were normalized by the basil C-81 steady state cyclin D level in each case. The error bars denote one standard deviation centered about the ensemble mean. The region denoted by the asterisks represents above-additive cyclin D expression. B: Species and interactions that demonstrated a positive (negative) synergy are shown as green (red) in the connectivity diagram. Species or interactions not effected are shown in grey. C: The full connectivity diagram qualitatively clustered in functional groups. Positive (negative) synergy are shown in green (red) in the connectivity diagram. Species or interactions not effected are shown in grey.