Mechanisms of Resistance to Intermittent Androgen Deprivation in Patients with Prostate Cancer Identified by a Novel Computational Method

In this work, we limit our attention to prostatic adenocarcinomas of non-neuroendocrine origin as this entity has by far the greatest clinical significance of all prostate tumors. Healthy prostatic epithelial cells and most malignant prostatic epithelial cells exhibit some degree of dependence on androgens, including testosterone and dihydrotestosterone (DHT), for proliferation and survival (1–5). The contribution androgens provide to cell proliferation and survival is mediated by the androgen receptor (AR) axis. It is well known that AR signaling regulates expression of genes involved in proliferation and survival and is generally accepted to be the primary driver of tumor progression in adenocarcinomas (1, 2, 4–7). AR signaling inhibition, therefore, regresses tumors (2, 4, 8, 9).

The most common and effective treatment for progressive prostate cancer is androgen-deprivation therapy (5, 6). The mechanism of action of this therapy is largely induction of apoptosis, contingent on prostate cancer cells being dependent on androgens, via the AR axis, for survival (2–5). Androgen deprivation can be accomplished by surgical castration (orchiectomy) or chemical castration (with drugs). The latter method is commonly carried out by “total androgen blockade,” which combines luteinizing hormone-releasing hormone (LHRH) analogues and antiandrogens (10). Chemical castration is almost universally effective at regressing tumors initially; however, almost all patients eventually develop resistance to treatment and advance to a more aggressive stage, for which there is no effective treatment (5, 6, 10, 11). This advanced stage has traditionally been referred to as androgen-independent prostate cancer, the rationale for this terminology being that if prostate cancer cells have no androgen access, but continue to survive and proliferate, then they must have bypassed the need for androgens. However, although chemical castration depletes blood serum testosterone levels by >90%, intraprostatic androgen concentrations remain at 20% to 50% (2, 4, 5, 12). These residual intratumoral androgen sources allow prostate cancer cells to survive because these cells are still dependent on androgens for survival (5). The recent clinical efficacy of abiraterone (which inhibits intraprostatic androgen synthesis) and enzalutamide (a second-generation AR antagonist), has further confirmed that these cells are not androgen independent (5, 13, 14). Therefore, this stage of cancer is now generally referred to as castration-resistant prostate cancer (CRPC; refs. 2, 4, 5, 12, 15).

Because continuous androgen deprivation (CAD) provides a prolonged stimulus selecting for castration-resistant clones, intermittent androgen deprivation (IAD) was proposed to delay the onset of CRPC (10). Although this basis for IAD remains controversial (16), IAD is cheaper and minimizes side effects relative to CAD (10). For detailed clinical examples of the IAD protocol, see refs. 9, 17, 18.

The gene coding for prostate-specific antigen (PSA) is exclusively activated by AR signaling in prostatic epithelial cells (1, 6). PSA is secreted by these cells into the blood serum (1). Because AR signaling is considered the primary driver of tumor progression, before and after androgen-deprivation failure, PSA is a commonly used predictive biomarker for monitoring androgen-deprivation efficacy (1, 9). Increasing PSA levels are characteristic of tumor progression, whereas decreasing PSA levels are characteristic of tumor regression. For IAD, patients are usually put on treatment until their PSA levels fall below some threshold value, usually <4 ng/mL, and taken off-treatment until PSA rises above some threshold level, usually >10 ng/mL (10).

Several mathematical models have been developed with the goal of fitting model output to clinical PSA data throughout the course of IAD (19–23). Here, we address how a mechanistic model of prostate cancer that takes an individual patient's serum androgen data as input, and uses data fitting to clinical PSA data to estimate parameters, might be used to inform clinicians of the inherent mechanism(s) of treatment resistance in individual patients.

Although most prostate cancer cells exhibit some degree of dependence on androgens for survival, as tumors continue to advance to a more progressive state, the degree of dependence tends to decrease. Several pathways to CRPC development have been identified with varying degrees of androgen dependence (see ref. 1 for detailed review of these mechanisms). Later, we attempt to quantify this variation in dependence and use it to identify which mechanism is developing in individual patients. Clinical identification of these mechanisms in individual tumors still suffers from a lack of maturity because CRPC tissue is difficult to obtain and study (11). Hence, mathematical models may be a potentially powerful clinical tool for analyzing mechanisms of CRPC development. Here, we highlight mechanisms that lead to CRPC, which we propose can be identified through our model.

Cells using this as their major mechanism have a lower androgen threshold for proliferation and survival but remain completely dependent on androgens. We consider three specific pathways to hypersensitivity: amplification and overexpression of AR, increased AR stability, and increased local androgen production. Overexpression of AR is characteristic of most CRPC cases (7, 15, 24). It has been demonstrated that overexpression of wild-type AR is adequate to fully transform a CS tumor into a CR tumor (24). AR amplification is one of the most common genetic aberrations to occur in CRPC (15). It has been suggested that 80% of CRPC cases exhibit elevated AR gene copy number whereas 30% exhibit high-level amplification (15). It has also been suggested that increased AR stability may be a mechanism leading to CRPC (1). Increased local androgen production can result from constitutive activation of 5α-reductase, which converts testosterone to DHT in prostate cells. For example, the V89L mutation in the gene coding for 5α-reductase leads to an increase in its activity (1). Furthermore, some CRPC cells may be able to synthesize the androgens themselves from cholesterol or other androgenic precursors (1, 5, 6).

Promiscuous ARs generally arise from missense mutations in the ligand-binding domain that decrease the specificity of the AR for a particular ligand. These mutations broaden the group of ligands that activate the AR, such that the AR obtains the ability to bind to and be activated by other nonandrogen agonists in addition to androgens (1, 4). Potential nonandrogen agonists include estrogen, progesterone, cortisol, antiandrogens, or adrenal androgens such as DHEA (4). For example, the most common promiscuous AR mutation is T877A, which is endogenously expressed in androgen-sensitive LNCaP cell lines in which it was first discovered (4, 12). Taplin and colleagues found the T877A mutation to be present in 30% of prostate cancer bone metastases (25). The T877A-mutant AR still binds to, and is, thus, activated by, testosterone and DHT. However, alanine replaces the threonine residue at position 877 and allows ligands other than testosterone and DHT to also fit in the binding pocket and activate the receptor (12). Promiscuous ARs with complete loss of responsiveness to androgens are extremely rare (4). Although mutations have been identified that reduce affinity for DHT relative to the wild-type AR, nearly all identified promiscuous ARs are still activated by DHT (4). We, therefore, assume that all CR cells with promiscuous ARs retain responsiveness to androgens; however, they are less dependent on these androgens than cells with hypersensitivity.

Outlaw receptors are steroid hormone receptors that are activated by ligand-independent mechanisms. No mutations in the AR gene coding region have been identified such that the AR is activated by ligand-independent mechanisms (1). Rather, this pathway apparently involves the wild-type receptor or AR splice variants, because the AR gene sequence is not changed, and may still exhibit responsiveness to androgens. There have been several AR splice variants identified such that the ligand-binding domain is absent and these ARs remain ligand independently constitutively active (6, 26). However, exclusive expression of AR splice variants is very rare. The vast majority of cells expressing AR splice variants also express wild-type AR (7). Therefore, most tumors developing resistance through outlaw AR splice variant expression should still exhibit responsiveness to androgens while maintaining active AR signaling in the absence of androgens. It has been shown the wild-type AR can be activated by growth factors, receptor tyrosine kinases, MAPK, AKT, and other activators in a ligand-independent fashion (1). If no ligands are available other than androgens in a cell with promiscuous ARs, the cell then depends on androgens. If no ligands are available in a cell with outlaw AR, AR signaling remains active. Thus, this pathway is less dependent on androgens than the promiscuous receptor pathway because it is completely ligand independent.

This mechanism results from constitutive activation of pathways that directly inhibit apoptosis in the absence of androgens (1). One mechanism may be the upregulation of BCL2, BCLX, or MCL1, genes known to directly inhibit apoptosis (1, 4, 10). Prostatic secretory epithelial cells normally do not produce BCL2, but many CRPC cells do (1, 10). However, besides prostatic sarcomas and neuroendocrine carcinomas, which our model does not consider, the complete bypass of AR signaling does not provide a selective advantage (4, 6). This pathway is the least dependent on androgens for survival.

We extend the prostate cancer IAD treatment model by Portz and colleagues (19) who consider tumors as two cell populations: cells that are castration sensitive (CS) and those that are castration resistant (CR). The mathematical model portrayed in the Quick Guide describes two populations of luminal secretory epithelial cells that secrete PSA (1). The CR population accounts for prostate cancer cells that have acquired at least one mechanism of treatment resistance. The static parameter values, free parameter value ranges, and their interpretations are listed in Table 1. Detailed model formulation and explanations can be found in (19) and are summarized in the Quick Guide above. What follows is a narrative summary of the equations and assumptions.

Here, we follow Portz and colleagues in assuming that proliferation of prostate carcinoma cells is more-or-less androgen dependent, this dependence can vary among cells in a given tumor, and cells can switch reversibly between CS and CR states. Our major extension of the Portz and colleagues' model is the addition of an androgen-dependent cell-death rate (CDR) function, |$D_i (Q_i)$| (Equation C). We retain the Droop formalism (27), adapted from chemostat and ecologic models, to govern androgen-dependent proliferation, as well as the forms used by Portz and colleagues for androgen-dependent switching between CS and CR states—i.e., λ₁ and λ₂ (Equations A, B, D, and E).

In this model, androgen enters cells in which it can bind to and dissociate from ARs in the intracellular compartment. Androgen may also be metabolized or transported out of the system. The dynamics of these processes are governed by the cell quota (Equation F). Bound ARs alter cell proliferation, death and state-switching rates. In general, proliferation rates increase and mortality rates decrease with increasing concentration of bound ARs. The total AR concentration, both bound and unbound, is regulated in the cell quota equation such that the proliferation terms in Equations A and B can never become negative. In low androgen environments, state switching by prostate cells is biased in favor of the CR state. This bias is reversed in high androgen environments. Both cell strains secrete PSA to the serum at the same basic rate. However, strains may vary in how their PSA secretion rate depends on bound AR (see equation G).

It is an open question whether castration resistance evolves from selection acting on a small population of CR cells that exist in the tumor before treatment (the selection hypothesis), or whether CR strains arise during treatment (the adaptation hypothesis), or both (28, 29, 30). The role accommodative changes play is also not well understood. Here, by accommodative changes, we mean reversible changes in cell behavior not associated with any genetic or genomic alterations—that is, plasticity in the cells' behavioral repertoire available to accommodate fluctuations in their environment. The possibility of accommodative switching between relative CS and CR states in a single clonal lineage was initially investigated in a model by Hirata and colleagues (20, 21). Portz and colleagues (19) assumed reversible switching between CS and CR states at rates λ₁ and λ₂, which both depend on the bound AR cell quota (see Equations A, B, D, and E). Although Portz and colleagues interpreted these as forward and backward mutations between CS and CR states, the form is much more general. Any mechanism of castration resistance, including mutation in coding and control regions of a variety of genes and epigenetic changes, can be accommodated by the model's form. Regardless of which process the mathematics truly reflects, the incorporation of this switch resulted in Portz and colleagues' model fitting data much more accurately than previous models from which it was derived. However, we interpret the switch as incorporating both adaptation and accommodation to cover all possible switching mechanisms. The adaptation label is supported by the fact that AR point mutations are rarely found in primary prostate cancer but are found in approximately 20% of CRPCs (11). Furthermore, one study showed that only two of 205 untreated prostate cancer cases had detectable AR amplification (31), and another study showed that 38% to 63% of metastatic CRPC cases exhibited high-level AR amplification (32, 33). This evidence strongly suggests that an oscillating androgen environment creates selective pressure inducing an adaptive switch. Although this evidence only shows evolution toward CRPC, it is possible that CR tumors could evolve back to CS if selective pressures change (e.g., androgen is reintroduced); although this return to castration sensitivity would be unlikely to be due to direct back mutation, we recognize that the phenotype can be regenerated by novel alterations of coding or control sequences to a variety of genes in a given resistance pathway. It is necessary to include accommodation into the interpretation because epigenetic modifications, which can lead to phenotypic switching, take place at a much higher rate than genetic modifications (34) and several epigenetic aberrations play a role in prostate cancer that can facilitate the phenotypic switch represented by λ₁ and λ₂ (35).

The CDR analysis method we propose here consists of analyzing the CDR levels and oscillations in response to intermittent treatment to determine which molecular basis is conferring treatment resistance in a given patient. The rationale behind CDR analysis is as follows:

1. The mechanism of action of IAD is largely induction of apoptosis, because all CS cells and most CR cells exhibit some degree of dependence on androgens for survival.

2. Because the treatment is intermittent, CDRs increase during on-treatment and decrease during off-treatment giving rise to CDR oscillations.

3. The amplitude of CDR oscillation in a given population reflects the degree to which that population depends on androgens for survival.

4. Each mechanism of treatment resistance varies in its dependence on androgens for survival.

5. In addition, CDR levels tell us whether or not there is a constitutively active antiapoptotic pathway inherent in the population (i.e., the bypass pathway). An exceedingly low CDR indicates there is.

6. We identify which major mechanism of treatment resistance is developing by analyzing the CDR levels and oscillations.

The data from Akakura and colleagues were obtained from seven patients with progressive prostate cancer undergoing IAD by total androgen blockade (9). Serum PSA and testosterone were assayed monthly to monitor therapy. Treatment was initiated and subsequently interrupted after serum PSA had been suppressed to adequately low levels for roughly six months. When PSA recovered past 15 to 20 ng/mL, treatment was reinitiated for a new cycle.

Treatment induced rapid decreases in serum testosterone, with levels reaching approximately 0.5 to 1.0 nmol/L after about one or two months. These low levels were sustained until interruption of treatment, after which testosterone recovered relatively quickly, with some variation depending on the patient. The coupling of these changes to those of PSA varied across not only patients but also cycles of individual patients. Akakura and his colleagues suggested that these dynamics were due both to changes in androgen sensitivity and to the “dual effect” of androgen on PSA; namely, the cellular level (PSA synthesis and secretion) and the population level (proliferation and apoptosis; ref. 9). The model captures these effects via the dependence of Q(t) (Equation F) on the serum androgen concentration A(t). CDR analysis quantifies this behavior across cycles and patients in a way that is both simple and biologically meaningful.

For each patient, CDR values were calculated by fitting the model to the PSA data. To implement the serum androgen function, A(t), testosterone data were interpolated using the same approach found in Portz and colleagues (19). The Nelder–Mead simplex algorithm (36) was then used to minimize the mean-square error (MSE) of the PSA data with the model solution P(t) by varying the parameters μ_m, σ₀, q_i, σ_i, d_i, R_i, and δ_i, |$i\, = \,1,2$|⁠. The remaining model parameters were held fixed across patients (see Table 1 for ranges and values).

To identify the major mechanism of resistance for each patient, we establish criteria based on F_i, the total CDR oscillation amplitude defined in Equation H. Because the amplitude of CDR oscillation reflects the degree to which a population depends on androgens for survival, and because each mechanism varies in its degree of dependence on androgens for survival, we propose that there exist threshold CDR oscillation amplitudes that distinguish one mechanism from another. Specifically, we define two threshold values, T₀ and T₁. For |$0\, < \,F_i \, \le \,T_0$|⁠, the CDR is said to be relatively constant and the population is, thus, relatively independent of androgens for survival. For |$F_i \, = \,0$|⁠, which occurs when |$d_i \, = \,0$| or |$R_i \, = \,0$| (⁠|$Q_{i,{\rm max}}^\alpha$| and |$Q_{i,{\rm min}}^\alpha$| are never equal), the CDR is constant and the population is completely independent of androgens for survival. For |$T_0 \, < \,F_i \, \le \,T_1$|⁠, the CDR is said to slightly oscillate with respect to oscillation in androgen concentration and the population, thus, retains some dependence on androgens for survival. For |$T_1 \, < \, F_i$|⁠, we say that there is relatively high oscillation in the CDR with respect to oscillation in androgen concentration and the population is, thus, completely dependent on androgens for survival. We also define a threshold L such that for |$\max \{D_i (Q_i)\} \, < \,L$|⁠, the CDR is said to be exceedingly low. An exceedingly low CDR would be characteristic of a constitutively active antiapoptotic pathway such that apoptosis is directly inhibited in the population, i.e., the bypass pathway.

As a novel approach for computational molecular diagnosis, our threshold values for |$T_0$|⁠, |$T_1$|⁠, and L are estimates based on existing hypotheses and data relevant to this work. Model parameters were constrained to biologically relevant ranges when applicable. When estimating these values, we considered consistency with patient-specific clinical observations and with the relative frequency of each mechanism reported in the literature. The threshold values were estimated to be |$L = 0.013\,{\rm day^{-1}}$|⁠, |$T_0 = 0.0003\,{\rm day^{-1}}$|⁠, and |$T_1 = 0.0013\,{\rm day^{-1}}$|⁠. Mutations in the AR gene that confer hypersensitivity are the most frequent mutations found in CRPC and there is strong evidence that AR overexpression is the main mechanism of CRPC development (15). Therefore, these threshold values were estimated such that 4 of 7 patients would be diagnosed with hypersensitivity. Thus, T₁ was estimated to be slightly below the CDR oscillation amplitude value of patient 7, |$F_2 \, = \,0.0014\,{\rm day^{-1}}$|⁠. Because patient 4 was the only patient proposed to have early onset of CRPC and had an exceedingly low CDR relative to others, L was estimated to be slightly above |$\max \{D_2 (Q_2)\} \, = \,0.0110\,{\rm day^{-1}}$|⁠, the maximum CDR of patient 4. T₀ was estimated to be slightly above the CDR oscillation amplitude value of patient 2, |$F_2 \, = \,0.0001\,{\rm day^{-1}}$|⁠, because patient 2 exhibited a more aggressive treatment resistance compared with others.

Although a given patient may develop CRPC due to a combination of mechanisms, we assume that it is primarily due to one major mechanism. If more than one mechanism is inherent in the population, the one that is the least dependent on androgens for survival is the major mechanism. Because the hypersensitivity pathway maintains the strongest dependence on androgens, we propose that patients developing resistance primarily through this pathway should exhibit a high oscillation amplitude F₂ with a maximum CDR that is not exceedingly low; i.e., |$T_1 \, < \,F_2$| and |$\max \{D_2 (Q_2)\} \, \ge \,L$|⁠. The promiscuous receptor pathway, which is less dependent on androgens than the hypersensitivity pathway, should exhibit a smaller oscillation amplitude F₂ while maintaining a maximum CDR that is not exceedingly low; i.e., |$T_0 \, < \,F_2 \, \le \,T_1$| and |$\max \{D_2 (Q_2)\} \, \ge \,L$|⁠. We propose that the outlaw receptor pathway is less dependent on androgens than the promiscuous receptor pathway because it is completely ligand independent and, therefore, exhibits an oscillation amplitude that is zero or relatively constant and the maximum CDR should remain not exceedingly low; i.e., |$0\, \le \,F_2 \, \le \,T_0$| and |$\max \{D_2 (Q_2)\} \, \ge \,L$|⁠. The bypass pathway exhibits an exceedingly low CDR and is completely androgen independent for survival and, thus, oscillates minimally; i.e., |$0\, \le \,F_2 \, \le \,T_0$| and |$\max \{D_2 (Q_2)\} \, < \,L$|⁠. We propose no mechanism for exceedingly low and oscillating CDRs as they are biologically unlikely to exist because the apoptotic program would be directly shut down. These criteria are outlined in Table 2.

For each patient, we computed the amplitudes of CDR oscillation F_i given in Equation H and the maximum CDRs for CDR analysis. We then diagnosed each patient using the outlined criteria (see Table 2). Figure 1 shows the PSA, cell quota, CDR, and proliferation-rate plots for patient 1. Figure 2 shows the PSA plots for patients 2 to 7. Figure 3 shows the CDR plots for patients 2 to 7. The proliferation-rate and cell-quota plots for patients 2 to 7 are found in the Supplementary Material, as are the MSE values of the fitted model solutions and PSA data. Table 3 contains the predictive diagnosis for each patient based on our criteria for CDR analysis. As expected, none of the patients have cell populations with exceedingly low and oscillating CDRs.

As expected, our results suggest that most CR populations retain at least some degree of dependence on androgens for survival. Patient 4 is the only patient whose PSA level never dropped below 4 ng/mL within the first 8 months of treatment (see Fig. 2). Akakura and colleagues (9) concluded that IAD is likely ineffective for patients whose PSA levels do not fall below this stable level within the first 8 months of treatment. They suggest that failure for PSA levels to normalize within this time period is a sign of early onset of CRPC (9). Our results show that patient 4 is the only patient that did not have a highly oscillating CS CDR, meaning that the CS population for patient 4 was not fully dependent on androgens for survival, consistent with the notion of early onset of CRPC. Patients whose CS CDR oscillations are inside the range of T₁, or |$F_1 \, < \,T_1$|⁠, may, thus, be developing early onset of CRPC. Furthermore, patient 4 was diagnosed with the most aggressive and least androgen-dependent pathway according to our criteria, which is also consistent with early onset of CRPC.

In future studies, it will be important to assess how much data are required for CDR analysis. In this study, diagnoses were based on the entire datasets from the clinical trial, which varied from 1.5 to 3.5 cycles per patient. At this stage of its development, we believe at least one full cycle of treatment is needed. This includes data from start until shortly after the second on-treatment has begun and PSA has regressed to physiologic range. Ideally, physicians will continually test patients' levels and update their diagnosis profile with subsequent targeted treatment.

Although the threshold values L, T₀, and T₁ remain to be validated clinically, even once these values are empirically determined it is possible that they will still be subject to variation among datasets, and among patients. As with most models used for clinical prediction, it is likely to require repeated updating and validation with new datasets (37). Furthermore, we expect such updates not only to involve these thresholds, but also to include modifications to the model itself. We recognize this study as one that is preliminary but that establishes the principles of CDR analysis—a novel computational method for molecular diagnosis—on which others may build.

Although we demonstrated how CDR analysis works using CRPC development as an example, this methodology may be applied to other types of cancer. In particular, other steroid hormone–driven malignancies, such as estrogen receptor–positive (ER⁺) breast cancers, share many similarities in terms of their progression. For example, certain ER⁺ breast cancers can be treated intermittently with LHRH analogues to suppress serum estrogen levels in combination with antiestrogens that directly inhibit ER signaling (38, 39), which is similar to the combination of LHRH analogues and antiandrogens in the total androgen blockade used for prostate cancer IAD.

Personalized therapy on the basis of identifying and targeting specific mechanisms of treatment resistance developing in individual patients is likely one of the most promising methods of future cancer therapy. CDR analysis is a method for identifying the mechanism, which physicians can then conceivably target, thus reestablishing treatment sensitivity in these patients. Sequential identification and targeting of resistance mechanisms may potentially control or eradicate cancers.

No potential conflicts of interest were disclosed.

Conception and design: J.D. Morken, A. Packer, R.A. Everett, Y. Kuang

Development of methodology: J.D. Morken, A. Packer, R.A. Everett, Y. Kuang

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): Y. Kuang

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J.D. Morken, A. Packer, R.A. Everett, J.D. Nagy, Y. Kuang

Writing, review, and/or revision of the manuscript: J.D. Morken, A. Packer, R.A. Everett, J.D. Nagy, Y. Kuang

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Packer

Study supervision: J.D. Nagy

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.