Estrogen Receptor Pathway Activity Score to Predict Clinical Response or Resistance to Neoadjuvant Endocrine Therapy in Primary Breast Cancer

Endocrine therapy is one of the mainstays in treatment of both early and metastatic breast cancer. Especially, the use of endocrine therapy has resulted in increased survival rates (1–3). Patients are currently selected for endocrine therapy using IHC analysis of estrogen receptor (ER) and progesterone receptor (PR) expression (4). Both the American Society for Clinical Oncology and the European Society for Medical Oncology advise a threshold of 1% ER-positive tumor cells (5, 6). In practice, many clinicians and countries choose a threshold of 10% (7). More quantifiable analyses like Allred scoring and H-score have been developed and suggested for clinical application, but are currently not routinely used (5).

Despite the success of endocrine therapy in ER-positive breast cancer, up to half of patients do not show the expected response (8, 9). In addition to cancer tissue heterogeneity, several mechanisms have been proposed to explain lack of therapy response, like emergence of ESR1-activating mutations or activation of other signal transduction pathways upon pharmacologic inhibition of the ER pathway (10, 11). Standard IHC analysis detects the presence of nuclear ER protein. However, to what extent positive nuclear ER staining indicates actual functional activity of the ER pathway has not been addressed satisfactorily. A test to predict response to endocrine therapy based on measuring functional activity of the ER pathway is expected to improve decision making regarding endocrine therapy and/or alternative therapies (12).

A number of data-driven RNA-based tests have been developed to assess recurrence risk and predict response to endocrine therapy (8, 12–14). However, no test is available to measure functional ER pathway activity. A knowledge-based Bayesian network computational model for the ER pathway has been developed to assess functional activity of this pathway in tumor tissues (15, 16). The model used mRNA expression levels of 27 high evidence target genes of the ER transcription factor to infer the activation state of the ER pathway. Earlier analysis using this ER pathway model showed that only part of ER-positive patients had an active pathway, which was associated with lower risk of relapse after adjuvant tamoxifen treatment (15). Here, we report results on use of this ER pathway model to predict and assess response to endocrine therapy in patients with ER-positive breast cancer treated in a neoadjuvant setting.

Stimulation experiments with 17β-estradiol (E2) and RNA isolation of a panel of breast cancer ER-positive cell lines were performed by BioDetection Systems (BDS). Cell lines were purchased by the ATCC, and maintained in DMEM/F12 supplemented with 10% FBS. Experiments were performed at passage n+4 (CAMA1 and BT474), n+12 (MCF7), and after thawing at an unknown passage for T47D. No extra tests, besides ATCC authentication, were performed. Two days before experiments, cells were seeded in 60-mm plates in DMEM/F12 with charcoal stripped phenol red-free serum supplemented with nonessential amino acids and penicillin/streptomycin, at a density such that 75% to 90% confluence was obtained by the time of harvesting. Medium was refreshed 24 hours before exposure to 1 nmol/L E2, 10 nmol/L E2, or DMSO (control) for 16 hours. Subsequently, cells were harvested and RNA was extracted using the NucleoSpin RNA Isolation Kit (Macherey-Nagel). Details on measurement of E2 concentrations in eight fresh-frozen ER-positive breast tissue samples supplied by Erasmus Medical Center (Rotterdam, the Netherlands) are given in Supplementary Material and in references (17, 18).

Affymetrix HGU133Plus2.0 GeneChip microarray analysis was performed on the samples of the cell line experiments as described before (15, 19). See Supplementary Methods for more details. Data were deposited in GEO (GSE127760). All Affymetrix (CEL) data (generated in this study or retrieved from GEO) were processed using R (https://www.R-project.org) BioConductor affy and frma (20, 21) packages.

Datasets GSE10281 (22) and GSE20181 (14, 23) contain data from a cohort of postmenopausal patients with large strongly ER-positive (Allred score 5 or more) breast cancers recruited at the Edinburgh Breast Unit. Patients underwent neoadjuvant treatment with 2.5 mg daily letrozole for three months. Both sets contain data from fresh-frozen biopsies with at least 20% tumor content. The Affymetrix HGU133Plus2.0 GSE10281 (AffyP2 Edinburgh) dataset concerns 18 patients with samples collected at baseline and after three months of treatment. The Affymetrix HGU133A GSE20181 (AffyA Edinburgh) dataset contains gene expression data from biopsies of 62 patients collected at baseline and after approximately two weeks and three months of treatment. Clinical response was defined as tumor volume reduction of at least 50%, assessed after three months using three-dimensional ultrasound (14). Data from six patients were included in both datasets (overlap).

Dataset GSE48905 (24) contains data from a multi-center phase II study, in which postmenopausal women with operable, locally advanced (T2, 3, 4b; N0-3; M0) ER-positive breast tumors were randomized for neoadjuvant treatment with 500 mg or 250 mg fulvestrant for 16 weeks.

Tumor core biopsies (fresh-frozen tissue) were obtained before start of treatment and Affymetrix HGU133Plus2.0 microarrays were performed. Tumor volumes were measured by three-dimensional ultrasound and clinical response was classified as complete response (CR, disappearance of all lesions), partial response (PR, at least 65% reduction in tumor volume), progressive disease (PD, at least 73% increase in tumor volume), or stable disease (SD, none of the previous).

The Dutch TEAM IIa trial is a completed neoadjuvant phase II prospective trial (25). Briefly, 102 patients with ER-positive (≥50% nuclear ER staining) breast cancer were randomized to receive 3 or 6-month neoadjuvant exemestane treatment. Because of slow accrual, the study changed to single arm design consisting of 6-month therapy. Standard clinicopathologic characteristics, including PR and HER2 status, were assessed. Pretreatment cancer tissue biopsies and post-neoadjuvant treatment tumor resection specimens were collected and analyzed retrospectively; all tissue samples were formalin-fixed paraffin-embedded (FFPE). A total of 49 biopsy and 48 resection samples were eligible for analysis, of which 28 samples were paired cases from the same patient (Supplementary Fig. S1). Primary clinical outcome was reduction in tumor size using RECIST 1.1 criteria assessed by manual palpation.

Cancer tissue was separated from surrounding normal and stromal tissue using laser capture microdissection (LCM) as described by Espina and colleagues (26). See Supplementary Methods for more details. In FFPE biopsy samples, low tumor content limited the amount of available mRNA.

The method used to develop the ER pathway model has been described before (15). The model uses mRNA expression levels of ER target genes measured in a tissue sample to infer the odds in favor of a transcriptionally active ER transcription factor and, as a consequence, odds in favor of an active ER pathway. A Bayesian network representing the ER pathway transcriptional program (Fig. 1; ref. 15) describes (i) how target gene regulation depends on ER transcription complex activity and (ii) how expression level intensities in turn depend on regulation of the respective target genes. The network consists of three types of nodes: (i) ER transcription complex activation node; (ii) target gene regulation node, with states “down” and “up”; and (iii) expression intensity level nodes, with states “low” and “high,” each corresponding to an ER target gene. The previously described ER pathway model was developed for analysis of Affymetrix HGU133Plus2.0 (AffyP2) data, that is, target gene expression level nodes correspond to AffyP2 probe sets (15). During ongoing development toward a diagnostic ER pathway test to further improve the sensitivity and specificity of the previously described ER pathway model, an additional informative ER target gene (PDZK1) was included (27, 28). Following measurements of breast cancer tissue estrogen levels, the model was adjusted to estradiol levels in breast tissue by recalibration on data from MCF7 cultures exposed to 1 nmol/L estradiol (E2) from GSE35428 (Table 1; ref. 29), instead of 25 nmol/L E2 used to calibrate the original model (Table 1; ref. 30). Agreement between the original and optimized AffyP2 model was assessed on a collection of 5,395 breast cancer tissue and cell line samples (Supplementary Table S1). Sensitivity and specificity were assessed on 421 breast cancer cell line samples with well-established ER pathway activation state (Supplementary Table S1). Subsequently, the sensitivity-optimized (1 nmol/L) ER pathway model was adapted for use on two other mRNA measurement platforms, that is, Affymetrix HGU133A (AffyA) and 1-step qRT-PCR.

Adaptation of the ER pathway Bayesian network model between AffyP2 and AffyA microarray platforms consisted of selecting probe sets that were identical on both platforms and recalibrating on samples measured on the respective platforms. For lack of probe availability in the AffyA platform, nine probe sets and the CDH26 gene were removed from the AffyA computational model. The resulting AffyA model retained 50 probe sets, representing 27 of the 28 ER target genes. The model was calibrated using Affymetrix HGU133A mRNA expression data from dataset GSE9936 (Table 1; ref. 31). For the qRT-PCR platform, the 28-target gene ER pathway Bayesian network was reduced to a smaller network containing the most informative target genes from the AffyP2 model (15). The main reason for reducing the number of target genes in the qRT-PCR–based model was intended use on FFPE samples. To avoid multiplexing of qRT-PCRs leading to bad qRT-PCR performance, multiple separate diagnostic grade qRT-PCRs needed to be developed. Because frequently very small samples are available from biopsy FFPE material (as in the TEAM IIa study), the amount of RNA isolated from such samples is limiting with respect to the number of qRT-PCRs that can be performed. For this reason, a smaller selection of target genes was identified based on evidence of the gene being a direct ER target gene, as well as their capacity to differentiate between ER pathway inactive and active MCF7 cell cultures. Subsequently, this selection was validated by comparing inferred ER pathway activities obtained with the 28-gene Affymetrix model to inferred activities obtained with a reduced Affymetrix model adapted to this limited target gene set. Following this successful initial validation step, qRT-PCR assays were developed and validated according to standard procedures for each of the selected target genes and a set of reference genes. A Bayesian computational qRT-PCR ER pathway model (TEAMII PCR ER pathway model) was constructed as described previously (15) and calibrated using the breast cancer cell line calibration set (Table 1). Target gene mRNA expression values from these cell line experiments were obtained as described below and used as input to calibrate and validate this TEAMII PCR model. A research use prototype qRT-PCR plate, containing reagents for performing the qRT-PCR reactions, was used in this TEAM IIa study (Molecular Pathway Dx, Philips, TEAMII PCR prototype plate). To facilitate comparison between results obtained on different measurement platforms, a normalized ER pathway activity score (ERPAS) was computed for all models: The inferred odds in favor of an active ER pathway were transformed into a (base 2) logarithmic scale and normalized to scores ranging from 0 to 100, where 0 corresponds to the lowest and 100 corresponds to the highest odds in favor of an active ER pathway that a specific model can infer.

In summary, the new 1 nmol/L E2 model for AffyU133P2.0 contains the PDKZ1 gene in addition to the gene list described before (15). The 6 nmol/L E2 AffyA model lost one gene for which no probe sets were available on the AffyA platform. For the qRT-PCR–based model, best performing target genes identified using AffyU133P2.0 results were selected for further development of the model toward a commercially available diagnostic qRT-PCR–based ER pathway test (www.Philips.com/Oncosignal).

Quantitative expression levels of target genes and reference genes were determined using one-step qRT-PCR. See Supplementary Methods for more details. Normalized quantification cycle (C_q) values were used as input for the TEAMII PCR ER pathway model. For the TEAM IIa trial samples, inferred ERPASs were calculated in a blinded manner at Philips and returned to Leiden University Medical Center (Leiden, the Netherlands) for correlation to therapy response.

Targeted sequencing was used to perform mutation analysis of ESR1 and PIK3CA genes on 28 available TEAM IIa resection samples with either a paired biopsy sample or an ERPAS on the resection sample. DNA was extracted from whole samples and, if possible, from macrodissected samples. See Supplementary Methods and ref. 32.

Agreement between ER pathway models was evaluated by correlating activity scores obtained by the different models on identical sample sets (Table 1). Linear relationship between inferred activity scores was assessed by orthogonal regression, assuming equal residual variance. Agreement was assessed by Pearson correlation and mean deviations (MD) from the fitted orthogonal regression line were used as indication of the measurement error between two models.

Statistical analysis was performed using R (https://www.R-project.org) and SPSS 23.0 (IBM). Paired t tests were used to assess mean decrease in ERPAS between baseline and posttreatment measurements. χ² tests were used to assess association between outcome categories and dichotomized ERPAS. Two sample t tests and ANOVA were used to assess ERPAS and decrease in ERPAS association with clinical outcome and baseline parameters. All tests assumed unequal variance. Unless indicated otherwise, two-sided P values and 95% confidence intervals (CI) are reported.

Optimization of the previously described ER pathway model included adjusting sensitivity to match estradiol (E2) levels present in breast cancer tissue. For this, E2 levels were measured in tissue samples of ER-positive breast cancers. Measured concentrations varied between 2.2 and 7 nmol/L E2 (Supplementary Table S2), in line with concentrations reported by others (33, 34). The new AffyP2 ER pathway model was calibrated on MCF7 breast cancer cells stimulated with 1 nmol/L estradiol, which lies between E2 concentrations in normal breast (<1 nmol/L) and in cancer tissue (7 nmol/L in this study). Results obtained with this 1 nmol/L–calibrated ER model strongly correlated with results obtained with the previously described 25 nmol/L–calibrated ER model (15), corr = 0.94, and showed the expected increase in sensitivity (from 83% to 91%), without loss in specificity (92%; Supplementary Figs. S2 and S3). The 1 nmol/L model was used for analysis of AffyP2 datasets in this study. Because the ER pathway models are measurement platform specific, the AffyP2 model was adapted for use on AffyA and qRT-PCR data, as described in the Materials and Methods section, and compared with the 1 nmol/L AffyP2 reference model on sets of technical replicate samples (Table 1; ref. 35). Results obtained with the AffyA model strongly correlated with results obtained with the 1 nmol/L E2 AffyP2 model on the dataset with 16 × 2 technical duplicates (corr = 0.80, Supplementary Fig. S4A). Differences in sample preparation protocols reduced the correlation (corr = 0.63, Supplementary Fig. S4B), illustrating the importance of standardized sample preparation procedures (36, 37). The TEAMII PCR model correlated well with the AffyP2 model on 11 technical duplicates (corr = 0.83, Supplementary Fig. S5). MDs from the fitted regression lines give an indication of the measurement error between two methods and were similar when comparing the AffyP2 1 nmol/L model versus the AffyA and TEAMII PCR models and the two AffyP2 models (MD = 3.1, 5.6, and 2.3, respectively).

For the AffyP2 Edinburgh cohort, ERPASs decreased significantly upon letrozole treatment (41.6 vs. 28.4), mean decrease of 13.2 (CI = 9.7–16.8; paired t test P < 0.001; Fig. 2A). ERPAS at the start of the treatment was positively correlated to decrease in ERPAS after three months of therapy, indicating that the higher the baseline ER pathway activity, the higher the decrease in ER pathway activity upon treatment (corr = 0.88; Fig. 2B). Clinical response data were not available for this cohort.

Similarly, in the AffyA dataset of the Edinburgh cohort (Table 2; Fig. 3A and B), ERPAS decreased significantly between baseline and two weeks (46.9 vs. 39, mean decrease of 7.9, paired t test P < 0.001) and between baseline and three months of letrozole (46.9 vs. 39.5, mean decrease of 7.5, paired t test P < 0.001). No significant difference in ERPAS was seen between two weeks and three months of treatment, indicating an early and maintained response to letrozole for the patient group as a whole. Again, baseline ERPAS and decrease in ERPAS were positively correlated, both after two weeks and three months of treatment (corr = 0.74 and 0.83, respectively, Fig. 3C and D). In this cohort, clinical response data were available and both baseline ERPAS and decrease in ERPAS after two weeks of letrozole treatment were significantly higher in clinical responders to letrozole than in nonresponders (P = 0.02, Table 2; Fig. 3B and E). Upon letrozole treatment, average ERPAS decreased in responders to the same low activity score level of nonresponders and remained at such level up to the 3-month measurement point (Fig. 3A, B, and F).

Using a threshold for ER pathway activity, as described before (15), enables analysis at individual patient level with respect to clinical response to letrozole. To evaluate the predictive value of the AffyA ER pathway activity model in an unbiased way, a threshold of 42.7 was set prior to analysis. This value corresponds to an odds of 1:1, the point where the inferred probability of the ER pathway being active is 50%. This preset value is very close to 41.5, the best cut-off point obtained by ROC analysis of the present dataset. Using the predefined 42.7 threshold, 69% of patients (n = 38) had an active pretreatment ER pathway. This proportion decreased to 24% (n = 13) after two weeks and slightly increased to 31% (n = 17) after three months of treatment (Table 2). At baseline, 77% (n = 34) of clinical responders had an active ER pathway, against only 36% (n = 4) of nonresponders (Table 2). Accordingly, a baseline ER pathway active state was significantly associated with response to therapy (χ² P = 0.009; Table 2). Univariate logistic regression analysis reported a protective OR of 0.17 (CI = 0.04–0.69; P = 0.01) on the dichotomous active/inactive ER pathway activity status and a corresponding protective OR of 0.35 (CI = 0.13–0.95; P = 0.02) for an increase in 10 points in the ERPAS. As reference, relation of response to established markers (baseline tumor grade and ER Allred score) was also analyzed (Supplementary Table S3). Although tumor grade was associated with response to therapy (χ² test P = 0.03), ER Allred score was not (χ² P = 0.2). Sample size was small for reliable multivariate analysis, but initial analysis indicated that the ERPAS remained significant after adjustment for grade and Allred score (Supplementary Table S4).

Analysis of this cohort not only illustrated the expected decrease in ER pathway activity upon effective letrozole treatment, but also that roughly one third of ER-positive patients had a functionally inactive ER pathway, associated with nonresponder status. Finally, in the first two weeks of letrozole treatment, ERPAS decreased in the majority of patients (85%, n = 47), while between two weeks and three months, ERPAS slightly increased in the majority of patients (60%, n = 33; Table 2). In line, the observed decrease in ERPAS from baseline to two weeks was higher in responders than nonresponders (9.3 vs. 2.7; P = 0.02), and remained higher in responders after three months of treatment (8.6 vs. 2.9, P = 0.13).

For the TEAM IIa cohort (Table 3), the average ERPAS was 42.2 at baseline (n = 49) and decreased to 33.9 (n = 48) after exemestane therapy. In the 28 paired samples (before and after therapy), mean decrease in ERPAS was 9.1 (paired t test P < 0.001). In this cohort, again a positive but weak correlation was found between baseline ERPAS and decrease in ERPAS during treatment (corr = 0.29; Supplementary Fig. S6).

To test the value of the TEAMII PCR ER pathway model to predict response to aromatase inhibitors (AI), activity scores were correlated to decrease in tumor size assessed by manual palpation (Fig. 4D; Supplementary Table S5). Patients with PD during exemestane treatment had a significantly lower ERPAS at baseline (28.4, n = 3) compared with patients with non-PD (44.6, n = 35); t test P = 0.03 (Supplementary Table S5, response assessed at last measurement). A relation between pretreatment ERPAS and clinical response was already present at three months exemestane (non-PD = 44.6, n = 27 vs. PD = 28.7, n = 3; P = 0.04, Supplementary Table S5, response assessed at 3 months).

Palpation-based response assessment was available for 23 patients with corresponding paired biopsy/resection activity scores (Fig. 4A and C). During exemestane treatment, ERPAS decreased in all six patients with complete remission (CR). In two of seven patients with partial response (PR) and in three of eight patients with stable disease (SD), an increase in ERPAS was measured after neoadjuvant therapy, indicating failure of blocking ER pathway activity at least at the time of surgical resection. For the two patients with PD with available paired samples, ERPAS was low before and had not changed after treatment. Decrease in ERPAS was higher in non-PD (9.9, n = 21) than PD (−0.16, n = 2), paired t test P = 0.003 (Fig. 4C). Although ERPAS did not differ significantly at baseline for CR, PR, and SD patients, a clear distinction was observed between CR/PR and PD patients (Fig. 4A). In addition, decrease in ERPAS was larger in CR patients, becoming gradually smaller in PR and SD patients (Fig. 4D).

For this cohort, available clinical data and sample size was too limited for extensive analysis. However, brief examination of the association of response to established risk markers did not indicate any of the available markers were associated to PD, and ERPAS remained associated to PD after adjustment.

For the NEWEST cohort, Affymetrix data were only available from biopsies prior to start of neoadjuvant fulvestrant. Data were analyzed using the AffyP2 ER model (Fig. 4B). Baseline ERPAS was 28.7 in the 2 PD samples versus 36.6 in the 37 non-PD patients (difference of 8.0 CI = 2.2–13.7; t test P = 0.02), similar to results of the TEAM IIa study. Again, there was no significant difference in ERPAS between PR and SD patients.

Combined ER/PR positive IHC staining is generally thought to be a better indicator of response to endocrine treatment than ER status alone. Therefore, combined positive staining might be a more accurate predictor of an active ER pathway. Indeed, the PR positive biopsy samples of the TEAM IIa cohort (n = 36, 73%) had significantly higher ERPAS than the ER⁺/PR⁻ samples (48.7 vs. 34.3, two sample t test P = 0.003) and the decrease in ERPAS after treatment was larger in the ER⁺/PR⁺ group, although this did not reach significance (Supplementary Table S6). However, baseline PR protein staining alone did not correlate with response to AI therapy, neither in the TEAM IIa cohort (this study) nor in the Edinburgh cohort (reported in ref. 38).

ESR1 and PIK3CA mutations could be assessed in 15 samples (13 patients) of the TEAM IIa study (Supplementary Figs. S7 and S8); remaining samples were of insufficient quality. Although nonsynonymous ERS1 SNVs were only found in patients with partial response or stable disease (Supplementary Table S7), no literature evidence as to their functionality was found. A known PIK3CA (H1047R) activating mutation was detected in one CR and one SD patient.

We previously described the development of a new method to quantitatively measure activity of the ER pathway in a tissue sample using a computational model that interprets measured ER target gene mRNA levels (15). In the current study, its use to predict response to neoadjuvant endocrine treatment in patients with ER-positive breast cancer was evaluated, analyzing two clinical studies from Edinburgh (14, 22), the Dutch TEAM IIa cohort (25), and the NEWEST study (24). For this purpose, the ER pathway model was adjusted to match sensitivity to E2 concentrations in cancer tissue, and adapted to the different mRNA measurement platforms used in these clinical studies. After freezing the models, analysis of the independent patient cohorts was performed to provide per sample a quantitative ERPAS.

In all three clinical studies for which pre- and posttreatment data were available, endocrine treatment induced a decrease in ERPAS, and the decrease was already maximal after two weeks of treatment. In the studies for which clinical response data were available, the ERPAS prior to treatment correlated with clinical response. In the AffyA Edinburgh cohort, a high ERPAS was associated with favorable clinical response to letrozole. In the TEAM IIa and NEWEST studies, a low ERPAS was associated with progressive disease under exemestane and fulvestrant treatment, respectively. In general, patients with higher activity scores responded better to endocrine treatment than patients with lower activity scores, and treatment-induced decrease in ERPAS was largest in clinical responders and smallest in progressive disease patients. These observations are in accordance with broad evidence that AI and fulvestrant therapies are effective in patients with an active ER pathway driving tumor growth, by interfering with production of estradiol, and by blocking of ER transcriptional activity and degrading ER, respectively (9, 39). An inactive or minimally active ER pathway is expected to be associated with primary resistance and progressive disease under endocrine therapy, as we observed. Results are also in line with our earlier finding that patients with ER-positive breast cancer with low ERPASs, treated with adjuvant tamoxifen, have a shorter relapse-free survival (15).

In the TEAM IIa study, ERPAS was not significantly different among the CR, PR, and SD groups. A contributing factor may have been inaccuracy in assessment of tumor size by manual palpation, which was performed by different MDs at nonfixed time points during treatment. Nevertheless, a detailed comparison of tumor size at resection versus estimated size by palpation, mammography, US, and MRI in the TEAM IIa cohort suggested palpation as most accurate measurement method (40). In addition, in the NEWEST study, three-dimensional ultrasound assessment of tumor volume was used to assess clinical response, and again, a low ERPAS was associated with PD, while ERPAS could not statistically separate PR and SD. However, in both clinical studies, best responder patients (respectively CR and PR) were clearly distinguishable from PD patients by higher activity scores.

Positive ER/PR IHC staining is generally thought to be a more reliable indicator of response to endocrine therapy than positive ER staining alone. However, although PR staining did correlate ERPAS, it did not correlate with response to AI, indicating that combined ER/PR IHC alone is insufficient to decide on an active ER pathway and endocrine treatment, in agreement with clinical experience (41). This is likely due to the indirect relationship between ER transcriptional activity and expression level of the actual PR protein (42).

The Allred score for ER IHC staining has also been described as a predictive marker for neoadjuvant hormonal therapy response, with improved response at Allred scores higher than 5 (1). In the AffyA Edinburgh cohort for which this score was available, no relation with response to therapy was found, probably due to intentional preselection bias toward high Allred scores to enrich the responders group (38).

In a subgroup of patients from the AffyA Edinburgh and TEAM IIa cohorts (16% and 25%, respectively), ERPAS had increased at the end of AI treatment, indicating failure of blocking ER pathway activity, at least at the time of surgical resection or last biopsy. Possible explanations for measuring persistent ER pathway activity under AI are lack of compliance, stopping therapy prior to surgical resection, tumor heterogeneity, and emerging endocrine resistance. In both studies, patients were followed closely and compliance was estimated as high by the treating oncologists; however, noncompliance cannot be excluded with certainty. Tumor heterogeneity with respect to ER active and ER inactive cancer cell clones is another explanation, as shown recently in a preliminary study (43). A side study comparing activity scores in different areas of resected tumor from six TEAM IIa patients showed a variation in the ERPAS of 9 points or more in four patients (Supplementary Fig. S9). Resistance mechanisms such as ER activating mutations and activation of other tumor driving signaling pathways are well-known consequences of endocrine therapy, and the incidence of ER-activating mutations lies around 12% after treatment with aromatase inhibitors (44, 45). The observed increase in ERPAS between two weeks and three months treatment in the AffyA Edinburgh study, especially in patients with an initial reduction in activity score, is suggestive of ER-activating ESR1 mutations, unfortunately, no mutation information was available. In the TEAM IIa study, due to lack of sufficient tissue, ESR1 mutations could only be assessed in a few samples by targeted sequencing and only ESR1 SNVs with unknown functionality and unrelated to AI resistance or clinical response were identified. Sequencing of PIK3CA mutations revealed two known PIK3CA activating mutations, also without relation to response. These results emphasize the challenge of interpreting mutations with respect to functional impact, in line with recent clinical findings (46).

An important reason for successful adaptation of the ER pathway model to other mRNA measurement platforms than the original AffymetrixU133Plus2.0, is the Bayesian network model approach based on causal relationships between the ER transcription factor and well-established ER target genes (15, 16). The differences between this knowledge-based approach to develop signaling pathway tests, and data-driven methods, have been described and discussed before (15, 19, 47, 48). In brief, data-driven efforts to generate pathway gene expression signatures, such as gene set enrichment analysis and DAVID, have been typically based on discovery of noncausal associations between expressed genes and a pathway, by analyzing datasets with more or less specified relationships to the pathway (49, 50). Such methods are prone to finding spurious associations and carry a high risk of overfitting, which interferes with biological validation and performance on independent datasets and with adaptation to other gene measurement platforms (51). In contrast, the Bayesian ER pathway model was biologically validated, can calculate the ERPAS of an individual sample despite variations in the expressed target gene subset, and its performance is relatively independent of mRNA measurement platform used. The latter is illustrated by the highly comparable results obtained by analysis of the different patient cohorts, where ERPAS was measured using different mRNA measurement platforms. For example in the TEAM IIa and NEWEST cohorts, a low ERPAS predicted progressive disease under endocrine therapy, irrespective of being measured on fresh-frozen samples with Affymetrix microarrays, or on FFPE samples with qRT-PCR.

Despite strongly positive ER staining, some patients fail to respond to neoadjuvant endocrine therapy, potentially resulting in progressive disease. By measuring ER pathway activity, we were able to identify such patients based on their low ERPAS. By providing additional evidence that ER-positive breast cancers do not always have high ER pathway activity, we provide a rational explanation for the clinical issue of nonresponders. Analysis of three independent clinical breast cancer cohorts confirmed the potential value of measuring functional ER pathway activity to predict response to endocrine therapy. The Bayesian computational models for ER pathway activity, adapted to different measurement platforms, revealed comparable clinical results on fresh-frozen tissue as well as FFPE samples. The ERPAS has clinical value as a continuous measure: the higher the score, the more likely that the ER pathway was active in the analyzed sample and the more likely that the patient responded to antihormonal treatment. In a clinical sample consisting of a mixture of cells, an absolute ER pathway activity threshold is arbitrary, because the activity score represents averaged ER pathway activity in the cell mixture, thus depending on percentage of cancer cells and their level of ER pathway activity. In the future, a clinically useful approach might be the use of tertiles of activity scores, allowing middle tertile scores to indicate uncertain treatment response.

The frequently extremely small FFPE samples from pretreatment biopsies reflect the routine diagnostic situation and are a point of concern. Current developments are directed toward improvement of signaling pathway analysis on such samples. Further clinical validation of the ER pathway test is ongoing, as well as elucidating the role of other signaling pathways, such as the PI3K pathway, in hormonal resistance (19, 48).

J. Martens reports receiving commercial and other research support from Philips Research, Pamgene, Therawis, and Sanofi. J.R. Kroep reports receiving commercial research support from Philips Research. W. Verhaegh, A. van de Stolpe, and M.A. Inda are listed as inventors on patent WO 2013/011479 (Patent applicant: Koninklijke Philips). No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.A. Inda, P.J.K. Kuppen, A. Charehbili, J.W.M. Martens, A. van de Stolpe

Development of methodology: M.A. Inda, E.C. den Biezen-Timmermans, B. van der Burg, W. Verhaegh, A. van de Stolpe

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P.J.K. Kuppen, A. Charehbili, A. van Brussel, S.E. Fruytier, S. Kloet, B. van der Burg, A.K. Turnbull, J.M. Dixon, J.R. Kroep, C.J.H. van de Velde

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.A. Inda, E.J. Blok, A. Charehbili, E.C. den Biezen-Timmermans, S. Kloet, J.W.M. Martens, A.H. Sims, C.J.H. van de Velde, A. van de Stolpe

Writing, review, and/or revision of the manuscript: M.A. Inda, E.J. Blok, P.J.K. Kuppen, E.C. den Biezen-Timmermans, A. van Brussel, E.M.-K. Kranenbarg, J.W.M. Martens, A.H. Sims, A.K. Turnbull, J.M. Dixon, J.R. Kroep, C.J.H. van de Velde, A. van de Stolpe

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): A. Charehbili, E.M.-K. Kranenbarg, S. Kloet, A.K. Turnbull, J.R. Kroep

Study supervision: P.J.K. Kuppen, J.R. Kroep, C.J.H. van de Velde, A. van de Stolpe

Other (data related to steroid measurement with bioassays): B. van der Burg

Other (PI of the clinical trial): C.J.H. van de Velde

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.