In silico identification of drug target pathways in breast cancer subtypes using pathway cross-talk inhibition

Background Despite great development in genome and proteome high-throughput methods, treatment failure is a critical point in the management of most solid cancers, including breast cancer (BC). Multiple alternative mechanisms upon drug treatment are involved to offset therapeutic effects, eventually causing drug resistance or treatment failure. Methods Here, we optimized a computational method to discover novel drug target pathways in cancer subtypes using pathway cross-talk inhibition (PCI). The in silico method is based on the detection and quantification of the pathway cross-talk for distinct cancer subtypes. From a BC data set of The Cancer Genome Atlas, we have identified different networks of cross-talking pathways for different BC subtypes, validated using an independent BC dataset from Gene Expression Omnibus. Then, we predicted in silico the effects of new or approved drugs on different BC subtypes by silencing individual or combined subtype-derived pathways with the aim to find new potential drugs or more effective synergistic combinations of drugs. Results Overall, we identified a set of new potential drug target pathways for distinct BC subtypes on which therapeutic agents could synergically act showing antitumour effects and impacting on cross-talk inhibition. Conclusions We believe that in silico methods based on PCI could offer valuable approaches to identifying more tailored and effective treatments in particular in heterogeneous cancer diseases. Electronic supplementary material The online version of this article (10.1186/s12967-018-1535-2) contains supplementary material, which is available to authorized users.


Background
Breast cancer (BC), the most invasive cancer in women worldwide, is a heterogeneous disease, characterized by different subtypes that lead to different clinical prognosis and responses to treatments [1].
The advent of genome-wide technologies has made possible the generation of new hypotheses about the role of genomics in the efficiency of drugs developed for cancer and the event of adverse responses to cancer therapy.
In this context, several studies examined the effects of drugs considering protein network approaches [2].
In particular, the analysis of network models revealed that the partial inhibition of a small number of proteins belonging to a network has a higher impact on the disease than the complete inhibition of a single protein: indeed, drugs that have only one target (single hit) often do not affect complex system networks of proteins in an effective way [2]. Moreover, experimental studies have shown that cancer cells are able to resist to drug treatments by creating and establishing news interactions in order to have an alternative signaling [3]. In network approaches, the effect of a drug treatment on some proteins, represented by the nodes of a network, is amplified by the interactions of these proteins with other proteins in the networks, being these connections represented by edges [2]. However, notwithstanding useful to assess the drug effects on proteins, these approaches have not still impacted in increasing the number of efficient drugs or in suggesting new potential targets for cancer treatment.
Hence, two critical points for drug discovery tools are: (i) to inhibit not single but multiple targets at the same time, and (ii) to prevent the formation of new interactions that could lead to phenomena of resistance or inefficacy of the drug.
Moreover, these two critical conditions can be further complicated by the interactions between pathways (pathway cross-talk) that can modify the effect of drugs.
For instance, in some cancer cells, rapamycin-like drugs inhibit mTORC1 complex but at the same time indirectly activate phosphatidylinositol 3-kinase (PI3kinase) and AKT, mitigating the inhibitor effect of drugs [3]. Similarly, in triple negative BC the inhibition of AKT as a consequence of a drug is indirectly damped by the activation of receptor tyrosine kinases (RTKs), reducing the efficacy of the drug [4]. Moreover, pathway cross-talk, i.e. the one existing between EGFR and HER2, is a possible evasive way for the cell to develop resistance to HER family receptor inhibitors [5]. These examples illustrate how drugs targeting an individual protein or pathway may not yield to the expected therapeutic effect due to activation of alternative pathways that avoid the barrier inflicted by the cancer drug.
Therefore, understanding how and where the pathway cross-talk can be inhibited by drug treatments during a disease process ideally could lead to more effective therapies, reducing the problem of drug resistance.
Once a map of pathway cross-talk specific of a disease is known, a potential solution hindering the formation of alternative signaling pathways created in response to therapy could be the administration of drugs able to act on both direct and indirect targets (generated by the inhibition of direct targets).
However, on this potential solution, the current opinions are conflicting. Indeed, multi-target drugs show lower affinity than one-target-drugs [2,3]. One of the most promising approaches is the use of synergistic drug combinations therapy able to act on both direct and indirect targets [2,3]. The use of drug combinations could overcome drug resistance issues associated with high doses of single-hit drugs, and their best efficacy could lead to use a lower drug concentrations reducing unwanted side-effect toxicity [2,3].
Understanding the effect of individual or combined drugs is critical in clinical studies. An important issue is the low reliability of the cell lines to predict the efficacy of drugs since cell lines did not show to be good models [3]. Thus, computational methods are demanded to deepen the potential role of a drug in a context of pathway cross-talk.
A recent work by Jaeger et al. [6] has proposed a computational method to simulate pathway cross-talk inhibition (PCI) given by individual or combined drugs in BC. In that study the authors considered the pathway crosstalk between two different pathways for shared protein interactions. The authors developed their computational algorithm (PCI index) considering all those KEGG pathways that contain any of the primary targets of the Food and Drug Administration (FDA) approved drugs in BC. No selections on pathways have been performed.
In our recent studies on BC [7,8], we generated a computational approach to select a network of pathways specific for distinct BC subtypes and quantified their cross-talk. In these works, we focused on a network of pathways composed of interactions among ten pathways; our aim was to study the role of miRNAs regulating pathway interactions in distinct BC subtypes.
In the present study we propose a computational method based on the work of Jaeger et al. [6] and optimized by taking advantage from results of our previous studies [7,8]. More precisely, here we describe a procedure to build a network of pathways de-regulated for different BC subtypes using gene expression data from The Cancer Genome Atlas (TCGA) and a list of pathways obtained by Ingenuity Pathway Analysis (IPA). We quantified pathway cross-talk with a dissimilarity measure and we assessed potential drug target pathways through PCI. We then applied PCI to quantify the effects of individual or combined drugs in distinct BC subtypes on our network of pathways. Finally, we speculated about the mode of action of FDA approved drugs and new potential drug targets that could decrease the activity of pathway-cross talk and therefore enhance clinical efficacy.

Datasets
We applied the computational approach on four BC subtypes with different diagnostic classification and prognosis [9]: "luminal A" tumors expressing hormone receptors, with a favorable prognosis; "luminal B" tumors expressing hormone receptors and high expression of proliferation genes with a good prognosis although with an increased risk of recurrence; "basal-like" tumors lacking the expression of hormone receptors and HER2 but increased levels of cytokeratin (myoepithelial) (CK 5/6 and CK 17), with shorter observed survival; and "HER2" tumors overexpressing HER2, with the worse survival.
We considered gene expression data from tissue samples studied by IlluminaHiSeq RNASeqV2
All the steps were applied to each of the four BC subtypes. The first four steps were applied 50 times in order to obtain a solid network of de-regulated cross-talking pathways for each BC subtype. More specifically, in order to perform a bootstrapping, we implemented a Monte Carlo cross-validation by randomly selecting some portions of the dataset (60%) to build the training dataset of the classifier and the rest of data (40%) as testing dataset.
Step 1, 2, and 3 were applied on the training data set.
Step 4 was applied on both training and testing dataset. To avoid problems of unbalanced classes, we randomly generated the same number of samples for each class of BC subtypes and NS.
Then, for each BC subtype, we generated a networkbased model of BC subtype pathways (5. Pathway network), and we simulated in silico the drug-induced inhibition of pathways (6. Pathway cross-talk inhibition) focusing on the change of network efficiency.
At the end of the sixth steps, we found potential drug target pathways for each BC subtype, which, if inhibited, can change significantly the network efficiency. These drug targets could be considered for future applications in drug discovery (7. Drug target pathway network). Figure 1 summarizes the proposed computational approach.
We validated our network-based models of BC subtype pathways and BC drug target pathways.
For the first validation task, we assessed the accuracy of using our subtype pathways in classifying four different BC subtypes using the GEO dataset, as data set independent from the TCGA dataset used for the identification of the subtype-derived pathways. For the second validation task, we evaluated the mechanism of action of 13 FDA approved drugs for BC on our BC subtype drug target pathway network (DTPN). We compared the decrease of pathway cross-talk within the DTCN followed by inhibition of drug target pathways.

Step 1: differential expression analysis
For each BC subtype, differential expression analysis (DEA) was applied with respect to NS using TCGAbiolinks [12]. In particular, we used the edgeR package from Bioconductor [13] to find differential expressed genes (DEGs) between each BC subtype and NS. For each DEG we calculated the log fold-change between the two conditions and corrected p-values using Benjamini-Hochberg procedure for multiple testing correction [14]. We defined DEGs if the absolute value of log fold change was > 1 and p value < 0.01. This step was applied 50 times on the 50 training datasets producing every time a different list of DEGs.

Step 2: pathway enrichment analysis
For each BC subtype, we identified a group of pathways significantly enriched with the list of subtype-derived DEGs. The original list of pathways (589) was obtained from IPA. Pathway Enrichment Analysis (PEA) was performed with a Fisher's test between DEGs and genes within IPA pathways [15]. We defined pathways enriched with DEGs if p-value of Fisher's test was < 0.01. p-values were corrected using Benjamini-Hochberg procedure for multiple-testing correction [14]. As for the step 1, the  Step 3: pathway cross-talk Pathway cross-talk between pathways enriched with subtype-derived DEGs was quantified using a discriminating score (DS) [7,8]. This score was defined by comparing the mean of the gene expression levels of each pair of pathways enriched with subtype-derived DEGs: where M x and S x represent the mean and the standard deviation of the gene expression levels in the pathway x, and M y and S y represent the mean and the standard deviation of the gene expression levels in the pathway y.
From this step we created a matrix for each sample (BC subtypes and NS) containing a DS value for each pair of pathways enriched with subtype DEGs.

Step 4: classification
In order to select the best discriminating pairs of crosstalking pathways, we implemented a Random Forest Classification using the R-package [16,17] to classify each BC subtype versus NS using the DS matrix obtained from step 3 as input of the classifier.
For each combination of pathways we estimated the Area Under Curve (AUC) values by cross-validation method (k-fold cross-validation, k = 10). We used the following parameters: mtry (number of variables randomly sampled as candidates at each split) = sqrt(p), p being the number of variables in the matrix of data; ntree (number of trees grown) = 500.
The classification was performed on the training dataset 50 times for each matrix obtained from the previous three steps. Every time we obtained the top-10 pairs of subtype cross-talking pathways with the best AUC value. Then we validated these pathways on the testing dataset.
In conclusion, for each BC subtype and for all 50 bootstraps we obtained 10 × 50 couples of pathways with the best AUC values validated on the testing dataset.

Step 5: pathway network
For each BC subtype, by univocally selecting the pairs of cross-talking pathways better discriminating each subtype, we then generated a network-based model of BC subtype pathways, where the nodes of each subtype network represent pathways and the links that connect nodes the pathway cross-talks.
Step 6: pathways cross-talk inhibition In a network-based model, the network efficiency is defined as the sum of the reciprocals of the shortest direct path lengths between all pairs of network elements. If N is the number of network elements and d i,j is the shortest direct path lengths of two elements i and j, NE can range from 0 to 1, where 1 means that all nodes interact directly with each other expressing the best efficiency of the network [2].
For each BC subtype, the network efficiency (NE) of the disease was calculated starting from the subtype pathway network, the number of the pathways (nodes) and the shortest direct path lengths between two pathways (e.g. i, j in N).
The efficient drug-induced inhibition of a single pathway can be modelled by the elimination of all direct interactions at the pathway. The corresponding drug effect on the network can be measured by NE as an index of network integrity reduction measuring the drug efficiency [2].
For each BC subtype, we simulated, in silico, the drug-induced inhibition of pathways cross-talk by eliminating, one-by-one, all direct cross-talks. We thus quantified the new NE value of the network, that we called nNE. nNE is thus a function of k, being k the pathway or the combination of the pairs of pathways inhibited in the network.
As effect of this operation, nNE was < NE or > NE, resulting in increasing or decreasing the drug efficiency, respectively.

Step 7: drug target pathway network
For each BC subtype, we selected those cross-talking pathways from the network, that, if inhibited, caused nNE < NE, thus building a potential DTPN.
Furthermore, from the equation: we quantified the percentage activity of the subtype pathway network that is inhibited, defined as PCI [6]. (1)

Validation of DTPN
From the independent GEO dataset GSE58212 (121 luminal A BC, 69 luminal B BC, 36 basal BC and 32 HER2-overexpressing) we considered the gene expression levels belonging to our BC subtype pathway networks. We tested their ability to classify each subtype. We then investigated the role of the top-10 genes that obtained the best classification performance in the perspective of DTPN.
We validated DTPN by assessing if FDA-approved BC drugs have the target pathways within DTPN, and evaluating how these drugs could reduce effectively the NE of the DTPN.
For each BC subtype, we evaluated the association of 13 FDA-approved BC drugs (from Matador database) and DTPN. For this purpose, we applied PEA with a Fisher's test between gene targets of the drug and genes within IPA pathways. We defined pathways enriched with drug target genes if p-value of Fisher's test was < 0.01.
For each BC subtype, we measured the effects of the considered drugs, when administered individually and in combination on the DTPN. We quantified this effect by PCI.

Luminal A
For the Luminal A we found a pathway network composed of 73 individual pathways and 157 pathway crosstalks (Fig. 3a). AUC for this pathway network was 0.93 in the TCGA training dataset and 0.88 in the TCGA testing dataset, respectively (Fig. 3b). The efficiency of the pathway network (NE) was 0.3218. nNE, calculated after the drug-induced in silico inhibition of each single pathway (among the 73) or each combination of pair of pathways (among the 2628) is shown in Fig. 3c and d, respectively. The inhibition of 34/73 individual pathways reduces NE (red line) of the pathway network with nNE values that ranges from 0.3043 to 0.3217 (Fig. 3c). The inhibition of 1388 combinations (among 2628) reduces the NE (red line) of the pathway network with values that range from 0.2927 to 0.3217 (Fig. 3d). Thus, the inhibition of couples of pathways of the cloud seems to reduce the efficiency more than individual pathway.
Additional file 1 shows nNE values in the luminal A network after inhibition of individual and pairs pathway(s). The top-10 pathways, which, if inhibited, have resulted in a better efficiency reduction of the network are: 'LXR/RXR Activation' , 'Extrinsic Prothrombin Activation Pathway' , 'Estrogen Receptor Signaling' , 'Human Embryonic Stem Cell Pluripotency' , 'Ethanol Degradation IV' , 'RAR Activation' , 'Fatty Acid oxidation' , 'Bladder Cancer Signaling' , 'Factors Promoting Cardiogenesis in Vertebrates' and 'Glioma Invasiveness Signaling' . 'LXR/RXR Activation' plays the most important role since its single inhibition reduces the network efficiency of 5.5% (PCI). The inhibition of both 'LXR/ RXR Activation Pathway' and 'Extrinsic and Prothrombin Activation Pathway' is the most effective reducing the efficiency of the network of 9% (PCI). In line with the use of synergistic drug combinations, drugs acting on both 'LXR/RXR Activation Pathway' and 'Extrinsic and Prothrombin Activation Pathway' could be more effective. DTPN   raloxifene, tamoxifen, toremifene and vinblastine) ( Table 1). In particular, the inhibition of the 'Estrogen Receptor Signaling' reduces the NE of 2.50%. This inhibition could be obtained using three FDA approved drug (fulvestrant, raloxifene and tamoxifen) that have a significant number of target genes belonging to this pathway.
The 'FXR/RXR Activation Pathway' if inhibited, reduces the NE of 0.75%, and it is target of three drugs (methotrexate, toremifene and vinblastine). The 'RAR Activation Pathway' , if inhibited, reduces the NE of 1.5% and it is target of two drugs (raloxifene and tamoxifen). The 'Regulation of Cellular Mechanics by Calpain Protease' pathway, if inhibited, reduces the NE of 0.6% and is target of two drugs (fulvestrant and gemcitabine).
Among the drug target pathways found in our DTPN, we found 5/34 pathways that were ranked as part of the top 10 pathways, playing the major role in the reduction of NE. These pathways are 'Estrogen Receptor' targeted by tamoxifen and raloxifene; 'Human Embryonic Stem Cell Pluripotency' targeted by tamoxifen; 'RAR Activation Pathway' targeted by raloxifene and tamoxifen; 'Bladder Cancer Signaling' targeted by fulvestrant; 'Factors Promoting Cardiogenesis in Vertebrates' targeted by tamoxifen. In particular, tamoxifen acting on 4 pathways reduces the NE of 8.59%, thus it is very effective.
FDA-approved drugs do not seem to act on 'LXR/ RXR Activation' , 'Extrinsic Prothrombin Activation'; 'Ethanol Degradation IV' , 'Fatty Acid oxidation' and 'Glioma Invasiveness Signaling' . In particular the first two pathways, according to our findings, could be potential dug targets for anticancer drug treatment.
Furthermore, considering FDA-approved drugs, our study confirms that tamoxifen is the drug with the best efficacy on the pathway network since it is able alone to inhibit four pathways ('Estrogen receptor' , 'Human Embryonic Stem Cell Pluripotency' , 'Factors Promoting Cardiogenesis in Vertebrates' and 'RAR activation').  (1); factors promoting cardiogenesis in vertebrates (2); human embryonic stem cell pluripotency (3); RAR activation (4) (

Luminal B
In luminal B we found a pathway network composed of 73 individual pathways and 129 pathway cross-talks (Fig. 4a). In the TCGA training dataset mean AUC was 0.98 and 0.96 in the TCGA testing dataset (Fig. 4b).
The efficiency of the pathway network was 0.3272. The inhibition of 27/73 individual pathways reduces the NE of the pathway network with values that range from 0.3123 to 0.3270 (Fig. 4c). The inhibition of 1199/2628 combinations of pairs of pathways reduces the NE with values that range from 0.2935 to 0.32720 (Fig. 4d). As for luminal A, also in luminal B the inhibition of pairs of pathways of the pathway network reduce the efficiency more than individual pathway. DTPN of luminal B was built considering the 27 pathways that inhibited reduced the NE. We obtained 4 BC drugs (among the 13 FDA approved) (docetaxel, fulvestrant, raloxifene, and tamoxifen) that interact with the DTPN of luminal B (Table 2).
In particular, the 'Estrogen-mediated S-phase Entry Pathway' is target of three approved drugs (fulvestrant, raloxifene and tamoxifen). FDA-approved drugs do not seem to act on 'Cell Cycle Control of Chromosomal Replication' , 'P2Y Purigenic Receptor Signaling Pathway' and 'Growth Hormone Signaling'; pathways that according our findings were the best potential drug targets. In particular, tamoxifen acting on three pathways reduces the NE of 3%.
Furthermore, considering FDA-approved drugs our study demonstrates that docetaxel is the drug with the best action on pathway network since is able to inhibit 'Epithelial Adherens Junction Signaling' .

HER2-overexpressing BC
We found a pathway network in HER2 BC composed of 100 individual pathways and 222 pathway cross-talks (Fig. 5a). Mean AUC was 0.98 and 0.91 in the TCGA training and testing dataset, respectively (Fig. 5b).
The efficiency of the pathway network in HER2 was 0.3492. Figure 5c and d show nNE in case of inhibition of individual pathways or each combination of pairs of pathways, respectively. The inhibition of 39/100 individual pathways reduces the NE of the pathway network with values that range from 0.3309 to 0.3492 (Fig. 5c). The inhibition of 2271/4950 combinations of pathways reduces the NE of the pathway network with values that ranges from 0.3189 to 0.3492 (Fig. 5d). It is confirmed that the inhibition of couples of pathways of the pathway network reduces the efficiency more than individual pathway.
Additional DTPN was built considering 39 pathways whose inhibition reduces the efficiency of the network. Starting from the 13 BC drugs, we obtained five drugs (fluoracil, docetaxel, fulvestrant, methotrexate, and tamoxifen) that interact with the HER 2 DTPN (Table 3).
In particular, the 'Salvage Pathways of Pyrimidine Deoxyribonucleotides' is target of two FDA-approved BC drugs (fluoracil, and methotrexate).  reduces the NE of 0.2% and it is target of fluoracil, and methotrexate. Among the drug target pathways found in our DTPN we found 1 pathway that was ranked as to be part of the top 10 pathways playing the major role in the reduction of NE ('Human Embryonic Stem Cell Pluripotency').
FDA-approved drugs do not seem to act on 'Role of BRCA1 in DNA Damage Response Pathway' and 'NAD biosynthesis II (from tryptophan)' , pathways that, according to our findings, were the best potential drug targets.
Furthermore, considering FDA-approved drugs, our study confirms that tamoxifen is the drug with the best action on pathway network since it is able to inhibit 'Human Embryonic Stem Cell Pluripotency' .

Basal BC
We found a pathway network in basal BC composed of 43 individual pathways and 74 pathway cross-talks (Fig. 6a). In the TCGA training dataset the mean AUC was 0.98 and 0.97 in the TCGA testing dataset (Fig. 6b).   16:154 The efficiency of the pathway network in basal BC was 0.3445. nNE in case of the single inhibition of the 43 individual pathways or each combination of pairs of pathways are shown in Fig. 6c and d, respectively. The inhibition of 19/43 of individual pathways reduces the NE of the pathway network with values that ranges from 0.2959 to 0.3443 (Fig. 6c). The inhibition of 416/903 combinations of pathways reduces the NE of the pathway network with values that ranges from 0.2148 to 0.3445 (Fig. 6d). The inhibition of couples of pathways is confirmed to reduce the efficiency more than individual pathway. DTPN of basal BC was built considering 19 pathways that reduced NE. We obtained 6 drugs (fluoracil, capecitabine, fulvestrant, methotrexate, raloxifene and tamoxifen) that interact with the DTPN in basal BC (Table 4).
In particular, the 'Salvage Pathways of Pyrimidine Deoxyribonucleotides' is target of two FDA-approved drugs (fluorouracil, and methotrexate) and ' Aryl Hydrocarbon Receptor Signaling' of 4 FDA-approved drugs (fluorouracil, fulvestrant, raloxifene, and tamoxifene). Among drug target pathways found in our DTPN we found 1 pathway of the top 5 pathways playing the major role in the reduction of NE (' Aryl Hydrocarbon Receptor Signaling'). In particular, fluorouracil acting on 2 pathways reduces the NE of 6.86%, thus it is very effective.
The considered FDA approved drugs for BC do not seem to act on 'EIF2 Signaling' and 'Mismatch Repair in Eukaryotes' , pathways that, according to our findings, could be interesting potential drug targets.

Validation of DTPN
To classify luminal A, we considered the expression levels of 2693 genes from the independent GEO dataset GSE58212 and belonging to our luminal A pathway network. We found PRC1 (AUC = 0. PLCB3, ERBB2, MYL5, and WNT3 belong to ' Axonal Guidance Signaling' pathway, whose inhibition can reduce the NE of HER2 of 0.8% (PCI). PPP1R10 belongs to 'Protein Kinase A Signaling' . Its inhibition can reduce the NE of 2.8% (PCI).
In basal BC we considered expression levels of 1196 genes from the GEO dataset belonging to 43 pathways of basal pathway network. We found RHOB (AUC = 0. RARA, ESR1, and CCNE1 belong to ' Aryl Hydrocarbon Receptor Signaling' pathway whose inhibition can reduce

DTPN pathways in BC subytpes
We have constructed a model of BC heterogeneity based on different subtype pathway networks. Table 5 shows de-regulated pathways in common in all subtypes, which can be considered as responsible for the initial stages of BC, and those de-regulated pathways present only in one subtype, specific of the behaviour of that subtype. We can observe in particular that: (i) 'LXR/RXR Activation Pathway' , that in our analyses emerged as a new potential drug target for luminal A, is specific for this subtype; (ii) 'NAD biosynthesis II (from tryptophan)' pathway, that according to our findings was found the best potential drug target in HER2, is specific for HER2 subtype; and (iii) ' Aryl Hydrocarbon Receptor Signaling' , that is targeted by 4 FDA-approved drugs (fluorouracil, fulvestrant, raloxifene, and tamoxifene) is specific of basal subtype.

Discussion
The formation of new alternative signaling pathways upon drug treatment is one of the main causes of inefficacy or development of drug resistance in cancer. The study of alternative signalling through pathway crosstalk can help drug-development strategies. Therefore, in this work, we optimized a computational method to investigate the role of pathways target of FDA approved drugs or new drugs that specifically addresses this issue. Our in silico method is based on the identification and quantification of the pathway cross-talk for distinct cancer subtypes. The effects of drugs inhibiting individual or combined pathways can be simulated and measured with the purpose to find new potential drug-targets or synergistic combination drugs.
Most works that study drug-development approaches use network models. Usually the most used methods to examine the effects of drugs involve protein-protein interactions, metabolic control analysis or neural network [18][19][20][21]. Nevertheless, most of these methods are focused on the detection of drugs that have only one target (single hit) and often are not able to reveal the pathway interactions in an effective way [2]. Furthermore, these methods show to depend on a large number of parameters [2]. Pathway cross-talk can be an ideal framework to assess pathway interactions that could originate phenomena of resistance or inefficacy of the drug. The investigation of drug-induced modifications on the different levels of pathway cross-talk can improve the knowledge of drug effects and potential drug targets.
Although based by the computational method of Jaeger et al. [6], our application presents several differences: (1) the Jaeger's work was focused on cross-talk due to the presence of overlapping genes between different pathways, we instead explore cross-talk as regulatory interactions among distinct pathways; (2) in Jaeger's work all KEGG pathways were considered, we instead made a selection on pathways based on their different activity in different cancer subtype vs normal tissues. In line with the last scenario, the computational method of Jaeger et al. is suitable to model static pathway cross-talks, while our method could be used to represent dynamic disease progression.
Using Monte-Carlo cross validation and classification, we built networks of cross-talking pathways able to classify with high performance cancer subtypes versus NS. We inhibited pathway cross-talks in such subtypederived networks and identified a number of pathways that could be potential targets of drugs by measuring drug ability in reducing the network efficiency.
We applied our network-based model of cross-talking pathways to the BC subtype gene expression level of an independent GEO dataset, to test independently the ability of the pathway networks to classify BC subtypes. We found that the top 10 genes achieving the best classification performance belong to the expected cross-talking pathways.
We then associated some drugs already approved by FDA for BC with the networks of cross-talking pathways. The association was done considering the number of known drug-target proteins mapped on pathways and selecting potential synergistic combination of drugs. In general, the inhibition of pairs of pathways of the pathway network reduces the network efficiency more than the inhibition of individual pathways. The considered FDA approved drugs for BC act on several pathways included in the networks of cross-talking pathways. However, they do not act on other pathways that, according to our findings, could be interesting potential drug targets in BC.
In luminal A, we found that the individual inhibition of 'LXR/RXR Activation' pathway reduces the efficiency of the network of 5.6%. LXRs are nuclear receptors (NRs) involved in cholesterol, glucose, fatty acid metabolism and inflammatory responses. NRs are a family of transcription factors that bind to respond to lipophilic signaling molecules (ligand) regulating downstream effectors. They also represent one of the most important drugdevelopment targets, since the designing of synthetic compounds that mimic the functions of ligands can selectively modulate the activity of NRs [22]. Selective estrogen receptor modulators (SERMs) and aromatase inhibitors (AI) are examples of synthetic compounds, which block the production of estrogen [23]. LXRs ligands on prostate cancer reported an effect on cell proliferation and cell cycle, acting on p27 and SKP2 [24]. In BC the effects of LXRs ligands seem to be slightly different, as there is a decrease of SKP2 but not of p27 [25]. However, there are conflicting reports about the potential role of 'LXR/RXR Activation' pathway in BC [25][26][27][28]. In our diagnostic pathway cloud, synergistic drug combinations acting on: (1) 'LXR/RXR Activation' and 'Extrinsic and Prothrombin Activation' Pathways; (2) 'LXR/RXR Activation' and 'Estrogen Receptor Signaling'; (3) 'LXR/ RXR Activation pathway' and 'RAR Activation' reduce of almost 9% the efficiency of the network in luminal A. To date, synergistic drug combinations performing on these pathways are not reported. Therefore, in luminal A we suggest LXR/RXR Activation as drug-pathway target combined with a drug acting on 'Extrinsic and Prothrombin Activation' Pathway, 'Estrogen Receptor Signaling' or 'RAR Activation' .
In luminal B subtype we found that the individual inhibition of 'Cell Cycle Control of Chromosomal Replication' reduces the efficiency of the network of 4.5% while, in our pathway network, synergistic drug combinations acting on: (1) 'Growth Hormone Signaling' and 'Regulation of the Epithelial-Mesenchymal Transition' Pathways; (2) 'Mitotic Roles of Polo-Like Kinase' and 'Cell Cycle Control of Chromosomal Replication'; (3) 'Cell Cycle Control of Chromosomal Replication' and 'Growth Hormone Signaling' reduce of almost 10% the efficiency of the network. To date, synergistic drug combinations targeting on these pathways are not presented. Regarding FDA approved drugs, we confirm docetaxel and tamoxifen as promising candidates for synergistic drug combinations since they can modulate 'Epithelial Adherens Junction Signaling' and 'Regulation of the Epithelial-Mesenchymal Transition Pathway' , respectively. Previous studies showed that a combination therapy of these drugs in BC cell lines increases the anti-proliferative effects of single agents [29,30].
In HER2 subtype we found that the inhibition of BRCA1 in 'DNA Damage Response' reduces the efficiency of 5.2% while, in our pathway network, synergistic drug combinations acting on 'NAD biosynthesis II (from tryptophan)' and 'Role of BRCA1 in DNA Damage Response' reduce the efficiency of 8%.
In basal subtype we found that the inhibition of EIF2 Signaling reduces the efficiency of 8.5%. While, in our diagnostic pathway cloud, synergistic drug combinations acting on 'Mismatch Repair in Eukaryotes' and 'Cell Cycle:G2/M DNA Damage Checkpoint Regulation' reduce the efficiency of the network of 37%. Regarding FDA-approved drugs, we suggest fulvestrant, or fluorouracil and/or capecitabine as promising candidates for synergistic drug combinations since they can modulate 'Epithelial Salvage Pathways of Pyrimidine Deoxyribonucleotides' , ' Aryl Hydrocarbon Receptor Signaling' and 'Triacylglycerol Degradation' .
Our approach presents some limitations. One limit is the dependence on expression data analysis to identify affected pathway, mainly based on the use of changes in gene expression that are largest in size or level. However, some studies [31] show that even small changes in expression levels, which seems not to be of greatest functional significance, can be relevant in terms of phenotypic difference. For example a yeast cell containing a mutation in a gene that confers temperature sensitivity, thus essential for survival at non-permissive temperature, can indeed growth at the normal or permissive temperature. The temperature-sensitive gene in this mutant cell usually has a point mutation that leads to subtle changes in its protein production, but has a high impact in the protein function and thus in the phenotype of the cell. Thus, one possible risk of our approach is to lose significant features from a biological point of view (e.g. significant pathways), even if not from a statistical point of view. Furthermore, in our work, as in common in many studies, there is also a dependency in the data analysis workflow based on the used public datasets. Computational methods that use biological information (as pathways) from prior knowledge could have a bias towards pathways or genes that are better known, since they are more present in the literature and databases [32][33][34]. Moreover, different pathway databases can give different list of genes of the same pathway leading to affect the results of the computational method. On the other hand, public experimental data could derive by different experimental conditions and/ or lack of standardizations of experimental designs. In addition, the details of cell-specific information or clinical data of samples are not always available [32][33][34]. However, we optimized our computational method based on the reported databases by reducing the bias on data from such databases performing a Monte Carlo cross-validation.
As overall result, our approach have generated promising hypotheses for identifying altered pathways as potential therapeutic targets, either by using synergistic drug combinations and new/approved drugs. Despite the number of drugs is increasing, the specific molecular mechanisms underlying drug combinations with therapeutic effect remain often unclear. We believe that the current efforts based on pathway cross-talk drug strategies will provide key information that are required to decipher molecular mechanism contributing to resistance or inefficacy of drugs.