Probabilistic graphical models relate immune status with response to neoadjuvant chemotherapy in breast cancer

Breast cancer is the most frequent tumor in women and its incidence is increasing. Neoadjuvant chemotherapy has become standard of care as a complement to surgery in locally advanced or poor-prognosis early stage disease. The achievement of a complete response to neoadjuvant chemotherapy correlates with prognosis but it is not possible to predict who will obtain an excellent response. The molecular analysis of the tumor offers a unique opportunity to unveil predictive factors. In this work, gene expression profiling in 279 tumor samples from patients receiving neoadjuvant chemotherapy was performed and probabilistic graphical models were used. This approach enables addressing biological and clinical questions from a Systems Biology perspective, allowing to deal with large gene expression data and their interactions. Tumors presenting complete response to neoadjuvant chemotherapy had a higher activity of immune related functions compared to resistant tumors. Similarly, samples from complete responders presented higher expression of lymphocyte cell lineage markers, immune-activating and immune-suppressive markers, which may correlate with tumor infiltration by lymphocytes (TILs). These results suggest that the patient’s immune system plays a key role in tumor response to neoadjuvant treatment. However, future studies with larger cohorts are necessary to validate these hypotheses.


INTRODUCTION
Breast cancer is the most common neoplasm and the fifth cause of cancer-associated death among women [1]. Estrogen receptor (ER), progesterone receptor (PR), and human epidermal growth factor receptor 2 (HER2) provide a system of classification and clinical diagnosis. Seventy percent of the tumors are hormonal receptor positive, and HER2 overexpression is observed in 15% of cases. ER+ and PR+ tumors respond to endocrine therapy, whereas tumors overexpressing HER2 respond to targeted therapies such as trastuzumab [2,3]. Tumors negative for

Research Paper: Immunology
Oncotarget 27587 www.oncotarget.com ER, PR and HER2 are known as Triple Negative Breast Cancer (TNBC) and do not respond to the aforementioned therapies.
A molecular classification of breast cancer defined four intrinsic subtypes [4]. Luminal A disease, which accounts for 67% of all tumors, shows high expression of genes related to hormone receptors and low expression of genes related to cell proliferation. Luminal B, HER2enriched and Basal-like subtypes have a more aggressive phenotype [5] [6,7].
Neoadjuvant chemotherapy has been increasingly administered to reduce the size of primary tumor, thus increasing the likelihood of breast conservation and enhancing survival [8]. Currently, there is no clinically useful molecular predictor of response to cytotoxic drugs in the neoadjuvant setting. Clinical parameters or the expression of single molecular markers (ie, Bcl-2, p53, MDR-1, and so on) show weak association with response and are not regimen-specific. Molecular subtyping may offer some help, as Luminal B and Basal-like tumors respond better than Luminal A tumors [9], but this is not accurate enough to make clinical decisions. As a consequence, many patients suffer the toxicity of useless neoadjuvant chemotherapy.
This study has been carried out using probabilistic graphical models, providing insights into the molecular biology of tumor response, allowing its use as a predictive model for response. These statistically inferred networks provide a deeper level of biological understanding in two main directions: giving support to previously identified biological observations, and giving new insights regarding novel biological interactions. Moreover, the transcriptional network approach has proven to be useful to unveil transcriptional regulation in breast cancer [10,11]. The objective of this study was to evaluate differences in gene expression patterns of breast cancer tumors from patients who had undergone neoadjuvant chemotherapy through a Systems Biology perspective.

Patient's characteristics
279 patients with histologically-confirmed primary non-metastatic breast adenocarcinoma from phase II trial (NCT00455533) [12] were included. They all had untreated tumors of at least 2 cm in size (T2-3, N0-3) regardless of hormone receptor or HER2 expression status. Clinical data were obtained from phase II trial (NCT00455533). Patient's clinical characteristics are provided in Table 1. On the basis of ER, PR and HER2 status, 111 tumors patients (39.78%) were ER+ or PR+ and Her2-(ER+ for now on), 28 (10.04%) were HER2+ and 140 (50.18%) were classified as TNBC. Patients received sequential neoadjuvant therapy starting with 4 cycles of doxorubicin/cyclophosphamide (AC), followed by 1:1 randomization to either ixabepilone or paclitaxel. All patients underwent definitive breast surgery 4 to 6

Breast cancer systems biology
Gene expression data from all tumor samples were used to build a probabilistic graphical model, with no other a priori information. The resulting graph was divided in eighteen branches (functional nodes) and a main function was assigned to each node by gene ontology analysis. The structure of the probabilistic graphical model clearly reflected different biological functions ( Figure  1) Functional node activities were then calculated as previously showed [10,15].

Functional structure of response to neoadjuvant chemotherapy
Patients were classified according to pathological response regardless of their tumor molecular subtype to study the response to neoadjuvant chemotherapy. Significant differences between functional node activities were observed in "Immune response (MHCII)" (node 9), "Immune response (B cell)" (node 11) and "Immune response (Interferon)" (node 12) nodes, in which, tumors attaining a complete response had higher activation ( Figure 2). Blown up pictures of the genes in the red boxes are provided in Supplementary Figures 1-5. A progressive decrease in the activity of immune functional nodes was seen depending on the response, being higher in tumors obtaining a CR and absent in those having a progression. Additionally, the relationship of immune nodes activities with the pathological response was evaluated using an ordinal logistic regression analysis. This analysis revealed that an increment of one unit in node 9, 11 and 12 activities increased the probability of a favorable response 1.739, 1.435 and 1.629 times respectively. By contrast, one unit increase in the activity of node 10 increased 0.519 times the probability of having an unfavorable response. On the other hand, PD tumors showed higher functional activity in "Cell cycle 1" (node 17) and "Cell cycle 2" (node18), followed by CR tumors.

Functional characterization of molecular subtypes
Patients in the network were further classified according to their molecular subtype (Basal-like, Luminal  Genes with an expression below 0 were represented in green; genes with an expression around 0 were represented in grey and genes with an expression above zero were represented in red. B. Functional node activities differences between pathological response groups: Box-and-whisker plots are Tukey boxplots. All p-values were two-sided and p < 0.05 was considered statistically significant. P-value < 0.05 (*); p-value < 0.01(**). A.U: arbitrary units.
In order to evaluate the functional implications between molecular subtypes and response to neoadjuvant therapy, data from patients of the same molecular subtype were mean centred and analysed independently. Luminal A group included no PD, whereas only one PD was found in Luminal B group, and was excluded from this analysis. Normal-like and HER2+ tumors were insufficient to perform subsequent analyses.
Concerning Luminal A and Luminal B subtypes, "Immune response (MHCII)" (node 9), "Immune response (chemotaxis)" (node 10), "Immune response (B cell)" (node 11) and "Immune response (Interferon)" Genes with an expression below 0 were represented in green; genes with an expression around 0 were represented in grey and genes with an expression above zero were represented in red. B. Functional node activities differences between molecular subtypes: Box-and-whisker plots are Tukey boxplots. All p-values were two-sided and p < 0.05 was considered statistically significant. P-value < 0.05 (*); p-value < 0.01(**); P-value< 0.001 (***); P-value <0.0001 (****). A.U: arbitrary units. www.oncotarget.com (node 12) functional nodes activities were higher in tumors attaining a CR, although differences were not statistically significant. As in the case of Basal-like tumors, a progressive decrease of activity in these nodes was observed from CR to SD.
Three separate probabilistic graphical models were built for Basal-like, Luminal-A and Luminal-B subtypes. As in the global network, tumors experiencing a CR had an increased activity of immune response-related nodes, although differences were not statistically significant.

Immunological markers
Markers of tumor-infiltrating lymphocytes were analysed to further characterize response according to the immune status. For that, patients were separated according to their pathological response and marker gene expression levels were compared between groups. Tumors obtaining a CR showed significantly higher expression levels of cell lineage markers (CD4, CD8 and CD20) as well as immune-activating (IGKC, CXCL9, CCL5, CXCL13) and immunosuppressive markers (IDO1, PD1) compared to the rest of tumors (Figure 4).
Oncotarget 27592 www.oncotarget.com immune functional nodes were related to pathological response to neoadjuvant chemotherapy. This correlation did not depend on the molecular subtype, indicating that a Systems Biology approach complements knowledge obtained from other research methods. Non-directed probabilistic graphical models allow managing large gene datasets and underscoring lots of gene interactions, many of which have not been previously described.
The activity of immune nodes was higher in tumors attaining a CR and decreased with the intensity of response. Tumors progressing on chemotherapy also showed increased activity in the nodes "Cell division 1" (node 17) and "Cell division 2" (node 18). These results suggest that the patient's immune system plays a crucial role in the response to neoadjuvant chemotherapy. Previous studies suggest that conventional therapies are effective in patients exhibiting some degree of immune activation in the tumor [16], supporting our findings. Chemotherapy may mediate the "immunogenic" death of tumor cells, thus facilitating an immune response against the disease [17].
As expected, all tumors attaining a CR-regardless of molecular subtype-showed significantly higher levels of cell lineage markers (CD4, CD8 and CD20) as well as immune-activating (IGKC, CXCL9, CCL5, CXCL13) and immunosuppressive markers (IDO1, PD1), suggesting a greater infiltration of immune cells High tumor-infiltrating lymphocytes (TILs) levels have been associated with increased CR rates in ER+ HER2+/-tumors [18] and also in TNBC [19]. However, high levels of PD-1+ TILs or Foxp3+ TILs have been related with poor prognosis [18]. Therefore, immune cell subpopulation profiles could help to predict response to neoadjuvant chemotherapy.
Basal-like and HER2+ subgroups have been associated with highest CR rates as opposed to Luminal and Normal-like tumors. However, the genes that were associated with CR in Basal-like subgroup were not associated with CR in the HER2+ subgroup, suggesting that different sets of genes are associated with CR in the different molecular subtypes [20,21]. In the present cohort, Basal-like tumors had the highest CR rate, as expected. However, the CR rate was poor in HER2+ tumors, probably because these patients did not receive anti-HER2 targeted therapy. On the other hand, BL1 tumors achieved a CR more commonly than other TNBC subtypes, as previously described [22]. Although node "Immune response (MHCII)" (node 9) had higher activity in Luminal A and Normal subtypes, the remaining nodes related to immune response showed increased activity in Basal-like tumors. This agrees with the fact that, in general, there are far fewer TILs in luminal disease than in HER2 or TNBC subtypes [23]. In fact, even though increased TILs concentrations are associated with increased frequency of response to neoadjuvant chemotherapy in all breast cancer subtypes, there is a different effect of TILs on survival in TNBC and luminal breast cancer. Increased TILs concentrations are associated with longer survival in TNBC and HER2+ disease, but not in luminal-HER2-negative tumours [24], suggesting again a differences in the biology of the immunological infiltrate across molecular subtypes.
One possible explanation of the higher "Immune response (MHCII)" (node 9) activity in Luminal A subtype could be the contribution of different immune cell types. Most types of immune cells were increased in TNBC compared with luminal-HER2-negative breast cancer. In TNBC, the presence of many immune cell subtypes, including B cells, T cells, and macrophages, were linked to improved survival [24]. By contrast, in luminal-HER2negative breast cancer, the presence of T cells was not prognostic for survival and the only cell types linked to improved prognosis were B and myeloid dendritic cells [24], which are MHCII presenting cells. Taking this into account, it would be interesting to perform further studies about MHCII presenting cells infiltration in luminal subtypes.
On the other hand, Basal-like tumors also had the highest activity in the node "Cell cycle" (nodes 17 and 18), which is in accordance with the fact proliferation renders tumor cells more sensitive to chemotherapy [6].
The neoadjuvant setting is appealing in the field of drug development because it allows early evaluation of efficacy. However, not all patients benefit from this approach, so markers predicting response to neoadjuvant chemotherapy are clearly necessary, as neoadjuvant therapy may have some drawbacks, such as promoting metastasis in some cases [25]. Our results suggest that immune activation in the tumor may identify responders. Although validation is needed, the use of these markers can help to determine the future use of neoadjuvant chemotherapy in breast cancer.

Patient's and samples origin and characteristics
A breast cancer tumor dataset, including gene expression and clinical data, was obtained from the Gene Expression Omnibus (GSE41998) and from a phase II trial (NCT00455533). Gene expression profiling was measured using an Affymetrix GeneChip, normalized and log2 transformed. Surgical specimens were evaluated by a pathologist at each study site. The pathological response was evaluated as the primary endpoint. A pathological complete response (CR) was defined by no histologic evidence of residual invasive adenocarcinoma in the breast and axillary lymph nodes, with or without the presence of ductal carcinoma in situ [12]. Responses were categorized as complete, partial, stabilization or progression.