Emergent dynamics of underlying regulatory network links EMT and androgen receptor-dependent resistance in prostate cancer

Abstract

Advanced prostate cancer patients initially respond to hormone therapy, be it in the form of androgen deprivation therapy or second-generation hormone therapies, such as abiraterone acetate or enzalutamide. However, most men with prostate cancer eventually develop hormone therapy resistance. This resistance can arise through multiple mechanisms, such as through genetic mutations, epigenetic mechanisms, or through non-genetic pathways, such as lineage plasticity along epithelial-mesenchymal or neuroendocrine-like axes. These mechanisms of hormone therapy resistance often co-exist within a single patient's tumor and can overlap within a single cell. There exists a growing need to better understand how phenotypic heterogeneity and plasticity results from emergent dynamics of the regulatory networks governing androgen independence. Here, we investigated the dynamics of a regulatory network connecting the drivers of androgen receptor (AR) splice variant-mediated androgen independence and those of epithelial-mesenchymal transition. Model simulations for this network revealed four possible phenotypes: epithelial-sensitive (ES), epithelial-resistant (ER), mesenchymal-resistant (MR) and mesenchymal-sensitive (MS), with the latter phenotype occurring rarely. We observed that well-coordinated "teams" of regulators working antagonistically within the network enable these phenotypes. These model predictions are supported by multiple transcriptomic datasets both at single-cell and bulk levels, including in vitro EMT induction models and clinical samples. Further, our simulations reveal spontaneous stochastic switching between the ES and MR states. Addition of the immune checkpoint molecule, PD-L1, to the network was able to capture the interactions between AR, PD-L1, and the mesenchymal marker SNAIL, which was also confirmed through quantitative experiments. This systems-level understanding of the driver of androgen independence and EMT could aid in understanding non-genetic transitions and progression of such cancers and help in identifying novel therapeutic strategies or targets.

Keywords: Androgen independence; Epithelial-Mesenchymal Transition; Multistability; Non-genetic heterogeneity; PD-L1; Phenotypic plasticity; Snail.

Graphical abstract

Fig. 1

Multi-stable dynamics of the coupled EMT–AR crosstalk network. A) (i) Gene regulatory network (GRN) showing crosstalk between EMT and Androgen Receptor (AR, AR-v7) signalling. Blue arrows represent activation links; red hammers represent inhibition. (ii) Heatmap of stable steady-state solutions for network shown in A (i), as obtained via RACIPE. B) (i) Kernel Density Estimate plots with histograms of z-score levels of individual nodes in a network, fitted to Gaussian distributions. Each panel legend shows the corresponding Bimodality Coefficient (BC). (ii) Pairwise correlation matrix in which each cell denotes the correlation coefficient between the corresponding set of genes. Red indicates positive correlation; blue indicates negative correlation; colormap indicates the strength and sign of correlation; cells with a cross (x) highlight non-significant correlation (p > = 0.05). C) (i) Histogram of Resistance score (AR + AR-V7) fit to a Gaussian. (ii) Histogram of epithelial-mesenchymal (EM) score (ZEB1 - miR200) fit to a Gaussian. (iii) Scatter plot showing the corresponding EM and resistance scores for all RACIPE solutions obtained. D) (i) UMAP (uniform manifold approximation and projection) dimensionality reduction plots for steady-state solutions colored by Resistance score; colormap indicates the Resistance score. (ii) Same as (i) but for EM score. (iii) UMAP based on EM and resistance scores together; Quadrant 1 represents mesenchymal resistant, Quadrant 2 represents epithelial resistant, Quadrant 3 represents epithelial sensitive and Quadrant 4 represents mesenchymal sensitive populations.

Fig. 2

Coordinated EMT and AR programs in network topology and transcriptomic data. A) Influence matrix for the biological network shown in Fig. 1A. B) Histogram of group strengths for 100 random (non-biological) networks generated by shuffling edges in the biological network. The group strength of the biological network us shown in red C) (i) Scree plot showing PCA components and their explained variance percentage (ii) PCA Correlation Circle of the RACIPE solutions of the biological network D) Scatter plots indicating Spearman's correlation between metrics of EMT (EMT KS score, EMT 76GS score, ZEB1 and SNAI2 expression levels) and ssGSEA scores for Wang Prostate Cancer Androgen Independent gene set, in TCGA (above) and GSE74685 (below). E) Pairwise correlation between different EMT scoring and androgen independence gene lists/metrics in different transcriptomic datasets (Left to right - GSE77959 (n = 30), GSE80042 (n = 44), GSE67681 (n = 6), TCGA (n = 551) and GSE22010 (n = 12)). Labels: 1 – ssGSEA scores for Wang Prostate Cancer Androgen Independent geneset (MsigDB), 2 – ZEB1 expression levels, 3 – KS score, 4 – ssGSEA scores for Hallmark EMT geneset (MsigDB), 5 – 76GS score. Color bar indicates Spearman's correlation coefficient. Crosses indicate corresponding p-value> = 0.05.

Fig. 3

Stochastic cell-state transitions along EMT and AR axes. A) (Top) Pie chart showing fraction of RACIPE parameter sets in terms of number of steady-state solutions enabled. (Bottom) Pie chart showing combinations of different phenotypes constituting bistable solutions. B) System dynamics for representing two biological EMT phenotypes (E, M) and Resistance (R, S) when starting from multiple different initial conditions. C) (Top) Stochastic cell-state transition from MR to ES phenotypes under the influence of gene expression noise. (Bottom) A zoomed-in version of the highlighted region. D) Landscape showing log (likelihood) on the Sensitivity and EM score planes, with the valleys representing the stable states possible in the system - epithelial sensitive, mesenchymal resistant E) EMT scores (KS scores, 76GS scores and ssGSEA scores for the Hallmark EMT gene set) for a single-cell RNA-seq dataset (GSE168668) comprising samples (cells) treated with DMSO, treated with enzalutamide for 48 h (ENZ48) and enzalutamide resistant (RES) cells. * ** *: p < 0.01 for Students’ t-test.

Fig. 4

Association of PD-L1 with EMT/AR crosstalk. A) RNA-Seq data from The Cancer Genome Atlas shows a significant positive correlation between SNAIL and PD-L1 expression in prostate cancer samples B) mRNA expression of PD-L1 with and without SNAIL in the LNCaP95-SNAIL inducible cell line C) Protein expression of SNAIL and PD-L1 in LNCaP95 with and without SNAIL activation D) Gene regulatory network showing the androgen receptor axis and EMT related genes, along with an added PD-L1 node E) Heatmap of stable steady-state solutions for the network shown in A, obtained via RACIPE F) Correlation matrix for the genes in the network; red indicates positive correlation and blue indicates negative correlation; x indicates correlations with p > 0.05 and the color of each cell indicates the strength of the correlation. G) Western blotting for Snail in paired enzalutamide-sensitive (S) and –resistant (R) cell lines. GAPDH was included as a loading control. H) PD-L1 protein level quantification based on phospho-protein arrays. * indicates p < 0.05.

Fig. 5

Association of PD-L1 signature with AR in simulation and transcriptomic data. A) (i) Scatter plots indicating the Pearson’s correlation between PD-L1 signature (calculated based on ssGSEA scores of 15 genes defined in [65]) and SNAI1 levels in various clinical data sets and (ii) Scatter plot indicating the Pearson’s correlation between PD-L1 and SNAI1 levels at steady state produced by simulations B. (i) Scatter plots indicating the Pearson’s correlation between PD-L1 and AR levels in various clinical data sets and (ii) Scatter plot indicating the Pearson’s correlation between PD-L1 and AR levels at steady state produced by simulations.

Fig. 6

Meta-analysis showing association between EMT, androgen independence and PD-L1 in bulk transcriptomics data. A) Volcano plots showing Spearman correlation coefficients (x-axis) and -log₁₀(p-value) (y-axis) for Epi vs. Mes (top) and Hallmark EMT vs. Mes (bottom). Significant correlations (R>± 0.3 and p < 0.05) shown as red (positive correlation), blue (negative correlation) datapoints. Each dot represents a unique data set. Numbers on each side of the volcano plots indicate the number of data sets that are positively correlated (red) or negatively correlated (blue). Same as A) but for B) Hallmark EMT vs. PD-L1 (top), Mes vs. PD-L1 (bottom), C) Hallmark EMT vs. Androgen Independence (top), Mes vs. Androgen independence (bottom).